From 2b0bb13fe0602fccbed594d51d7e1a4724c1b94a Mon Sep 17 00:00:00 2001 From: Peter La Follette Date: Mon, 25 Mar 2024 14:15:50 -0700 Subject: [PATCH] tested LASAM over approximately 40k parameter sets. Resolved mass balance errors and infinite loop errors that became apparent when LASAM was run with many parameter sets. Also updated field capacity formulation to use a head rather than a relative water content, and updated caluclation for capillary drive G to use an adaptive integral, both as in LGARTO. Cleaned up some commenting and redundancy after edits were made. --- src/aet.cxx | 3 -- src/lgar.cxx | 30 +++++---------- src/soil_funcs.cxx | 91 +--------------------------------------------- 3 files changed, 10 insertions(+), 114 deletions(-) diff --git a/src/aet.cxx b/src/aet.cxx index 1df0688..60cb8ab 100644 --- a/src/aet.cxx +++ b/src/aet.cxx @@ -28,8 +28,6 @@ extern double calc_aet(double PET_timestep_cm, double time_step_h, double wiltin struct wetting_front *current; double theta_wp; - - // double relative_moisture_at_which_PET_equals_AET = 0.75; double Se,theta_e,theta_r; double vg_a, vg_m, vg_n; @@ -46,7 +44,6 @@ extern double calc_aet(double PET_timestep_cm, double time_step_h, double wiltin vg_n = soil_properties[soil_num].vg_n; // compute theta field capacity - // double theta_fc = (theta_e - theta_r) * relative_moisture_at_which_PET_equals_AET + theta_r; double head_at_which_PET_equals_AET_cm = 340.9;//*10/33; //340.9 is 0.33 atm, expressed in water depth, which is a good field capacity for most soils. //Coarser soils like sand will have a field capacity of 0.1 atm or so. double theta_fc = calc_theta_from_h(head_at_which_PET_equals_AET_cm, vg_a,vg_m, vg_n, theta_e, theta_r); diff --git a/src/lgar.cxx b/src/lgar.cxx index ceea916..c402e79 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -1154,7 +1154,7 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf the layer depth */ current->depth_cm = fmin(current->depth_cm, column_depth); - double theta_old = current->theta; + double theta_old = current->theta; //might not be necessary if (current->dzdt_cm_per_h == 0.0 && current->to_bottom == FALSE) // a new front was just created, so don't update it. current->theta = current->theta; else @@ -1337,11 +1337,10 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf // loop to adjust the depth for mass balance int iter = 0; - bool is_WF_deeper_than_domain = FALSE; while (fabs(mass_balance_error - tolerance) > 1.E-10) { iter++; if (iter>10000000) { - printf("mass balance closure not possible within 10000000 iterations. case 1. Timeout \n"); + printf("mass balance closure not possible within 10000000 iterations. Timeout \n"); // printf("depth_new: %lf \n", depth_new); printf("states: \n"); listPrint(*head); @@ -1366,12 +1365,6 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf current_mass = lgar_calc_mass_bal(cum_layer_thickness_cm, *head); mass_balance_error = fabs(current_mass - mass_timestep); - if (depth_new > cum_layer_thickness_cm[num_layers]){ - is_WF_deeper_than_domain = TRUE; - // *volin_cm = *volin_cm - mass_balance_error; - // break; - } - } } @@ -1401,9 +1394,10 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf all wetting fronts, and then a new loop does lower boundary crossing for all wetting fronts. this is a thorough way to deal with these scenarios. */ - //As of 20 March 2024, we approach wetting front merging / layer boundary crossing / lower boundary crossing in a new way. + //As of 20 March 2024, we have added a new method by which wetting fronts merge or cross a layer or the model lower boundary, although it is commented out for now. //we now have a while loop that detects if these scenarios are necessary and then only does them if so, instead of relying on - //a fixed order of merging and crossing (which was previously merge->cross layers->merge->cross lower bdy). + //a fixed order of merging and crossing (which was previously merge->cross layers->merge->cross lower bdy), it can automatically be detected which mergeing / crossing events are necessary. + //however, thius code slowed down some simulations by orders of magnitude. Since the old version of relying on a fixed, predetermined order of merging and crossing seems to work, I'd say don't fix it. /* Note: we check for dry over wet case before we call merge_wetting_fronts to avoid negative wetting fronts due to unknown corner/rare cases */ @@ -1509,7 +1503,7 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf // make sure all psi values are updated // there is a rare error where, after a wetting front crosses a layer boundary to a soil layer with a much more sensitive soil water retention curve, and the wetting front is near but not at saturation, // the barely unsatruated wetting front in the new layer will update its psi value as 0. This causes unequal psi values among adjacent wetting fronts in different layers and then a mass balance error. - // For example, consider the following: (1.0/(pow(1.0+pow(0.001039*0.00296620706633,2.938410),0.659680))*(0.618625-0.052623)+0.052623). This expression, using specific values in calc_theta_from_h, will yield a theta of theta_e (0.618625), while having a psi value that is slightly greater than 0, which is 0.00296620706633. + // For example, consider the following: (1.0/(pow(1.0+pow(0.001039*0.00296620706633,2.938410),0.659680))*(0.618625-0.052623)+0.052623). This expression, using specific values in calc_theta_from_h, will yield a theta of theta_e (0.618625), while having a psi value that is slightly greater than 0, which is 0.00296620706633 in this case. // While that is ok for this soil layer in particular, adjacent wetting fronts above this one with a less sensitive soil water retention curve will yield a non-theta_e value for the psi value that is slightly above 0. // A solution is to either not run the following code, or to not run it when the wetting front is very close to saturation with a very sensitive soil water retention curve. Adding code that runs the following only for psi>1. @@ -2076,7 +2070,7 @@ extern double lgar_insert_water(bool use_closed_form_G, int nint, double timeste // if ( (current_mass)/max_storage > 0.99 ){ // printf("warning: vadose zone is 99 percent full or greater. If you are using the model in an environment with more precipitation than PET, LGAR is not an appropriate model because its lower boundary condition is no flow. \n "); - // } //turning off for now + // } //turning off for now; should really print like once per model run or so, not every time step double ponded_depth_temp = *ponded_depth_cm; @@ -2115,12 +2109,6 @@ extern double lgar_insert_water(bool use_closed_form_G, int nint, double timeste // order is important here; assign zero to ponded depth once we compute volume in and runoff *volin_this_timestep = fmin(*ponded_depth_cm, fp_cm); // runoff = *ponded_depth_cm < fp_cm ? 0.0 : (*ponded_depth_cm - *volin_this_timestep); - // if (head->psi_cm<1e-8 && head->depth_cm==(cum_layer_thickness_cm[-1])){// ok this works for 1 layer, but generally for multilayer you want to check if the storage is saturated - //replaced this code with code that reduces f_p after if's calculated if there is not enough room in the soil (currently line 1915) - // runoff = *ponded_depth_cm; - // *volin_this_timestep = 0; - // listPrint(head); - // } *ponded_depth_cm = 0.0; } @@ -2577,7 +2565,7 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm iter++; // if (iter>10000000) { - // printf("mass balance closure not possible within 10000000 iterations. case 2. Timeout \n"); + // printf("mass balance closure not possible within 10000000 iterations. Timeout \n"); // printf("psi_cm_loc: %lf \n", psi_cm_loc); // listPrint(*head); // abort(); @@ -2693,7 +2681,7 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm } - +//while this function may or may not be used in LASAM, it will be used in LGARTO extern int lgar_correction_type(int num_layers, double* cum_layer_thickness_cm, struct wetting_front* head){ int correction_type = 0; struct wetting_front *current = head; diff --git a/src/soil_funcs.cxx b/src/soil_funcs.cxx index b252282..ad50367 100755 --- a/src/soil_funcs.cxx +++ b/src/soil_funcs.cxx @@ -26,100 +26,11 @@ /* used for redistribution following the equation published by Ogden and Saghafian 1995. */ /* to compile: "gcc one_block.c -lm" */ /* */ -/* author: Fred Ogden, June, 2021, edited by Peter La Follette in 2023 */ +/* author: Fred Ogden, June, 2021, edited by Peter La Follette in 2023 for adaptive integral */ //updated to calculate G using an adaptive integral method, rather than always using a fixed number of trapezoids. //This really helps when the limits of integration in terms of psi are quite far apart, and also helps to save runtime. /***********************************************************************************************/ -// extern double calc_Geff(bool use_closed_form_G, double theta1, double theta2, double theta_e, double theta_r, -// double vg_alpha, double vg_n, double vg_m, double h_min, double Ksat, int nint, -// double lambda, double bc_psib_cm) -// { -// double Geff; // this is the result to be returned. - -// if (use_closed_form_G==false) { -// // local variables -// // note: units of h in cm. units of K in cm/s -// double h_i,h_f,Se_i,Se_f; // variables to store initial and final values -// double Se; -// double h2; // the head at the right-hand side of the trapezoid being integrated [m] -// double dh; // the delta h over which integration is performed [m] -// double Se1,Se2; // the scaled moisture content on left- and right-hand side of trapezoid -// double K1,K2; // the K(h) values on the left and right of the region dh integrated [m] - - -// Se_i = calc_Se_from_theta(theta1,theta_e,theta_r); // scaled initial water content (0-1) [-] -// Se_f = calc_Se_from_theta(theta2,theta_e,theta_r); // scaled final water content (0-1) [-] - -// h_i = calc_h_from_Se(Se_i,vg_alpha,vg_m,vg_n); // capillary head associated with Se_i [cm] -// h_f = calc_h_from_Se(Se_f,vg_alpha,vg_m,vg_n); // capillary head associated with Se_f [cm] - -// if(h_i < h_min) {/* if the lower limit of integration is less than h_min FIXME */ -// //return h_min; // commenting out as this is not used in the Python version -// } - -// if (verbosity.compare("high") == 0) { -// // debug statements to see if calc_Se_from_h function is working properly -// Se = calc_Se_from_h(h_i,vg_alpha,vg_m,vg_n); -// printf("Se_i = %8.6lf, Se_inverse = %8.6lf\n", Se_i, Se); - -// Se = calc_Se_from_h(h_f,vg_alpha,vg_m,vg_n); -// printf("Se_f = %8.6lf, Se_inverse = %8.6lf\n", Se_f, Se); -// } - -// // nint = number of "dh" intervals to integrate over using trapezoidal rule -// dh = (h_f-h_i)/(double)nint; -// Geff = 0.0; - -// // integrate k(h) dh from h_i to h_f, using trapezoidal rule, with subscript -// // 1 denoting the left-hand side of the trapezoid, and 2 denoting the right-hand side - -// Se1 = Se_i; // could just use Se_i in next statement. Done 4 completeness. -// K1 = calc_K_from_Se(Se1, Ksat, vg_m); -// h2 = h_i + dh; - -// for(int i=1; i<=nint; i++) { - -// Se2 = calc_Se_from_h(h2, vg_alpha, vg_m, vg_n); -// K2 = calc_K_from_Se(Se2, Ksat, vg_m); -// Geff += (K1+K2)*dh/2.0; // trapezoidal rule -// // reset for next time through loop -// K1 = K2; -// h2 += dh; -// } - -// Geff = fabs(Geff/Ksat); // by convention Geff is a positive quantity - -// if (verbosity.compare("high") == 0){ -// printf ("Capillary suction (G) = %8.6lf \n", Geff); -// } - -// return Geff; - -// } -// else { - -// double Se_f = calc_Se_from_theta(theta1,theta_e,theta_r); // the scaled moisture content of the wetting front -// double Se_i = calc_Se_from_theta(theta2,theta_e,theta_r); // the scaled moisture content below the wetting front -// double H_c = bc_psib_cm*((2+3*lambda) / (1+3*lambda)); // Green ampt capillary drive parameter, which can be used in the approximation of G with the Brooks-Corey model (See Ogden and Saghafian, 1997) - -// Geff = H_c*(pow(Se_i,(3+1/lambda))-pow(Se_f,(3+1/lambda)))/(1-pow(Se_f,(3+1/lambda))); - -// if (isinf(Geff)) { -// Geff = H_c; -// } -// if (isnan(Geff)) { -// Geff = H_c; -// } - -// return Geff; - -// } - -// } - - - extern double calc_Geff(bool use_closed_form_G, double theta1, double theta2, double theta_e, double theta_r, double vg_alpha, double vg_n, double vg_m, double h_min, double Ksat, int nint, double lambda, double bc_psib_cm)