Skip to content

Commit

Permalink
Merge pull request #189 from EcoExtreML/fix_issue_181
Browse files Browse the repository at this point in the history 
Fix issue 181
  • Loading branch information
@SarahAlidoost
SarahAlidoost authored Sep 1, 2023
2 parents 0e0af05 + 911928c commit 51a0056
Show file tree
Hide file tree
Showing 37 changed files with 1,114 additions and 877 deletions.
9 changes: 9 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,14 @@
# Unreleased

**Changed:**

- All `cond*_` functions are refcatored and moved to `+conductivity` folder in
[#189](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/189)

**Fixed:**

- issue [#181](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/181)


<a name="1.3.0"></a>
# [1.3.0](https://github.com/EcoExtreML/STEMMUS_SCOPE/releases/tag/1.3.0) - 22 Jun 2023
Expand Down
Binary file modified run_model_on_snellius/exe/STEMMUS_SCOPE
View file
Binary file not shown.
43 changes: 43 additions & 0 deletions src/+conductivity/+hydraulicConductivity/calculateDTheta.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
function dtheta = calculateDTheta(subRoutine, heat_term, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m)

% get soil constants
SoilConstants = io.getSoilConstants();

switch subRoutine
case 0
% for DTheta_LLh, heat_term = hh
% for DTheta_UUh, heat_term = hh + hh_frez
A = (-theta_r) / (abs(heat_term) * log(abs(SoilConstants.hd / SoilConstants.hm))) * (1 - (1 + abs(alpha * heat_term)^n)^(-m));
B = alpha * n * m * (theta_m - gamma_hh * theta_r);
C = ((1 + abs(alpha * heat_term)^n)^(-m - 1));
D = (abs(alpha * heat_term)^(n - 1));
dtheta = A - B * C * D;
if (hh + hh_frez) <= SoilConstants.hd
dtheta = 0;
end
case 1
% heat_term = hh or hh + hh_frez
A = (theta_m - theta_r) * alpha * n;
B = abs(alpha * heat_term)^(n - 1) * (-m);
C = (1 + abs(alpha * heat_term)^n)^(-m - 1);
dtheta = A * B * C;
case 2
A = theta_ll - theta_l;
B = hh + hh_frez - h - h_frez;
dtheta = A / B;
case 3
A = theta_uu - theta_u;
B = hh - h;
dtheta = A / B;
case 4
dtheta = theta_s / phi_s * (hh / phi_s)^(-1 * lamda - 1);
case 5
% this case differs with case 1 only regarding `theta_m` and `theta_s`
% see issue 181, item 7
% heat_term = hh or hh + hh_frez
A = (theta_s - theta_r) * alpha * n;
B = abs(alpha * heat_term)^(n - 1) * (-m);
C = (1 + abs(alpha * heat_term)^n)^(-m - 1);
dtheta = A * B * C;
end
end
65 changes: 65 additions & 0 deletions src/+conductivity/+hydraulicConductivity/calculateDTheta_UUh.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,65 @@
function dtheta_uuh = calculateDTheta_UUh(theta_uu, theta_m, theta_ll, gamma_hh, SoilVariables, VanGenuchten)
theta_s = VanGenuchten.Theta_s;
theta_r = VanGenuchten.Theta_r;
alpha = VanGenuchten.Alpha;
n = VanGenuchten.n;
m = VanGenuchten.m;

theta_u = SoilVariables.Theta_U;
theta_l = SoilVariables.Theta_L;
hh = SoilVariables.hh;
h = SoilVariables.h;
h_frez = SoilVariables.h_frez;
hh_frez = SoilVariables.hh_frez;
phi_s = SoilVariables.Phi_s;
lamda = SoilVariables.Lamda;

% get model settings
ModelSettings = io.getModelSettings();

heatTerm = hh + hh_frez;
if ModelSettings.SWCC == 1
if ModelSettings.SFCC == 1
if hh >= -1
if hh_frez >= 0
dtheta_uuh = 0;
else
if ModelSettings.Thmrlefc
subRoutine = 0;
dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
elseif abs(hh - h) < 1e-3
subRoutine = 1;
dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
else
subRoutine = 2;
dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
end
end
elseif ModelSettings.Thmrlefc
subRoutine = 0;
dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
end
elseif hh >= -1e-6 || hh <= -1e7
dtheta_uuh = 0;
elseif ModelSettings.Thmrlefc || abs(hh - h) < 1e-3
subRoutine = 1;
dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
else
subRoutine = 3;
dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
end
else
if hh >= phi_s || hh <= -1e7
dtheta_uuh = 0;
elseif ModelSettings.Thmrlefc
subRoutine = 4;
dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
elseif abs(hh - h) < 1e-3
subRoutine = 5;
dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
else
subRoutine = 3;
dtheta_uuh = conductivity.hydraulicConductivity.calculateDTheta(subRoutine, heatTerm, theta_s, theta_r, theta_m, gamma_hh, theta_ll, theta_l, theta_uu, theta_u, hh, h, hh_frez, h_frez, phi_s, lamda, alpha, n, m);
end
end
end
12 changes: 12 additions & 0 deletions src/+conductivity/+hydraulicConductivity/calculateGamma_hh.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
function gamma_hh = calculateGamma_hh(hh)
% get soil constants
SoilConstants = io.getSoilConstants();

if abs(hh) >= abs(SoilConstants.hd)
gamma_hh = 0;
elseif abs(hh) >= abs(SoilConstants.hm)
gamma_hh = log(abs(SoilConstants.hd) / abs(hh)) / log(abs(SoilConstants.hd) / abs(SoilConstants.hm));
else
gamma_hh = 1;
end
end
23 changes: 23 additions & 0 deletions src/+conductivity/+hydraulicConductivity/calculateKL_h.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
function kl_h = calculateKL_h(mu_w, se, Ks, m)

% load Constants
Constants = io.define_constants();

MU_WN = Constants.MU_W0 * exp(Constants.MU1 / (8.31441 * (20 + 133.3)));
CKT = MU_WN / mu_w;
if se == 0
kl_h = 0;
else
kl_h = CKT * Ks * (se^(0.5)) * (1 - (1 - se^(1 / m))^m)^2;
end
if kl_h <= 1E-20
kl_h = 1E-20;
end
if isnan(kl_h) == 1
kl_h = 0;
warning('\n case "isnan(kl_h) == 1", set "kl_h = 0" \r');
end
if ~isreal(kl_h)
warning('\n case "~isreal(kl_h)", dont know what to do! \r');
end
end
57 changes: 57 additions & 0 deletions src/+conductivity/+hydraulicConductivity/calculateSe.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,57 @@
function se = calculateSe(theta_ll, gamma_hh, SoilVariables)
% get model settings
ModelSettings = io.getModelSettings();

% get soil constants
SoilConstants = io.getSoilConstants();

hh_frez = SoilVariables.hh_frez;
POR = SoilVariables.POR;
hh = SoilVariables.hh;
phi_s = SoilVariables.Phi_s;
se = SoilVariables.Se;
seCalculated = theta_ll / POR;

if ModelSettings.SWCC == 1
if ModelSettings.SFCC == 1
if hh >= -1
if hh_frez >= 0
se = 1;
elseif Thmrlefc
if (hh + hh_frez) <= SoilConstants.hd
se = 0;
else
se = seCalculated;
end
end
elseif ModelSettings.Thmrlefc
if (hh + hh_frez) <= SoilConstants.hd
se = 0;
else
se = seCalculated;
end
end
elseif theta_ll <= 0.06 || hh <= -1e7
se = 0;
else
se = seCalculated;
end
else
if hh + hh_frez <= -1e7 || hh <= -1e7
se = 0;
elseif hh + hh_frez >= phi_s
se = 1;
else
se = seCalculated;
end
end

if se > 1
se = 1;
elseif se < 0
se = 0;
end
if isnan(se) == 1
warning('\n case "isnan(se) == 1" happens. Dont know what to do! \r');
end
end
25 changes: 25 additions & 0 deletions src/+conductivity/+hydraulicConductivity/calculateTheta.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
function theta = calculateTheta(subRoutine, theta_m, heat_term, gamma_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m)
% get soil constants
SoilConstants = io.getSoilConstants();
hd = SoilConstants.hd;

switch subRoutine
case 0
% heat_term = hh or hh + hh_frez
theta = gamma_hh * theta_r + (theta_m - gamma_hh * theta_r) / (1 + abs(alpha * heat_term)^n)^m;
if heat_term <= hd
theta = 0;
end
case 1
% heat_term = hh or hh + hh_frez
theta = theta_s * (heat_term / phi_s)^(-1 * lamda);
if heat_term <= -1e7
theta = Theta_r;
elseif heat_term >= phi_s
theta = theta_s;
end
case 2
% heat_term = hh or hh + hh_frez
theta = Theta_r + (theta_s - Theta_r) / (1 + abs(alpha * heat_term)^n)^m;
end
end
23 changes: 23 additions & 0 deletions src/+conductivity/+hydraulicConductivity/calculateTheta_II.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
function theta_ii = calculateTheta_II(tt, xcap, hh, Theta_II)

% get model settings
ModelSettings = io.getModelSettings();

Tf1 = 273.15 + 1;
Tf2 = 273.15 - 3;

theta_ii = Theta_II;

if hh <= -1e7
theta_ii = 0;
end
if ModelSettings.SWCC == 1 && ModelSettings.SFCC ~= 1
if tt + 273.15 > Tf1
theta_ii = 0;
elseif tt + 273.15 >= Tf2
theta_ii = 0.5 * (1 - sin(pi() * (tt + 273.15 - 0.5 * Tf1 - 0.5 * Tf2) / (Tf1 - Tf2))) * xcap;
else
theta_ii = xcap;
end
end
end
54 changes: 54 additions & 0 deletions src/+conductivity/+hydraulicConductivity/calculateTheta_LL.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,54 @@
function theta_ll = calculateTheta_LL(theta_uu, theta_ii, theta_m, gamma_hh, SoilVariables, VanGenuchten)
% get model settings
ModelSettings = io.getModelSettings();

% load Constants
Constants = io.define_constants();

hh = SoilVariables.hh;
hh_frez = SoilVariables.hh_frez;
phi_s = SoilVariables.Phi_s;
lamda = SoilVariables.Lamda;

theta_s = VanGenuchten.Theta_s;
theta_r = VanGenuchten.Theta_r;
alpha = VanGenuchten.Alpha;
n = VanGenuchten.n;
m = VanGenuchten.m;

heatTerm = hh + hh_frez;

% calculate theta_ll
theta_ll = SoilVariables.Theta_LL;
if ModelSettings.SWCC == 1
if ModelSettings.SFCC == 1
if hh >= -1 && hh_frez >= 0
theta_ll = theta_s;
elseif ModelSettings.Thmrlefc
subRoutine = 0;
theta_ll = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, heatTerm, gamma_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m);
end
else
if hh >= -1e-6 && TT + 273.15 > (273.15 + 1)
theta_ll = theta_s;
elseif hh <= -1e7
theta_ll = theta_r;
elseif TT + 273.15 > (273.15 + 1)
subRoutine = 2;
theta_ll = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, heatTerm, gamma_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m);
else
theta_ll = theta_uu - theta_ii * Constants.RHOI / Constants.RHOL;
end
if theta_ll <= 0.06
theta_ll = 0.06;
end
end
else
if hh >= phi_s
subRoutine = 1;
theta_ll = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, heatTerm, gamma_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m);
elseif hh <= -1e7
theta_ll = theta_r;
end
end
end
43 changes: 43 additions & 0 deletions src/+conductivity/+hydraulicConductivity/calculateTheta_UU.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
function theta_uu = calculateTheta_UU(theta_m, gamma_hh, SoilVariables, VanGenuchten)

hh = SoilVariables.hh;
phi_s = SoilVariables.Phi_s;
lamda = SoilVariables.Lamda;

theta_s = VanGenuchten.Theta_s;
theta_r = VanGenuchten.Theta_r;
alpha = VanGenuchten.Alpha;
n = VanGenuchten.n;
m = VanGenuchten.m;

% get model settings
ModelSettings = io.getModelSettings();

% calculate theta_uu
if ModelSettings.SWCC == 1
if ModelSettings.SFCC == 1
if hh >= -1
theta_uu = theta_s;
elseif ModelSettings.Thmrlefc
if gamma_hh == 0
theta_uu = 0;
else
subRoutine = 0;
theta_uu = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, hh, gamma_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m);
end
end
else
if hh >= -1e-6
theta_uu = theta_s;
elseif hh <= -1e7
theta_uu = theta_r;
else
subRoutine = 2;
theta_uu = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, hh, gamma_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m);
end
end
else
subRoutine = 1;
theta_uu = conductivity.hydraulicConductivity.calculateTheta(subRoutine, theta_m, hh, gamma_hh, theta_s, theta_r, lamda, phi_s, alpha, n, m);
end
end
15 changes: 15 additions & 0 deletions src/+conductivity/+hydraulicConductivity/calculateTheta_m.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,15 @@
function theta_m = calculateTheta_m(gamma_hh, VanGenuchten, POR)

theta_s = VanGenuchten.Theta_s;
theta_r = VanGenuchten.Theta_r;
alpha = VanGenuchten.Alpha;
n = VanGenuchten.n;
m = VanGenuchten.m;

theta_m = gamma_hh * theta_r + (theta_s - gamma_hh * theta_r) * (1 + abs(alpha * (-1))^n)^m;
if theta_m >= POR
theta_m = POR;
elseif theta_m <= theta_s
theta_m = theta_s;
end
end
Loading

0 comments on commit 51a0056

Please sign in to comment.