From f754860f5af1d321a1e97efae76463ff249c9add Mon Sep 17 00:00:00 2001 From: Magnus Dahler Norling Date: Thu, 31 Oct 2024 11:57:06 +0100 Subject: [PATCH] Model reorganization --- dev_notes/todo.txt | 2 + ..._model_test.txt => magic_model_simple.txt} | 10 +-- .../modules/microplastic/mp_soil_column.txt | 61 ++++++++++++++++--- models/modules/nivafjord/basin.txt | 2 +- models/modules/nivafjord/point_sources.txt | 28 ++++++--- models/mp_soil_model.txt | 2 +- ...nivafjord_model.txt => nivafjord_base.txt} | 0 models/nivafjord_chem_base.txt | 4 +- models/nivafjord_isolated_basin_fpv_model.txt | 2 +- ...lake_model.txt => nivafjord_lake_base.txt} | 0 models/nivafjord_lake_simplytox_model.txt | 6 +- src/model_declaration.cpp | 7 +++ 12 files changed, 91 insertions(+), 33 deletions(-) rename models/{magic_model_test.txt => magic_model_simple.txt} (89%) rename models/{nivafjord_model.txt => nivafjord_base.txt} (100%) rename models/{nivafjord_lake_model.txt => nivafjord_lake_base.txt} (100%) diff --git a/dev_notes/todo.txt b/dev_notes/todo.txt index e0cf3d0c..99be6a98 100644 --- a/dev_notes/todo.txt +++ b/dev_notes/todo.txt @@ -56,6 +56,8 @@ NIVAFjord + Find a way to unify nivafjord_base and nivafjord_lake_base into one model declaration (using options) + Need a complete rethink of horizontal fluxes. They cause too much numerical instability when there is large in/outflow for small layers. There should really be an error if there is an opening to a neighboring basin on a depth below the max depth of that basin. diff --git a/models/magic_model_test.txt b/models/magic_model_simple.txt similarity index 89% rename from models/magic_model_test.txt rename to models/magic_model_simple.txt index 6354c219..adabf9db 100644 --- a/models/magic_model_test.txt +++ b/models/magic_model_simple.txt @@ -1,6 +1,6 @@ -model("MAGIC (test)") { +model("MAGIC (simple)") { ci : index_set("Compartment index") @@ -27,7 +27,6 @@ model("MAGIC (test)") { f : quantity("F") po4 : quantity("PO4") - ts : property("Time step size") concn : property("Ionic concentration") conc_all : property("Total concentration") @@ -46,7 +45,10 @@ model("MAGIC (test)") { load("modules/magic/magic_forest_cnp.txt", module("MAGIC-Forest CNP", comp, nh4, no3, po4, ts, ext_in, bal)) - # TODO: parametrize: - sol : solver("MAGIC solver", euler, [1/10, month], 0.01) + par_group("Solver resolution") { + step : par_real("Solver step size", [month], 0.1, 0.0001, 1) + } + + sol : solver("MAGIC solver", euler, step) solve(sol, comp.ca, comp.mg, comp.na, comp.k, comp.nh4, comp.so4, comp.cl, comp.no3, comp.f, comp.po4) } \ No newline at end of file diff --git a/models/modules/microplastic/mp_soil_column.txt b/models/modules/microplastic/mp_soil_column.txt index 8ed7b4bf..d79348a7 100644 --- a/models/modules/microplastic/mp_soil_column.txt +++ b/models/modules/microplastic/mp_soil_column.txt @@ -18,11 +18,20 @@ This is a module for transport of microplastic in a single *experimentally contr sph : par_real("Soil grain sphericity", [], 1, 0, 1) dp : par_real("Soil grain diameter", [m m], 0.5, 0.001, 10) - grad : par_real("Hydraulic gradient", [], 0.05, 0, 1, "When the box is full") + grad : par_real("Hydraulic gradient", [], 0.05, 0, 1, "When the box is full") #TODO: Do this differently + + roots : par_bool("Roots are present", false) #TODO: May also need the fraction of volume this takes up? } - par_group("General transport") { + par_group("General soil properties") { + fcs : par_real("Field capacity saturation", [], 0.2, 0.1, 0.4, "0.1=sand, 0.4=clay") + pwps : par_real("Permanent wilting point saturation", [], 0.09, 0.04, 0.22, "0.04=sand, 0.22=clay") + # Tabular values: + # https://extension.okstate.edu/fact-sheets/understanding-soil-water-content-and-thresholds-for-irrigation-management.html + } + + par_group("General transport") { transp_f_base : par_real("Base transport factor", [], 1e-5, 0, 1) bio_down : par_real("Bioturbation down rate", [day-1], 0.01, 0, 1) bio_up : par_real("Bioturbation up rate", [day-1], 0.001, 0, 1) @@ -40,11 +49,33 @@ This is a module for transport of microplastic in a single *experimentally contr load("stdlib/basic_math.txt", library("Response"), library("Basic")) precip : property("Precipitation") + pet : property("Potential evapotranspiration") theta : property("Porosity") perm : property("Soil permeability") cond : property("Soil hydraulic conductivity") + etp : property("Evapotranspiration") + etp_above : property("Evapotranspiration above") + + var(air.precip, [m m, day-1]) + + var(air.pet, [m m, day-1]) + + # TODO: May not be that realistic that etp works from top down since roots will always take some?? + var(layer.etp, [m m, day-1]) { + can_etp := is_at[vert.top] | roots, + { + maxcap := dz*theta, + pot := max(0, air.pet-etp_above), # Remaining unsatisfied pet + s_response(water, maxcap*pwps, maxcap*fcs, 0, pot) # Turn ET off as saturation goes below field capacity and approaches PWP. + } if can_etp, + 0 otherwise + } + + var(layer.etp_above, [m m, day-1]) { + etp[vert.above] + etp_above[vert.above] + } @no_store - var(air.precip, [m m, day-1], "Precipitation") + flux(layer.water, out, [m m, day-1], "Evapotranspiration flux") { etp } @no_store var(layer.theta, [], "Layer porosity") { theta0 } # Eventually allow change over time @@ -62,12 +93,13 @@ This is a module for transport of microplastic in a single *experimentally contr perc : property("Percolation") # Hmm this is done a bit awkward.. - water_tot : property("Water above") - var(layer.water_tot, [m m]) { water + water_tot[vert.above] } @no_store + water_tot : property("Water above") #Not actually above, but above including this. + var(layer.water_tot, [m m]) { max(0, water - dz*theta*fcs) + water_tot[vert.above] } @no_store z : property("Effective depth") var(layer.z, [m m]) { dz*theta + z[vert.above] } @no_store var(layer.water.side, [m m, day-1]) { + # TODO: This is wrong... Should compute gradient based on z and, z of outside and distance. agrad := grad*(water_tot/z), agrad*cond->> } @@ -77,13 +109,19 @@ This is a module for transport of microplastic in a single *experimentally contr maxcap_b := dz[vert.below] * theta[vert.below], maxrate := cond->>, + mincap := maxcap*fcs, + + maxperc := s_response(water, mincap, min(1.1*mincap, maxcap), 0, maxrate), + #s_response(water[vert.below], maxcap_b, 0.9*maxcap_b, 0, maxperc) + #1.1*mincap, mincap, 0, maxrate), # Cut flow off when water saturation reaches field capacity + s_response(water[vert.below], 0.9*maxcap_b, maxcap_b, maxperc, 0) # Cut flow off when water saturation below fills up - # TODO: This is just a test. - maxperc := s_response(water, 0, 0.1*maxcap, 0, maxrate), - s_response(water[vert.below], 0.8*maxcap_b, maxcap_b, maxperc, 0) + # TODO: Capillary action when water below is above fc and this layer is below fc? + # We could also compute forces on the water, but probably not necessary + # https://onlinelibrary.wiley.com/doi/epdf/10.1155/2020/6671479 } - var(layer.water, [m m], "Soil layer pore water") @initial { 0[m m] } # Assume it is dry. + var(layer.water, [m m], "Soil layer pore water") @initial { dz*theta*fcs } flux(out, layer.water[vert.top], [m m, day-1], "Precipitation to soil") { air.precip } @@ -91,7 +129,10 @@ This is a module for transport of microplastic in a single *experimentally contr flux(layer.water, vert, [m m, day-1], "Layer water percolation") { perc } - # TODO: Capillarity, retention of water in a layer, evapotranspiration. + + + # Need, min saturation (=field capacity), max saturation (=1?) + # Also permanent wilting point (etp cutoff). diff --git a/models/modules/nivafjord/basin.txt b/models/modules/nivafjord/basin.txt index 77df6751..338a7728 100644 --- a/models/modules/nivafjord/basin.txt +++ b/models/modules/nivafjord/basin.txt @@ -6,7 +6,7 @@ preamble("NIVAFjord dimensions", version(0, 0, 0), par_group("Layer thickness", all_layer) { dz0 : par_real("Stable layer thickness", [m], 1, 0.1, 5) - } @fully_distribute + } @fully_distribute # This is because an external_computation uses this parameter, and it can't change its indexing scheme. par_group("Layer area", layer) { A : par_real("Layer area", [m 2], 10000, 0, 1e6) diff --git a/models/modules/nivafjord/point_sources.txt b/models/modules/nivafjord/point_sources.txt index 22ec499c..555e38be 100644 --- a/models/modules/nivafjord/point_sources.txt +++ b/models/modules/nivafjord/point_sources.txt @@ -3,8 +3,10 @@ module("NIVAFjord point sources", version(0, 0, 0), layer : compartment, water : quantity, din : quantity, + on : quantity, phos : quantity, sed : quantity, + oc : quantity, o2 : quantity, heat : quantity, temp : property @@ -14,36 +16,42 @@ This is a module for point source inputs to NIVAFjord. """ par_group("Point sources", layer) { - eff_q : par_real("Water point source", [m 3, day-1], 0) - eff_n : par_real("DIN point source", [k g, day-1], 0) - eff_p : par_real("DIP point source", [k g, day-1], 0) - eff_ss : par_real("SS point source", [k g, day-1], 0) + eff_q : par_real("Water point source", [m 3, year-1], 0) + eff_in : par_real("DIN point source", [M g, year-1], 0) #Maybe we should introduce the shorthand ton for M g. + eff_on : par_real("DON point source", [M g, year-1], 0) + eff_p : par_real("DIP point source", [M g, year-1], 0) + eff_ss : par_real("SS point source", [M g, year-1], 0) + ss_c : par_real("SS point source C fraction", [g, k g-1], 0) o2sat : par_real("Point source water O2 saturation", [], 1) } load("stdlib/seawater.txt", library("Sea oxygen")) load("stdlib/physiochemistry.txt", library("Chemistry")) - flux(out, layer.water, [m 3, day-1], "Point source water to layer") { eff_q } + flux(out, layer.water, [m 3, day-1], "Point source water to layer") { eff_q/time.days_this_year } flux(out, layer.water.heat, [J, day-1], "Point source heat to layer") { - conc(layer.water.heat)*eff_q # Assume same temperature in effluent and bay. TODO: This is not always the case. + conc(layer.water.heat)*eff_q/time.days_this_year # Assume same temperature in effluent and bay. TODO: This is not always the case. } flux(out, layer.water.o2, [k g, day-1], "Point source O2 to layer") { eff_temp := layer.water.temp, eff_salin := 0, o2c := o2sat*o2_mol_mass*o2_saturation(eff_temp, eff_salin), - o2c*eff_q->> + o2c*eff_q/time.days_this_year->> } # TODO: CO2, CH4 content of effluent water? - flux(out, layer.water.din, [k g, day-1], "Point source DIN to layer") { eff_n } + flux(out, layer.water.din, [k g, day-1], "Point source DIN to layer") { eff_in/time.days_this_year->> } - flux(out, layer.water.phos, [k g, day-1], "Point source DIP to layer") { eff_p } + flux(out, layer.water.on, [k g, day-1], "Point source DON to layer") { eff_on/time.days_this_year->> } - flux(out, layer.water.sed, [k g, day-1], "Point source suspended solids to layer") { eff_ss } + flux(out, layer.water.phos, [k g, day-1], "Point source DIP to layer") { eff_p/time.days_this_year->> } + + flux(out, layer.water.sed, [k g, day-1], "Point source suspended solids to layer") { eff_ss/time.days_this_year->> } + + flux(out, layer.water.sed.oc, [k g, day-1], "Point source POC to layer") { eff_ss*ss_c/time.days_this_year->> } # TODO: What about carbon content of the effluent suspended solids } \ No newline at end of file diff --git a/models/mp_soil_model.txt b/models/mp_soil_model.txt index 24ae21d4..b91df99a 100644 --- a/models/mp_soil_model.txt +++ b/models/mp_soil_model.txt @@ -21,7 +21,7 @@ This is a model for simulating microplastic movement in a single experimentally module("Microplastic soil column", air, layer, water, mp, vert) ) - sol : solver("Soil solver", inca_dascru, [10, min], 1e-2) + sol : solver("Soil solver", inca_dascru, [10, min], 1e-5) solve(sol, layer.water, layer.mp) } \ No newline at end of file diff --git a/models/nivafjord_model.txt b/models/nivafjord_base.txt similarity index 100% rename from models/nivafjord_model.txt rename to models/nivafjord_base.txt diff --git a/models/nivafjord_chem_base.txt b/models/nivafjord_chem_base.txt index d936ca03..38b22a15 100644 --- a/models/nivafjord_chem_base.txt +++ b/models/nivafjord_chem_base.txt @@ -6,7 +6,7 @@ model("NIVAFjord-Chem-Base") { This won't run as its own model, it is just a base for model extension. """ - extend("nivafjord_model.txt") + extend("nivafjord_base.txt") phyt_type : index_set("Phytoplankton type") @@ -42,7 +42,7 @@ This won't run as its own model, it is just a base for model extension. module("NIVAFjord boundary chemistry", bnd_layer, water, o2, oc, din, on, phos, op, fjord_phyt, chl_a, sed, chem_par)) load("modules/nivafjord/point_sources.txt", - module("NIVAFjord point sources", layer, water, din, phos, sed, o2, heat, temp)) + module("NIVAFjord point sources", layer, water, din, on, phos, sed, oc, o2, heat, temp)) load("modules/easychem.txt", module("Simple River O₂", river, water, o2, temp)) diff --git a/models/nivafjord_isolated_basin_fpv_model.txt b/models/nivafjord_isolated_basin_fpv_model.txt index 4de80566..963d97ee 100644 --- a/models/nivafjord_isolated_basin_fpv_model.txt +++ b/models/nivafjord_isolated_basin_fpv_model.txt @@ -9,7 +9,7 @@ It is also currently without horizontal exchange between basins, so it only work # TODO: FPV cover does not dynamically affect wind mixing (not sure to how large a degree that would happen). - extend("nivafjord_lake_model.txt") + extend("nivafjord_lake_base.txt") tox_index : index_set("Contaminant type") diff --git a/models/nivafjord_lake_model.txt b/models/nivafjord_lake_base.txt similarity index 100% rename from models/nivafjord_lake_model.txt rename to models/nivafjord_lake_base.txt diff --git a/models/nivafjord_lake_simplytox_model.txt b/models/nivafjord_lake_simplytox_model.txt index f398605e..a5410cd6 100644 --- a/models/nivafjord_lake_simplytox_model.txt +++ b/models/nivafjord_lake_simplytox_model.txt @@ -2,9 +2,7 @@ model("NIVAFjord lake SimplyTox") { """ -TODO: This one is incomplete. - -Maaybe we want a simpler biochemistry model for this one, not the entire NIVAFjord biochemistry... +TODO: Maybe we want a simpler biochemistry model for this one, not the entire NIVAFjord biochemistry... """ extend("simplytox_model.txt") @exclude( @@ -20,7 +18,7 @@ Maaybe we want a simpler biochemistry model for this one, not the entire NIVAFjo sol : solver ) - extend("nivafjord_lake_model.txt") + extend("nivafjord_lake_base.txt") downstream : connection("Downstream") @directed_graph { river+ } @no_cycles diff --git a/src/model_declaration.cpp b/src/model_declaration.cpp index 50b69cfc..b75afc7e 100644 --- a/src/model_declaration.cpp +++ b/src/model_declaration.cpp @@ -1304,6 +1304,13 @@ Solver_Registration::process_declaration(Catalog *catalog) { solver_fun = scope->resolve_argument(Reg_Type::solver_function, decl->args[1]); hmin = 0.01; + // TODO: There is a weakness with the matcher that it can't know the type of an identifier, so an identifier matches any Decl_Type. + // Hence this hack.. + auto h_tmp = scope->resolve_argument(Reg_Type::unrecognized, decl->args[2]); + if(h_tmp.reg_type == Reg_Type::parameter) + which += 2; + + if(which == 0 || which == 1) { h_unit = scope->resolve_argument(Reg_Type::unit, decl->args[2]); if(which == 1)