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.

Sign up for GitHub

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

handle_invalid_soil_type #31

Merged
merged 8 commits into from
Sep 12, 2024
Merged

handle_invalid_soil_type #31

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
37bd2e3
modification to handle invalid soil_type such as waterbody, glacier, …
ajkhattak Sep 11, 2024
11e33f2
update tests/config files, and replace max_soil_type with max_valid_s…
ajkhattak Sep 11, 2024
03490b7
updated unittest config file.
ajkhattak Sep 11, 2024
1558b70
Update src/bmi_lgar.cxx
ajkhattak Sep 11, 2024
116affa
set max limit on the max_valid_soil_types (25 soils), removed +=0.
ajkhattak Sep 11, 2024
bdf904a
changed stod -> stoi for num_soil_types, removed 25 that was added as…
ajkhattak Sep 12, 2024
32690c3
unit correction (for invalid soiltype outputs).
ajkhattak Sep 12, 2024
470fda5
set MAX_NUM_SOIL_TYPE as both the upper limit the default value for n…
ajkhattak Sep 12, 2024
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
2 changes: 1 addition & 1 deletion configs/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ A detailed description of the parameters for model configuration (i.e., initiali
| forcing_resolution | double (scalar)| - | sec/min/hr | temporal resolution | - | timestep of the forcing data. Recommended value of 3600 seconds. |
| endtime | double (scalar)| >0 | sec, min, hr, d | simulation duration | - | time at which model simulation ends |
| layer_soil_type | int (1D array) | - | - | state variable | - | layer soil type (read from the database file soil_params_file) |
| max_soil_types | int | >1 | - | - | - | maximum number of soil types read from the file soil_params_file (default is set to 15) |
| max_valid_soil_types | int | >1 | - | - | - | maximum number of valid soil types read from the file soil_params_file (default is set to 12; model not valid for soil_type = waterbody, lava, glacier, etc.) |
ajkhattak marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It should be noted that if this value is more than the constant defined here

LGAR-C/include/all.hxx

Line 39 in 5120df4

#define MAX_NUM_SOIL_TYPES 16

then an error/exit condition exists. I don't think this the number impacts the implementation details of this PR, but it is a somewhat odd out-of-band validation that is relevant.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not seeing if MAX_NUM_SOIL_TYPES is used anywhere, I am going to push a change to put this as an upper limit.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I realized that it isn't used in any of the BMI code paths, for sure, but in the src/main.cxx there is this check:

LGAR-C/src/main.cxx

Line 258 in 5120df4

if(num_soil_types > (MAX_NUM_SOIL_TYPES -1))

| wilting_point_psi | double (scalar) | - | cm | state variable | - | wilting point (the amount of water not available for plants) used in computing AET. Suggested value is 15495.0 cm, corresponding to 15 atm. |
| field_capacity_psi | double (scalar) | - | cm | state variable | - | capillary head corresponding to volumetric water content at which gravity drainage becomes slower, used in computing AET. Suggested value is 340.9 cm for most soils, corresponding to 1/3 atm, and 103.3 cm for sands, corresponding to 1/10 atm. |
| use_closed_form_G | bool | true or false | - | - | - | determines whether the numeric integral or closed form for G is used; a value of true will use the closed form. This defaults to false. |
Expand Down
4 changes: 2 additions & 2 deletions configs/config_lasam_Bushland.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,8 @@ forcing_resolution=3600[sec]
ponded_depth_max=0[cm]
use_closed_form_G=false
layer_soil_type=16,17,18
max_soil_types=25
max_valid_soil_types=25
wilting_point_psi=15495.0[cm]
field_capacity_psi=340.9[cm]
giuh_ordinates=0.06,0.51,0.28,0.12,0.03
adaptive_timestep=true
adaptive_timestep=true
2 changes: 1 addition & 1 deletion configs/config_lasam_Phillipsburg.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ forcing_resolution=3600[sec]
ponded_depth_max=2[cm]
use_closed_form_G=false
layer_soil_type=13,14,15
max_soil_types=25
max_valid_soil_types=25
wilting_point_psi=15495.0[cm]
field_capacity_psi=340.9[cm]
giuh_ordinates=0.06,0.51,0.28,0.12,0.03
Expand Down
2 changes: 1 addition & 1 deletion configs/config_lasam_sft_ngen.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ forcing_resolution=3600[sec]
ponded_depth_max=2[cm]
use_closed_form_G=false
layer_soil_type=13,14,15
max_soil_types=25
max_valid_soil_types=25
wilting_point_psi=15495.0[cm]
field_capacity_psi=340.9[cm]
giuh_ordinates=0.06,0.51,0.28,0.12,0.03
Expand Down
3 changes: 2 additions & 1 deletion include/all.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -146,7 +146,8 @@ struct lgar_bmi_parameters
double *giuh_ordinates; // geomorphological instantaneous unit hydrograph
int num_giuh_ordinates; // number of giuh ordinates

int calib_params_flag = 0; // flag for calibratable parameters; if true, then calibratable params are updated otherwise not
int calib_params_flag = 0; // flag for calibratable parameters; if true, then calibratable params are updated otherwise not
bool is_invalid_soil_type = true; // checks if the provided soil type is valid for the model, if not then return Q_out = precip
};

// Define a data structure for local (timestep) and global mass balance parameters
Expand Down
36 changes: 34 additions & 2 deletions src/bmi_lgar.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,40 @@ Update()
std::cerr<<"|****************** LASAM BMI Update... ******************|\n";
std::cerr<<"---------------------------------------------------------\n";
}


double mm_to_cm = 0.1; // unit conversion

if (state->lgar_bmi_params.is_invalid_soil_type) {
// add to mass balance accumulated variables
state->lgar_mass_balance.volprecip_cm += state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm;
state->lgar_mass_balance.volin_cm += 0.0;
ajkhattak marked this conversation as resolved.
Show resolved Hide resolved
state->lgar_mass_balance.volon_cm = 0.0;
state->lgar_mass_balance.volend_cm = state->lgar_mass_balance.volstart_cm;;
ajkhattak marked this conversation as resolved.
Show resolved Hide resolved
state->lgar_mass_balance.volAET_cm += 0.0;
state->lgar_mass_balance.volrech_cm += 0.0;
state->lgar_mass_balance.volrunoff_cm += state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is volrunoff and volQ the same thing? If not, this looks like a mass balance logical error since the precipitation is being accounted for twice.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

currently they are the same, once we merge LGARTO, they will be different

state->lgar_mass_balance.volQ_cm += state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm;
state->lgar_mass_balance.volQ_gw_cm += 0.0;
state->lgar_mass_balance.volPET_cm += 0.0;
state->lgar_mass_balance.volrunoff_giuh_cm += 0.0;
state->lgar_mass_balance.volchange_calib_cm += 0.0;

// converted values, a struct local to the BMI and has bmi output variables
bmi_unit_conv.mass_balance_m = state->lgar_mass_balance.local_mass_balance * state->units.cm_to_m;
bmi_unit_conv.volprecip_timestep_m = state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm;
bmi_unit_conv.volin_timestep_m = 0.0;
bmi_unit_conv.volend_timestep_m = 0.0;
bmi_unit_conv.volAET_timestep_m = 0.0;
bmi_unit_conv.volrech_timestep_m = 0.0;
bmi_unit_conv.volrunoff_timestep_m = state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm;
bmi_unit_conv.volQ_timestep_m = state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm;
ajkhattak marked this conversation as resolved.
Show resolved Hide resolved
bmi_unit_conv.volQ_gw_timestep_m = 0.0;
bmi_unit_conv.volPET_timestep_m = 0.0;
bmi_unit_conv.volrunoff_giuh_timestep_m = 0.0;

return;
}

// if lasam is coupled to soil freeze-thaw, frozen fraction module is called
if (state->lgar_bmi_params.sft_coupled)
frozen_factor_hydraulic_conductivity(state->lgar_bmi_params);
Expand All @@ -95,7 +128,6 @@ Update()
state->lgar_bmi_params.calib_params_flag = false;
}

double mm_to_cm = 0.1; // unit conversion

// local variables for readibility
int subcycles;
Expand Down
24 changes: 16 additions & 8 deletions src/lgar.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -224,7 +224,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
bool is_wilting_point_psi_cm_set = false;
bool is_field_capacity_psi_cm_set = false;
bool is_soil_params_file_set = false;
bool is_max_soil_types_set = false;
bool is_max_valid_soil_types_set = false;
bool is_giuh_ordinates_set = false;
bool is_soil_z_set = false;
bool is_ponded_depth_max_cm_set = false;
Expand Down Expand Up @@ -347,9 +347,9 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)

continue;
}
else if (param_key == "max_soil_types") {
else if (param_key == "max_valid_soil_types") {
state->lgar_bmi_params.num_soil_types = stod(param_value);
ajkhattak marked this conversation as resolved.
Show resolved Hide resolved
is_max_soil_types_set = true;
is_max_valid_soil_types_set = true;
continue;
}
else if (param_key == "soil_params_file") {
Expand Down Expand Up @@ -533,8 +533,8 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
std::cerr<<" ***** \n";
}

if(!is_max_soil_types_set)
state->lgar_bmi_params.num_soil_types = 15; // maximum number of soil types defaults to 15
if(!is_max_valid_soil_types_set)
state->lgar_bmi_params.num_soil_types = 12; // maximum number of soil types defaults to 12

if (verbosity.compare("high") == 0) {
std::cerr<<"Maximum number of soil types: "<<state->lgar_bmi_params.num_soil_types<<"\n";
Expand All @@ -560,9 +560,17 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
wilting_point_psi_cm, state->soil_properties);

// check if soil layers provided are within the range
state->lgar_bmi_params.is_invalid_soil_type = false; // model not valid for soil types = waterbody, glacier, lava, etc.
for (int layer=1; layer <= state->lgar_bmi_params.num_layers; layer++) {
assert (state->lgar_bmi_params.layer_soil_type[layer] <= state->lgar_bmi_params.num_soil_types);
assert (state->lgar_bmi_params.layer_soil_type[layer] <= max_num_soil_in_file);
//assert (state->lgar_bmi_params.layer_soil_type[layer] <= state->lgar_bmi_params.num_soil_types);
//assert (state->lgar_bmi_params.layer_soil_type[layer] <= max_num_soil_in_file);
if (state->lgar_bmi_params.layer_soil_type[layer] > state->lgar_bmi_params.num_soil_types) {
ajkhattak marked this conversation as resolved.
Show resolved Hide resolved
state->lgar_bmi_params.is_invalid_soil_type = true;
std::cerr << "Invalid soil type: "
<< state->lgar_bmi_params.layer_soil_type[layer]
<<". Model returns input_precip = ouput_Qout. \n";
break;
}
}

if (verbosity.compare("high") == 0) {
Expand Down Expand Up @@ -2646,4 +2654,4 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm

}

#endif
#endif
4 changes: 2 additions & 2 deletions tests/configs/config_lasam_synth_0.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,8 @@ timestep=300[sec]
endtime=12[h]
forcing_resolution=3600[sec]
layer_soil_type=13,14,15
max_soil_types=15
max_valid_soil_types=15
wilting_point_psi=15495.0[cm]
field_capacity_psi=340.9[cm]
giuh_ordinates=0.06,0.51,0.28,0.12,0.03
ponded_depth_max=1.0[cm]
ponded_depth_max=1.0[cm]
2 changes: 1 addition & 1 deletion tests/configs/config_lasam_synth_1.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ timestep=300[sec]
endtime=12[hr]
forcing_resolution=300[sec]
layer_soil_type=13,14,15
max_soil_types=25
max_valid_soil_types=25
wilting_point_psi=15495.0[cm]
field_capacity_psi=340.9[cm]
giuh_ordinates=0.0,0.0,0.0,0.0,0.0
2 changes: 1 addition & 1 deletion tests/configs/config_lasam_synth_2.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ timestep=300[sec]
endtime=12[hr]
forcing_resolution=300[sec]
layer_soil_type=13,14,15
max_soil_types=25
max_valid_soil_types=25
wilting_point_psi=15495.0[cm]
field_capacity_psi=340.9[cm]
giuh_ordinates=0.0,0.0,0.0,0.0,0.0
2 changes: 1 addition & 1 deletion tests/configs/unittest.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ endtime=1.0[hr]
forcing_resolution=3600[sec]
ponded_depth_max=0[cm]
layer_soil_type=13,14,15
max_soil_types=15
max_valid_soil_types=15
wilting_point_psi=15495.0[cm]
field_capacity_psi=340.9[cm]
giuh_ordinates=0.06,0.51,0.28,0.12,0.03
Expand Down
Loading