Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

updating LASAM to take field capacity as a parameter in the config file #22

Merged
merged 15 commits into from
Apr 16, 2024
Merged
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
updated unit test to test calibratable parameters for field capacity …
…and maximum ponded head
Peter La Follette authored and Peter La Follette committed Apr 12, 2024
commit a846fca7ed4d6be109c8510360dd3e1adeb10ce0
4 changes: 2 additions & 2 deletions include/all.hxx
Original file line number Diff line number Diff line change
@@ -192,8 +192,8 @@ struct lgar_calib_parameters
double *vg_n; // Van Genuchten n [-]
double *vg_alpha; // Van Genuchten alpha [1/cm]
double *Ksat; // Hydraulic conductivity [cm/hr]
double *field_capacity_psi; // field capacity in capillary head [cm]
double *ponded_depth_max; // maximum ponded depth of surface water [cm]
double field_capacity_psi; // field capacity in capillary head [cm]
double ponded_depth_max; // maximum ponded depth of surface water [cm]

};

12 changes: 11 additions & 1 deletion src/bmi_lgar.cxx
Original file line number Diff line number Diff line change
@@ -524,6 +524,8 @@ update_calibratable_parameters()
<<", vg_n = " << state->soil_properties[soil].vg_n
<<", vg_alpha = " << state->soil_properties[soil].vg_alpha_per_cm
<<", Ksat = " << state->soil_properties[soil].Ksat_cm_per_h
<<", field_capacity_psi = " << state->lgar_bmi_params.field_capacity_psi_cm
<<", ponded_depth_max = " << state->lgar_bmi_params.ponded_depth_max_cm
<<", theta = " << current->theta <<"\n";
}

@@ -533,6 +535,8 @@ update_calibratable_parameters()
state->soil_properties[soil].vg_m = 1.0 - 1.0/state->soil_properties[soil].vg_n;
state->soil_properties[soil].vg_alpha_per_cm = state->lgar_calib_params.vg_alpha[layer_num-1];
state->soil_properties[soil].Ksat_cm_per_h = state->lgar_calib_params.Ksat[layer_num-1];
state->lgar_bmi_params.field_capacity_psi_cm = state->lgar_calib_params.field_capacity_psi;
state->lgar_bmi_params.ponded_depth_max_cm = state->lgar_calib_params.ponded_depth_max;

current->theta = calc_theta_from_h(current->psi_cm, state->soil_properties[soil].vg_alpha_per_cm,
state->soil_properties[soil].vg_m, state->soil_properties[soil].vg_n,
@@ -546,6 +550,8 @@ update_calibratable_parameters()
<<", vg_n = " << state->soil_properties[soil].vg_n
<<", vg_alpha = " << state->soil_properties[soil].vg_alpha_per_cm
<<", Ksat = " << state->soil_properties[soil].Ksat_cm_per_h
<<", field_capacity_psi = " << state->lgar_bmi_params.field_capacity_psi_cm
<<", ponded_depth_max = " << state->lgar_bmi_params.ponded_depth_max_cm
<<", theta = " << current->theta <<"\n";
}

@@ -582,7 +588,7 @@ GetVarGrid(std::string name)
|| name.compare("actual_evapotranspiration") == 0) // double
return 1;
else if (name.compare("surface_runoff") == 0 || name.compare("giuh_runoff") == 0
|| name.compare("soil_storage") == 0) // double
|| name.compare("soil_storage") == 0 || name.compare("field_capacity") == 0 || name.compare("ponded_depth_max") == 0)// double
return 1;
else if (name.compare("total_discharge") == 0 || name.compare("infiltration") == 0
|| name.compare("percolation") == 0 || name.compare("groundwater_to_stream_recharge") == 0) // double
@@ -813,6 +819,10 @@ GetValuePtr (std::string name)
return (void*)this->state->lgar_calib_params.vg_alpha;
else if (name.compare("hydraulic_conductivity") == 0)
return (void*)this->state->lgar_calib_params.Ksat;
else if (name.compare("ponded_depth_max") == 0)
return (void*)&this->state->lgar_calib_params.field_capacity_psi;
else if (name.compare("field_capacity") == 0)
return (void*)&this->state->lgar_calib_params.ponded_depth_max;
else {
std::stringstream errMsg;
errMsg << "variable "<< name << " does not exist";
1 change: 1 addition & 0 deletions tests/configs/unittest.txt
Original file line number Diff line number Diff line change
@@ -6,6 +6,7 @@ initial_psi=2000.0[cm]
timestep=300[sec]
endtime=1.0[hr]
forcing_resolution=3600[sec]
ponded_depth_max=0[cm]
layer_soil_type=13,14,15
max_soil_types=15
wilting_point_psi=15495.0[cm]
39 changes: 33 additions & 6 deletions tests/main_unit_test_bmi.cxx
Original file line number Diff line number Diff line change
@@ -453,8 +453,8 @@ int main(int argc, char *argv[])

// Benchmark values of wetting fronts depth and moisture (b is for benchmark)
//std::vector<double> depth_wf_b = {1.873813, 44.00,175.0, 200.0}; // in cm
std::vector<double> depth_wf_b = {4.55355, 44.00,175.0, 200.0}; // in cm
std::vector<double> theta_wf_b = {0.213716, 0.172703948143618, 0.252113867764474, 0.179593529195751};
std::vector<double> depth_wf_b = {4.55332255489760840, 44.00,175.0, 200.0}; // in cm
std::vector<double> theta_wf_b = {0.21414942004659104, 0.172703948143618, 0.252113867764474, 0.179593529195751};

int m_to_cm = 100;
int m_to_mm = 1000;
@@ -518,7 +518,7 @@ int main(int argc, char *argv[])

// check total infiltration, AET, and PET.
double infiltration_check_mm = 1.896; // in mm
double AET_check_mm = 0.029801; // in mm
double AET_check_mm = 0.00922002885520525; // in mm
double PET_check_mm = 0.104; // in mm
double infiltration_computed = 0.0;
double PET_computed = 0.0;
@@ -562,8 +562,8 @@ int main(int argc, char *argv[])
std::stringstream errMsg;
errMsg << "Error between benchmark and simulated AET is "<< fabs(AET_check_mm - AET_computed * m_to_mm)
<< " which is unexpected. \n";
printf("benchmark AET: %lf \n", AET_check_mm);
printf("computed AET: %lf \n", AET_computed*m_to_mm);
printf("benchmark AET: %.17lf \n", AET_check_mm);
printf("computed AET: %.17lf \n", AET_computed*m_to_mm);

throw std::runtime_error(errMsg.str());
}
@@ -582,36 +582,61 @@ int main(int argc, char *argv[])
double *vg_n = new double[num_layers];
double *vg_alpha = new double[num_layers];
double *Ksat = new double[num_layers];
double field_capacity;
double ponded_depth_max;

double smcmax_set[] = {0.3513, 0.3773, 0.3617};
double vg_n_set[] = {1.44260592334, 1.14965918354, 1.39051695249};
double vg_alpha_set[] = {0.0021297, 0.0073272, 0.0027454};
double Ksat_set[] = {0.446, 0.0743, 0.415};
double field_capacity_set = 340.9;
double ponded_depth_max_set = 0.0;

// Get the initial values set through the config file
model_calib.GetValue("smcmax", &smcmax[0]);
model_calib.GetValue("van_genuchten_n", &vg_n[0]);
model_calib.GetValue("van_genuchten_alpha", &vg_alpha[0]);
model_calib.GetValue("hydraulic_conductivity", &Ksat[0]);
model_calib.GetValue("field_capacity", &field_capacity);
model_calib.GetValue("ponded_depth_max", &ponded_depth_max);

for (int i=0; i < num_layers; i++)
std::cout<<"| Initial values: layer = "<< i+1 <<", smcmax = "<< smcmax[i]
<<", vg_n = "<< vg_n[i] <<", vg_alpha = " << vg_alpha[i]
<<", Ksat = "<< Ksat[i] <<"\n";
printf("field_capacity: %lf \n", field_capacity);
printf("ponded_depth_max: %lf \n", ponded_depth_max);

// set the new values
model_calib.SetValue("smcmax", &smcmax_set[0]);
model_calib.SetValue("van_genuchten_n", &vg_n_set[0]);
model_calib.SetValue("van_genuchten_alpha", &vg_alpha_set[0]);
model_calib.SetValue("hydraulic_conductivity", &Ksat_set[0]);
model_calib.SetValue("field_capacity", &field_capacity_set);
model_calib.SetValue("ponded_depth_max", &ponded_depth_max_set);

// get the new/updated values
model_calib.GetValue("smcmax", &smcmax[0]);
model_calib.GetValue("van_genuchten_n", &vg_n[0]);
model_calib.GetValue("van_genuchten_alpha", &vg_alpha[0]);
model_calib.GetValue("hydraulic_conductivity", &Ksat[0]);
model_calib.GetValue("field_capacity", &field_capacity);
model_calib.GetValue("ponded_depth_max", &ponded_depth_max);


if (fabs(ponded_depth_max - ponded_depth_max_set) > 1.E-5) {
std::stringstream errMsg;
errMsg << "Mismatch between ponded_depth_max calibrated values set and get "<< ponded_depth_max_set<<" "<< ponded_depth_max
<< " which is unexpected. \n";
throw std::runtime_error(errMsg.str());
}

if (fabs(field_capacity - field_capacity_set) > 1.E-5) {
std::stringstream errMsg;
errMsg << "Mismatch between field_capacity calibrated values set and get "<< field_capacity_set<<" "<< field_capacity
<< " which is unexpected. \n";
throw std::runtime_error(errMsg.str());
}

for (int i=0; i < num_layers; i++) {

if (fabs(smcmax[i] - smcmax_set[i]) > 1.E-5) {
@@ -649,6 +674,8 @@ int main(int argc, char *argv[])
std::cout<<"| Calib. values: layer = "<< i+1 <<", smcmax = "<< smcmax[i]
<<", vg_n = "<< vg_n[i] <<", vg_alpha = " << vg_alpha[i]
<<", Ksat = "<< Ksat[i] <<"\n";
printf("field_capacity = %lf \n", field_capacity);
printf("ponded_depth_max = %lf \n", ponded_depth_max);
std::cout<<"| *************************************** \n";
std::cout<<"| LASAM Calibration test passed? YES \n";
std::cout<<RESET<<"\n";