diff --git a/.gitignore b/.gitignore index ae2f9d700..dc3262012 100644 --- a/.gitignore +++ b/.gitignore @@ -64,3 +64,4 @@ tutorial/jupyter/.ipynb_checkpoints/ # images *.png .vscode/settings.json + diff --git a/experimental/analytical_model_perirhizal/test_perirhizal.py b/experimental/analytical_model_perirhizal/test_perirhizal.py new file mode 100644 index 000000000..38ed2d65b --- /dev/null +++ b/experimental/analytical_model_perirhizal/test_perirhizal.py @@ -0,0 +1,870 @@ + +# this is a test file by Erik Kopp to test the alternative implementation of the perirhizal resistances (water) and solute flow + +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt +import time + +import sys; sys.path.append("../.."); sys.path.append("../../src/") +import plantbox as pb +#from plantbox import Perirhizal +from plantbox.functional.Perirhizal import PerirhizalPython +import plantbox.functional.van_genuchten as vg +from plantbox.visualisation import figure_style +#import Perirhizal + + +# dumux rosi imports +#sys.path.append("../../dumux-rosi/build-cmake/cpp/python_binding/") +#sys.path.append("../../../dumuxtest/dumux/dumux-rosi/build-cmake/cpp/python_binding/") +from rosi.richards_flat import RichardsFlatWrapper as RichardsWrapper # Python part, macroscopic soil model +#from richards import RichardsWrapper # Python part +from rosi.richards_no_mpi import RichardsNoMPIWrapper # Python part of cylindrcial +from rosi.rosi_richardsnc_cyl import RichardsNCCylFoam # C++ part (Dumux binding), macroscopic soil model +#from rosi_richardsnc_cyl import RichardsNCCylFoam as RichardsNC_cyl # C++ part (Dumux binding) +#from rosi_richards22c import RichardsNCSPILU as RichardsNCSP #test + + + +#numerics +import math +from scipy.optimize import fsolve, root_scalar +from numpy import linalg as LA + +# run the dumux implementation of root water and nitrate uptake, later compare it to the analytical approximation + +n_tests = 1 #try everything here for this many random parameter sets +do_computation = False #should the computation be run or take the data from a saved file +showuniform = True #should the uniform implementation be shown +showTiina = True #should Tiinas implementation be shown + +# general parameters +max_time = 10 # d +n_times = 1000+1 # number of time slots, -1 for intervals +r_prhiz = 0.6 # perirhizal radius[cm], computed for a RLD above 1cm/cm3 +r_root = 0.02 # root radius [cm] +NC = 10 # number of spatial discretisations +length = 1 #default length of the segment, will not change the outcpme as all variables are assumed constant in this direction [cm] + +#three scenarios will be computed: one without inflow, one with an advective flow (Dirichlet BC), another with a Dirichlet BC +n_scenarios = 3 + +#initial conditions +initial_waterpotential = -100 +initial_soluteconcentration = 1.6e-6#mol/cm3, 103mg/L of NO3 (one of the Tereno measurements in 2015, TODO: look for another source) leads to slightly above 1.6 + +#space for the outputs +# number of tests, number of scenarios, number of timesteps, (rootuptake, inflow, discretisation) +#water outputs +#simulations in dumux +watercontent_dumux = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #watercontent in cm3/cm3, (radial) uptake and inflow in cm2/d +waterpotential_dumux = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #waterpotential in cm, (radial) uptake and inflow in cm2/d +#analytical approximations +watercontent_sr = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #watercontent in cm3/cm3, (radial) uptake and inflow in cm2/d +waterpotential_sr = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #waterpotential in cm, (radial) uptake and inflow in cm2/d + +#solute outputs +#simulations in dumux +solutes_dumux = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) +#analytical approximations +#steady state approximations +solutes_ss = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) +#simple steady rate approximations +solutes_sr_simp = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) +#steady rate approximations +solutes_sr = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) +#general steady rate with far field approximation, not used here as there is no scenario for it +#solutes_ff = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) +#dirichlet outer BC +solutes_d = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) +#assume a uniform solute concentration +solutes_u = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) +#Tiina Roose +solutes_TR = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) + + +#discretisation +lb = 0.5 +points = np.logspace(np.log(r_root) / np.log(lb), np.log(r_prhiz) / np.log(lb), NC+1, base = lb) +CC = np.array([(points[i] + points[i+1])/2 for i in range(NC)]) +volumes = np.array([(points[i+1]**2 - points[i]**2)*3.14 for i in range(NC)]) + +soilVG = [0.078, 0.43, 0.036, 1.56, 24.96] # hydrus loam soil +sp = vg.Parameters(soilVG) + +def run_perirhizal_test(max_time, n_times, r_prhiz, r_root, NC, points, CC, volumes, length, n_scenarios, initial_waterpotential, initial_soluteconcentration): + + #space for the outputs + watercontent_dumux = np.zeros((n_scenarios, n_times, NC+2)) #watercontent in cm3/cm3, (radial) uptake and inflow in cm/d + waterpotential_dumux = np.zeros((n_scenarios, n_times, NC+2)) #waterpotential in cm, (radial) uptake and inflow in cm/d + watercontent_sr = np.zeros((n_scenarios, n_times, NC+2)) #watercontent in cm3/cm3, (radial) uptake and inflow in cm/d + waterpotential_sr = np.zeros((n_scenarios, n_times, NC+2)) #waterpotential in cm, (radial) uptake and inflow in cm/d + solutes_dumux = np.zeros((n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, (radial) uptake and inflow in mol/(cm2d) + solutes_ss = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) + solutes_sr_simp = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) + solutes_sr = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) + solutes_ff = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) + solutes_d = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) + solutes_u = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) + solutes_TR = np.zeros((n_tests, n_scenarios, n_times, NC+2)) #solute concentration in mol/cm3, uptake and inflow in mol/(cm2d) + + simtimes = np.linspace(0,max_time,n_times)#[1:] + dt = max_time / n_times + rho = r_prhiz / r_root + + #soil parameters + soilVG = [0.078, 0.43, 0.036, 1.56, 24.96] # hydrus loam soil + + # root conductivity and solute uptake parameters, constant throughout the entire simulation time + molarMassWater = 18 #g/mol + molarMassSolute = 62 #g/mol, NO3 + root_conductivity = 1e-4 #1/d + inner_kr = root_conductivity * r_root * 2 * 3.14 + waterdemand = -0.05 #cm/d + radial_waterdemand = 2*3.14*r_root * waterdemand #cm2/d + Vmax = 4.0e-11 * (2*3.14*r_root) * (24*3600) #mol/(cm2 s) * cm * cm * (s/d) -> mol / (cm d) + Vmax_per_area = Vmax / (2*3.14*r_root) #mol / (cm d) /cm2 = mol/(cm2d) + Km = 1.5e-7 #mol/cm3 + + #diffusion coefficient of nitrate + Ds = 1.902e-5 * 24 * 3600 #cm2/s -> cm2/d + + + # load the perirhizal model + peri = PerirhizalPython() + sp = vg.Parameters(soilVG) + peri.set_soil(sp) + #no lookup tables are used here as there arent many simulations + + + + + # initialise the dumux models for the scenarios + s_nf = RichardsWrapper(RichardsNCCylFoam()) #no flux outer BC + s_af = RichardsWrapper(RichardsNCCylFoam()) #dirichlet BC for water, advective flow for solutes + s_d = RichardsWrapper(RichardsNCCylFoam()) #Dirichlet outer BC + + for s in [s_nf, s_af, s_d]: + s.initialize() + s.createGrid1d(points, length = length/100) # [m] -> [cm] + s.setVGParameters([soilVG]) + s.setHomogeneousIC(initial_waterpotential) # cm pressure head + + s.setBotBC("constantFluxCyl",waterdemand) # "noFlux")# Flux in cm/d + if s == s_nf: + s.setTopBC("constantFluxCyl",0.0) # [cm/day] "noFlux")#default, will be changed for one scenario + s.setParameter("Soil.BC.Top.SType", "3") # constantFluxCyl=3 (SType = Solute Type) + s.setParameter("Soil.BC.Top.CValue", "0.0") + else: + if s == s_af: + s.setParameter("Soil.BC.Top.Type", "11") #dummy Dirichlet BC so dumux does not get mad + s.setParameter("Soil.BC.Top.Value", str(initial_waterpotential)) + s.setParameter("Soil.BC.Top.SType", "10") #advective flow + s.setParameter("Soil.BC.Top.CValue", str(0)) + else: + s.setParameter("Soil.BC.Top.Type", "11") #dummy Dirichlet BC so dumux does not get mad + s.setParameter("Soil.BC.Top.Value", str(initial_waterpotential)) + s.setParameter("Soil.BC.Top.SType", "11") #dummy Dirichlet BC so dumux does not get mad + s.setParameter("Soil.BC.Top.CValue", str(initial_soluteconcentration*molarMassSolute)) + s.setParameter("Soil.BC.Bot.SType", "8") # michaelisMenten=8 (SType = Solute Type) + s.setParameter("Soil.BC.Bot.CValue", "0.0") #should not matter + s.setParameter("RootSystem.Uptake.Vmax", str(Vmax_per_area*molarMassSolute)) #mol/(cm2d) -> g/(cm2 d) + #s.setParameter("RootSystem.Uptake.Vmax", str(Vmax*molarMassSolute)) #mol/d -> g/d #TODO: Vmax or Vmax per area? + s.setParameter("RootSystem.Uptake.Km", str(Km*molarMassSolute)) # mol/cm3 -> g/cm3 + s.setParameter("Soil.IC.C", str(initial_soluteconcentration*molarMassSolute)) # g / cm3 + s.setParameter("Component.MolarMass", str(molarMassWater/1000)) #g/mol -> kg/mol water + s.setParameter("1.Component.MolarMass", str(molarMassSolute/1000)) #g/mol -> kg/mol nitrate + s.setParameter("1.Component.LiquidDiffusionCoefficient", str(Ds / 1.e4 / (24*3600))) #cm^2/s -> m^2/s + s.initializeProblem(maxDt = 0.01) + s.ddt = 1.e-4 # days + + cellVolumes = s_d.getCellSurfacesCyl() * length # cm3 + + #initial solute concentration for the analytical approximations + #simplified steady rate no flux (1d lookup table) + solutes_sr_simp = np.zeros((n_scenarios, n_times, NC+2)) + solutes_sr_simp[0,0,1]= initial_soluteconcentration + solutes_sr_simp[1,0,1]= initial_soluteconcentration + solutes_sr_simp[2,0,1]= initial_soluteconcentration + + #steady state solute flow + solutes_ss = np.zeros((n_scenarios, n_times, NC+2)) + solutes_ss[1,0,1]= initial_soluteconcentration + solutes_ss[2,0,1]= initial_soluteconcentration + + #steady rate no flux outer BC solute flow in the general water flow + solutes_sr = np.zeros((n_scenarios, n_times, NC+2)) + solutes_sr[0,0,1]= initial_soluteconcentration + solutes_sr[1,0,1]= initial_soluteconcentration + solutes_sr[2,0,1]= initial_soluteconcentration + + #Dirichlet BC + solutes_d = np.zeros((n_scenarios, n_times, NC+2)) + solutes_d[2,0,1]= initial_soluteconcentration + + #uniform concentration (= no analytical approximation) + solutes_u = np.zeros((n_scenarios, n_times, NC+2)) + solutes_u[0,0,1]= initial_soluteconcentration + solutes_u[1,0,1]= initial_soluteconcentration + solutes_u[2,0,1]= initial_soluteconcentration + + #Tiina Roose approximation + solutes_TR = np.zeros((n_scenarios, n_times, NC+2)) + + #variables for the no flux case + mean_watercontent_nf = 0.1 #cm3/cm3 + mean_waterpotential_nf = -100 #cm + mean_solutecontent_nf = 0 #mol/cm3 + Phi_soil_nf = 1.0 #mfp of the mean soil + Phi_root_nf = 1.0 #mfp next to the root + #mfp parameters for the perirhizal model + Phi_outer_nf = 1.0 #mfp at the outer perirhizal radius + rootuptake_w_nf = 1.0 #water uptake of the root, cm / d + inflow_w_nf = 0.0 #water uptake from outside the perirhizal zone, cm/d + rootuptake_s_nf = 0.0 #solute uptake of the root, mol / cm2d + inflow_s_nf = 0.0 #solute uptake from outside the perirhizal zone, mol / cm2d + + #variables for the general case + mean_watercontent_af = 0.1 #cm3/cm3 + mean_waterpotential_af = -100 #cm + mean_solutecontent_af = 0 #mol/cm3 + Phi_soil_af = 1.0 #mfp of the mean soil + Phi_root_af = 1.0 #mfp next to the root + #mfp parameters for the perirhizal model + Phi_outer_af = 1.0 #mfp at the outer perirhizal radius + rootuptake_w_af = 1.0 #water uptake of the root, cm / d + inflow_w_af = 0.0 #water uptake from outside the perirhizal zone, cm/d + rootuptake_s_af = 0.0 #solute uptake of the root, mol / cm2d + inflow_s_af = 0.0 #solute uptake from outside the perirhizal zone, mol / cm2d + + #variables for the Dirichlet case + mean_watercontent_d = 0.1 #cm3/cm3 + mean_waterpotential_d = -100 #cm + mean_solutecontent_d = 0 #mol/cm3 + Phi_soil_d = 1.0 #mfp of the mean soil + Phi_root_d = 1.0 #mfp next to the root + #mfp parameters for the perirhizal model + Phi_outer_d = 1.0 #mfp at the outer perirhizal radius + rootuptake_w_d = 1.0 #water uptake of the root, cm / d + inflow_w_d = 0.0 #water uptake from outside the perirhizal zone, cm/d + rootuptake_s_d = 0.0 #solute uptake of the root, mol / cm2d + inflow_s_d = 0.0 #solute uptake from outside the perirhizal zone, mol / cm2d + + + for r in range(1,len(simtimes)): + + dt = simtimes[r] + if r>0: + dt = simtimes[r] - simtimes[r-1] + + print('time',simtimes[r]) + print('no flux BC') + print("*****", "#", r, "external time step", dt, " d, simulation time", s_nf.simTime, "d, internal time step", s_nf.ddt, "d") + print('advective flow') + print("*****", "#", r, "external time step", dt, " d, simulation time", s_af.simTime, "d, internal time step", s_af.ddt, "d") + print('Dirichlet BC') + print("*****", "#", r, "external time step", dt, " d, simulation time", s_d.simTime, "d, internal time step", s_d.ddt, "d") + + #one timestep + s_nf.solve(dt, saveInnerFluxes_ = True) + s_af.solve(dt, saveInnerFluxes_ = True) + s_d.solve(dt, saveInnerFluxes_ = True) #TODO: this has numerical problems + + #watercontent and solute content, discretised + watercontent_nf = s_nf.getWaterContent() # cm3 + waterpotential_nf = np.array([vg.pressure_head(watercontent_nf[i],peri.sp) for i in range(NC)]) # cm + watercontent_af = s_af.getWaterContent() # cm3 + waterpotential_af = np.array([vg.pressure_head(watercontent_af[i],peri.sp) for i in range(NC)]) # cm + watercontent_d = s_d.getWaterContent() # cm3 + waterpotential_d = np.array([vg.pressure_head(watercontent_d[i],peri.sp) for i in range(NC)]) # cm + solutecontents_nf = s_nf.getSolution(1) / molarMassSolute # mol/cm3 + solutecontents_af = s_af.getSolution(1) / molarMassSolute # mol/cm3 + solutecontents_d = s_d.getSolution(1) / molarMassSolute # mol/cm3 + + #inflow and outflow + rootuptake_w_nf = s_nf.getInnerFlow(0, length) /(length*(2*np.pi*r_root)) # cm /d + rootuptake_s_nf = s_nf.getInnerFlow(1, length) / molarMassSolute /(length*(2*np.pi*r_root)) # mol / (cm2d) + inflow_w_nf = s_nf.getOuterFlow(0, length) /(length*(2*np.pi*r_prhiz)) # cm /d + inflow_s_nf = s_nf.getOuterFlow(1, length) / molarMassSolute /(length*(2*np.pi*r_prhiz)) # mol / (cm2d) + + rootuptake_w_af = s_af.getInnerFlow(0, length) /(length*(2*np.pi*r_root)) # cm /d + rootuptake_s_af = s_af.getInnerFlow(1, length) / molarMassSolute /(length*(2*np.pi*r_root)) # mol / (cm2d) + inflow_w_af = s_af.getOuterFlow(0, length) /(length*(2*np.pi*r_prhiz)) # cm /d + inflow_s_af = s_af.getOuterFlow(1, length) / molarMassSolute /(length*(2*np.pi*r_prhiz)) # mol / (cm2d) + + rootuptake_w_d = s_d.getInnerFlow(0, length) /(length*(2*np.pi*r_root)) # cm /d + rootuptake_s_d = s_d.getInnerFlow(1, length) / molarMassSolute /(length*(2*np.pi*r_root)) # mol / (cm2d) + inflow_w_d = s_d.getOuterFlow(0, length) /(length*(2*np.pi*r_prhiz)) # cm /d + inflow_s_d = s_d.getOuterFlow(1, length) / molarMassSolute /(length*(2*np.pi*r_prhiz)) # mol / (cm2d) + + rootuptake_w_nf = rootuptake_w_nf[0] + rootuptake_s_nf = rootuptake_s_nf[0] + inflow_w_nf = abs(inflow_w_nf[0]) + inflow_s_nf = inflow_s_nf[0] + rootuptake_w_af = rootuptake_w_af[0] + rootuptake_s_af = rootuptake_s_af[0] + inflow_w_af = abs(inflow_w_af[0]) + inflow_s_af = inflow_s_af[0] + rootuptake_w_d = rootuptake_w_d[0] + rootuptake_s_d = rootuptake_s_d[0] + inflow_w_d = abs(inflow_w_d[0]) + inflow_s_d = inflow_s_d[0] + + #store the dumux outputs + #scenario no flux outer BC + watercontent_dumux[0,r,0]=rootuptake_w_nf + watercontent_dumux[0,r,1]=inflow_w_nf + watercontent_dumux[0,r,2:]=watercontent_nf + waterpotential_dumux[0,r,0]=rootuptake_w_nf + waterpotential_dumux[0,r,1]=inflow_w_nf + waterpotential_dumux[0,r,2:]=waterpotential_nf + solutes_dumux[0,r,0]=rootuptake_s_nf + solutes_dumux[0,r,1]=inflow_s_nf + solutes_dumux[0,r,2:]=solutecontents_nf + #advective flow scenario + watercontent_dumux[1,r,0]=rootuptake_w_af + watercontent_dumux[1,r,1]=inflow_w_af + watercontent_dumux[1,r,2:]=watercontent_af + waterpotential_dumux[1,r,0]=rootuptake_w_af + waterpotential_dumux[1,r,1]=inflow_w_af + waterpotential_dumux[1,r,2:]=waterpotential_af + solutes_dumux[1,r,0]=rootuptake_s_af + solutes_dumux[1,r,1]=inflow_s_af + solutes_dumux[1,r,2:]=solutecontents_af + #Dirichlet BC + watercontent_dumux[2,r,0]=rootuptake_w_d + watercontent_dumux[2,r,1]=inflow_w_d + watercontent_dumux[2,r,2:]=watercontent_d + waterpotential_dumux[2,r,0]=rootuptake_w_d + waterpotential_dumux[2,r,1]=inflow_w_d + waterpotential_dumux[2,r,2:]=waterpotential_d + solutes_dumux[2,r,0]=rootuptake_s_d + solutes_dumux[2,r,1]=inflow_s_d + solutes_dumux[2,r,2:]=solutecontents_d + + #determine means + #for the explicit Euler method of the timesteps + f_root = (2*np.pi*r_root)/(np.pi*(r_prhiz**2-r_root**2)) + f_prhiz = Ds* pow(0.25,13/3)/(0.4**2)*(2*np.pi*r_prhiz)/(np.pi*(r_prhiz**2-r_root**2)) #only have diffusion on the outside for now, assume water content of 0.25, saturated watercontent of 0.4 + IC = initial_soluteconcentration + #no flux outer BC + mean_watercontent_nf = np.average(watercontent_nf, weights=volumes) + mean_waterpotential_nf = vg.pressure_head(mean_watercontent_nf, peri.sp) + mean_soluteconcent_nf = np.average(solutecontents_nf, weights=np.multiply(mean_watercontent_nf, volumes)) + mean_soluteconcent_sr_simp_nf = solutes_sr_simp[0,r-1,1]+(f_root*solutes_sr_simp[0,r-1,0])*dt + solutes_sr_simp[0,r,1] = mean_soluteconcent_sr_simp_nf + #advective flow + mean_watercontent_af = np.average(watercontent_af, weights=volumes) + mean_waterpotential_af = vg.pressure_head(mean_watercontent_af, peri.sp) + mean_soluteconcent_af = np.average(solutecontents_af, weights=np.multiply(mean_watercontent_af, volumes)) + mean_soluteconcent_sr_simp_af = solutes_sr_simp[1,r-1,1]+(f_root*solutes_sr_simp[1,r-1,0])*dt + solutes_sr_simp[1,r,1] = mean_soluteconcent_sr_simp_af + #Dirichlet BC + mean_watercontent_d = np.average(watercontent_d, weights=volumes) + mean_waterpotential_d = vg.pressure_head(mean_watercontent_d, peri.sp) + mean_soluteconcent_d = np.average(solutecontents_d, weights=np.multiply(mean_watercontent_d, volumes)) + mean_soluteconcent_sr_simp_d = solutes_sr_simp[2,r-1,1]+(f_root*solutes_sr_simp[2,r-1,0]+f_prhiz*(IC-solutes_sr_simp[2,r-1,-1]))*dt + solutes_sr_simp[2,r,1] = mean_soluteconcent_sr_simp_d + #determine coefficients for the analytical approximation + #equation [4] in Schroeder2008 doi:10.2136/vzj2007.0114 + #Phi(r)=Phi_outer + (q_root*r_root-q_out*r_prhiz)*(rho**2)/(1-rho**2)*(((r/r_prhiz)**2-1)/2-ln(r/r_prhiz))+q_out*r_prhiz*ln(r/r_prhiz) + #Assume Phi(0.53r)=Phi(mean water potential) + #both q_root and q_out are assumed to have positive signs if the water flows to the root + #no flux outer BC + Phi_soil_nf = vg.fast_mfp[peri.sp](mean_waterpotential_nf) #mfp of the mean soil + Phi_outer_nf = Phi_soil_nf - (rootuptake_w_nf*r_root-inflow_w_nf*r_prhiz)*(rho**2)/(1-rho**2)*(((0.53)**2-1)/2-np.log(0.53))-(inflow_w_nf*r_prhiz)*r_prhiz*np.log(0.53) #mfp at the outer perirhizal radius + Phi_nf = lambda r: Phi_outer_nf + (rootuptake_w_nf*r_root-inflow_w_nf*r_prhiz)*(rho**2)/(1-rho**2)*(((r/r_prhiz)**2-1)/2-np.log(r/r_prhiz))+(inflow_w_nf*r_prhiz)*r_prhiz*np.log(r/r_prhiz)#mfp function depending on radius + Phi_root_nf = Phi_nf(r_root)#mfp next to the root + #advective flow, Dirichlet for water transport + Phi_soil_af = vg.fast_mfp[peri.sp](mean_waterpotential_af) #mfp of the mean soil + Phi_outer_af = Phi_soil_af - (rootuptake_w_af*r_root-inflow_w_af*r_prhiz)*(rho**2)/(1-rho**2)*(((0.53)**2-1)/2-np.log(0.53))-(inflow_w_af*r_prhiz)*np.log(0.53) #mfp at the outer perirhizal radius + Phi_af = lambda r: Phi_outer_af + (rootuptake_w_af*r_root-inflow_w_af*r_prhiz)*(rho**2)/(1-rho**2)*(((r/r_prhiz)**2-1)/2-np.log(r/r_prhiz))+(inflow_w_af*r_prhiz)*np.log(r/r_prhiz)#mfp function depending on radius + Phi_root_af = Phi_af(r_root)#mfp next to the root + #Dirichlet BC + Phi_soil_d = vg.fast_mfp[peri.sp](mean_waterpotential_d) #mfp of the mean soil + Phi_outer_d = Phi_soil_d - (rootuptake_w_d*r_root-inflow_w_d*r_prhiz)*(rho**2)/(1-rho**2)*(((0.53)**2-1)/2-np.log(0.53))-(inflow_w_d*r_prhiz)*np.log(0.53) #mfp at the outer perirhizal radius + Phi_d = lambda r: Phi_outer_d + (rootuptake_w_d*r_root-inflow_w_d*r_prhiz)*(rho**2)/(1-rho**2)*(((r/r_prhiz)**2-1)/2-np.log(r/r_prhiz))+(inflow_w_d*r_prhiz)*np.log(r/r_prhiz)#mfp function depending on radius + Phi_root_d = Phi_d(r_root)#mfp next to the root + + #write the steady rate approximations for the water as outputs + waterpotential_sr[0,r,0] = rootuptake_w_nf + waterpotential_sr[0,r,1] = inflow_w_nf + waterpotential_sr[0,r,2:] = np.array([vg.fast_imfp[peri.sp](Phi_nf(CC[i])) for i in range(NC)]) + watercontent_sr[0,r,0] = rootuptake_w_nf + watercontent_sr[0,r,1] = inflow_w_nf + watercontent_sr[0,r,2:] = np.array([vg.water_content(waterpotential_sr[0,r,2+i],peri.sp) for i in range(NC)]) + + waterpotential_sr[1,r,0] = rootuptake_w_af + waterpotential_sr[1,r,1] = inflow_w_af + waterpotential_sr[1,r,2:] = np.array([vg.fast_imfp[peri.sp](Phi_af(CC[i])) for i in range(NC)]) + watercontent_sr[1,r,0] = rootuptake_w_af + watercontent_sr[1,r,1] = inflow_w_af + watercontent_sr[1,r,2:] = np.array([vg.water_content(waterpotential_sr[1,r,2+i],peri.sp) for i in range(NC)]) + + waterpotential_sr[2,r,0] = rootuptake_w_d + waterpotential_sr[2,r,1] = inflow_w_d + waterpotential_sr[2,r,2:] = np.array([vg.fast_imfp[peri.sp](Phi_d(CC[i])) for i in range(NC)]) + watercontent_sr[2,r,0] = rootuptake_w_d + watercontent_sr[2,r,1] = inflow_w_d + watercontent_sr[2,r,2:] = np.array([vg.water_content(waterpotential_sr[2,r,2+i],peri.sp) for i in range(NC)]) + + + + #compute the analytical approximations for the solute uptake + #case of dumux no flux outer BC + result_solutes_sr_nf = peri.soil_root_solutes_steadyrate_simplified_([Phi_root_nf], [Phi_soil_nf], [r_root], [r_prhiz], [mean_soluteconcent_sr_simp_nf], [Vmax_per_area], [Km], Ds, [waterdemand], peri.sp, n_approx = 5) + result_solutes_sr_af = peri.soil_root_solutes_steadyrate_simplified_([Phi_root_af], [Phi_soil_af], [r_root], [r_prhiz], [mean_soluteconcent_sr_simp_af], [Vmax_per_area], [Km], Ds, [waterdemand], peri.sp, n_approx = 5) + result_solutes_sr_d = peri.soil_root_solutes_steadyrate_simplified_([Phi_root_d], [Phi_soil_d], [r_root], [r_prhiz], [mean_soluteconcent_sr_simp_d], [Vmax_per_area], [Km], Ds, [waterdemand], peri.sp, n_approx = 5) + + #safe the results + #(here there is no computed inflow, so index 1 doesn't do anything) + result_solutes_sr_nf = result_solutes_sr_nf[0] + solutes_sr_simp[0,r,0]=-Vmax_per_area * result_solutes_sr_nf / (Km + result_solutes_sr_nf) + result_solutes_sr_af = result_solutes_sr_af[0] + solutes_sr_simp[1,r,0]=-Vmax_per_area * result_solutes_sr_af / (Km + result_solutes_sr_af) + result_solutes_sr_d = result_solutes_sr_d[0] + solutes_sr_simp[2,r,0]=-Vmax_per_area * result_solutes_sr_d / (Km + result_solutes_sr_d) + + F0_nf = peri.integral_AdvectionDiffusion_(Phi_nf(CC[0]),peri.sp) + F0_af = peri.integral_AdvectionDiffusion_(Phi_af(CC[0]),peri.sp) + F0_d = peri.integral_AdvectionDiffusion_(Phi_d(CC[0]),peri.sp) + D_tilde = 1/Ds/math.pow(sp.theta_S-sp.theta_R,13/3)*(sp.theta_S*sp.theta_S) + + # the ratio of waterflow and soluteflow is assumed to remain constant throughout the perirhizal zone + for j in range(NC): + r_current = CC[j] + F_nf = peri.integral_AdvectionDiffusion_(Phi_nf(r_current),peri.sp)-F0_nf + F_af = peri.integral_AdvectionDiffusion_(Phi_af(r_current),peri.sp)-F0_af + F_d = peri.integral_AdvectionDiffusion_(Phi_d(r_current),peri.sp)-F0_d + F_tilde_nf=math.exp(-D_tilde*F_nf) + F_tilde_af=math.exp(-D_tilde*F_af) + F_tilde_d=math.exp(-D_tilde*F_d) + solutes_sr_simp[0,r,2+j] = result_solutes_sr_nf / F_tilde_nf - (1-1/F_tilde_nf) * solutes_sr_simp[0,r,0] / (waterdemand) #rewrite this using water solutes disc + solutes_sr_simp[1,r,2+j] = result_solutes_sr_af / F_tilde_af - (1-1/F_tilde_af) * solutes_sr_simp[1,r,0] / (waterdemand) + solutes_sr_simp[2,r,2+j] = result_solutes_sr_d / F_tilde_d - (1-1/F_tilde_d) * solutes_sr_simp[2,r,0] / (waterdemand) + + #Tiina Roose approximation + E = 0 #minimal solute uptake + DS_TR = Ds * math.pow(mean_watercontent_nf,10/3)/(sp.theta_S**2) + waterflow_TR = 2*np.pi*r_root*waterdemand + waterflow_TR = waterdemand + rsc = peri.solutesuptake_convdiff_([mean_watercontent_nf], [initial_soluteconcentration], [Vmax_per_area], [Km], DS_TR, [waterflow_TR], [r_root], [E], [simtimes[r]], sp) + solutes_TR[0,r,0]=-Vmax_per_area*rsc[0]/(Km + rsc[0]) + solutes_TR[0,r,2]=rsc[0] + + DS_TR = Ds * math.pow(mean_watercontent_af,10/3)/(sp.theta_S**2) + rsc = peri.solutesuptake_convdiff_([mean_watercontent_af], [initial_soluteconcentration], [Vmax_per_area], [Km], DS_TR, [waterflow_TR], [r_root], [E], [simtimes[r]], sp) + solutes_TR[1,r,0]=-Vmax_per_area*rsc[0]/(Km + rsc[0]) + solutes_TR[1,r,2]=rsc[0] + + DS_TR = Ds * math.pow(mean_watercontent_d,10/3)/(sp.theta_S**2) + rsc = peri.solutesuptake_convdiff_([mean_watercontent_d], [initial_soluteconcentration], [Vmax_per_area], [Km], DS_TR, [waterflow_TR], [r_root], [E], [simtimes[r]], sp) + solutes_TR[2,r,0]=-Vmax_per_area*rsc[0]/(Km + rsc[0]) + solutes_TR[2,r,2]=rsc[0] + + # case of general steady rate water uptake + # safe the means, they are computed via the explicit Euler timestepping scheme + + mean_soluteconcent_ss_af = solutes_ss[1,r-1,1]+f_root*solutes_ss[1,r-1,0]*dt + mean_soluteconcent_ss_d = solutes_ss[2,r-1,1]+(f_root*solutes_ss[2,r-1,0]+f_prhiz*(IC-solutes_ss[2,r-1,-1]))*dt + mean_soluteconcent_sr_nf = solutes_sr[0,r-1,1]+f_root*solutes_sr[0,r-1,0]*dt + mean_soluteconcent_sr_af = solutes_sr[1,r-1,1]+f_root*solutes_sr[1,r-1,0]*dt + mean_soluteconcent_sr_d = solutes_sr[2,r-1,1]+f_root*solutes_sr[2,r-1,0]*dt + mean_soluteconcent_d_d = solutes_d[2,r-1,1]+(f_root*solutes_d[2,r-1,0]+f_prhiz*(IC-solutes_d[2,r-1,-1]))*dt + mean_soluteconcent_u_nf = solutes_u[0,r-1,1]+f_root*solutes_u[0,r-1,0]*dt + mean_soluteconcent_u_af = solutes_u[1,r-1,1]+f_root*solutes_u[1,r-1,0]*dt + mean_soluteconcent_u_d = solutes_u[2,r-1,1]+(f_root*solutes_u[2,r-1,0]+f_prhiz*(IC-solutes_u[2,r-1,-1]))*dt + + solutes_ss[1,r,1] = mean_soluteconcent_ss_af + solutes_ss[2,r,1] = mean_soluteconcent_ss_d + solutes_sr[0,r,1] = mean_soluteconcent_sr_nf + solutes_sr[1,r,1] = mean_soluteconcent_sr_af + solutes_sr[2,r,1] = mean_soluteconcent_sr_d + solutes_d[2,r,1] = mean_soluteconcent_d_d + solutes_u[0,r,1] = mean_soluteconcent_u_nf + solutes_u[1,r,1] = mean_soluteconcent_u_af + solutes_u[2,r,1] = mean_soluteconcent_u_d + + #for the steady state take again the outer concentration + #general steady state + rsc, Uptake, quadratic_flow, c_noflux = peri.soil_root_solutes_sr([Phi_outer_af], [rootuptake_w_af*2*np.pi*r_root], [inflow_w_af*2*np.pi*r_prhiz], [CC[0]], [CC[-1]], [mean_soluteconcent_af], [initial_soluteconcentration], [Vmax], [Km], [Ds], peri.sp, mode = "ss") + _, _, soluteconcentration, soluteconcentration_mean = peri.watersolutes_disc(Phi_outer_af, 2*np.pi * rootuptake_w_af*r_root, 2*np.pi * inflow_w_af*r_prhiz, CC[0], CC[-1], CC, Ds, Uptake, quadratic_flow, c_noflux, peri.sp) + solutes_ss[1,r,0] = -(Uptake[0]+r_root**2*quadratic_flow[0]) / (2 * np.pi * r_root) + #solutes_ss[1,r,1] = -(Uptake[0] + r_prhiz**2 * quadratic_flow[0]) + solutes_ss[1,r,2:] = soluteconcentration[:] + rsc, Uptake, quadratic_flow, c_noflux = peri.soil_root_solutes_sr([Phi_outer_d], [rootuptake_w_d*2*np.pi*r_root], [inflow_w_d*2*np.pi*r_prhiz], [CC[0]], [CC[-1]], [mean_soluteconcent_d], [initial_soluteconcentration], [Vmax], [Km], [Ds], peri.sp, mode = "ss") + _, _, soluteconcentration, soluteconcentration_mean = peri.watersolutes_disc(Phi_outer_d, 2*np.pi * rootuptake_w_d*r_root, 2*np.pi * inflow_w_d*r_prhiz, CC[0], CC[-1], CC, Ds, Uptake, quadratic_flow, c_noflux, peri.sp) + solutes_ss[2,r,0] = -(Uptake[0]+r_root**2*quadratic_flow[0]) / (2 * np.pi * r_root) + #solutes_ss[2,r,1] = -(Uptake[0] + r_prhiz**2 * quadratic_flow[0]) + solutes_ss[2,r,2:] = soluteconcentration[:] + #general steady rate no flux outer BC + rsc, Uptake, quadratic_flow, c_noflux = peri.soil_root_solutes_sr([Phi_outer_nf], [rootuptake_w_nf*2*np.pi*r_root], [inflow_w_nf*2*np.pi*r_prhiz], [CC[0]], [CC[-1]], [mean_soluteconcent_nf], [initial_soluteconcentration], [Vmax], [Km], [Ds], peri.sp, mode = "sr") + _, _, soluteconcentration, soluteconcentration_mean = peri.watersolutes_disc(Phi_outer_nf, 2*np.pi * rootuptake_w_nf*r_root, 2*np.pi * inflow_w_nf*r_prhiz, CC[0], CC[-1], CC, Ds, Uptake, quadratic_flow, c_noflux, peri.sp) + #print("rsc_sr_nf", rsc, Uptake, quadratic_flow, c_noflux) + solutes_sr[0,r,0] = -(Uptake[0]+r_root**2*quadratic_flow[0]) / (2 * np.pi * r_root) + #solutes_sr[0,r,1] = -(Uptake[0] + r_prhiz**2 * quadratic_flow[0]) + solutes_sr[0,r,2:] = soluteconcentration[:] + rsc, Uptake, quadratic_flow, c_noflux = peri.soil_root_solutes_sr([Phi_outer_af], [rootuptake_w_af*2*np.pi*r_root], [inflow_w_af*2*np.pi*r_prhiz], [CC[0]], [CC[-1]], [mean_soluteconcent_af], [initial_soluteconcentration], [Vmax], [Km], [Ds], peri.sp, mode = "sr") + _, _, soluteconcentration, soluteconcentration_mean = peri.watersolutes_disc(Phi_outer_af, 2*np.pi * rootuptake_w_af*r_root, 2*np.pi * inflow_w_af*r_prhiz, CC[0], CC[-1], CC, Ds, Uptake, quadratic_flow, c_noflux, peri.sp) + #print("rsc_sr_af", rsc, Uptake, quadratic_flow, c_noflux) + solutes_sr[1,r,0] = -(Uptake[0]+r_root**2*quadratic_flow[0]) / (2 * np.pi * r_root) + #solutes_sr[1,r,1] = -(Uptake[0] + r_prhiz**2 * quadratic_flow[0]) + solutes_sr[1,r,2:] = soluteconcentration[:] + rsc, Uptake, quadratic_flow, c_noflux = peri.soil_root_solutes_sr([Phi_outer_d], [rootuptake_w_d*2*np.pi*r_root], [inflow_w_d*2*np.pi*r_prhiz], [CC[0]], [CC[-1]], [mean_soluteconcent_d], [initial_soluteconcentration], [Vmax], [Km], [Ds], peri.sp, mode = "sr") + _, _, soluteconcentration, soluteconcentration_mean = peri.watersolutes_disc(Phi_outer_d, 2*np.pi * rootuptake_w_d*r_root, 2*np.pi * inflow_w_d*r_prhiz, CC[0], CC[-1], CC, Ds, Uptake, quadratic_flow, c_noflux, peri.sp) + #print("rsc_sr_d", rsc, Uptake, quadratic_flow, c_noflux) + solutes_sr[2,r,0] = -(Uptake[0]+r_root**2*quadratic_flow[0]) / (2 * np.pi * r_root) + #solutes_sr[2,r,1] = -(Uptake[0] + r_prhiz**2 * quadratic_flow[0]) + solutes_sr[2,r,2:] = soluteconcentration[:] + #Dirichlet BC + rsc, Uptake, quadratic_flow, c_noflux = peri.soil_root_solutes_sr([Phi_outer_d], [rootuptake_w_d*2*np.pi*r_root], [inflow_w_d*2*np.pi*r_prhiz], [CC[0]], [CC[-1]], [mean_soluteconcent_d], [initial_soluteconcentration], [Vmax], [Km], [Ds], peri.sp, mode = "dirichlet") + _, _, soluteconcentration, soluteconcentration_mean = peri.watersolutes_disc(Phi_outer_d, 2*np.pi * rootuptake_w_d*r_root, 2*np.pi * inflow_w_d*r_prhiz, CC[0], CC[-1], CC, Ds, Uptake, quadratic_flow, c_noflux, peri.sp) + solutes_d[2,r,0] = -(Uptake[0]+r_root**2*quadratic_flow[0]) / (2 * np.pi * r_root) + #solutes_d[2,r,1] = -(Uptake[0] + r_prhiz**2 * quadratic_flow[0]) + solutes_d[2,r,2:] = soluteconcentration[:] + #uniform concentration + solutes_u[0,r,0] = -Vmax_per_area*mean_soluteconcent_u_nf/(Km + mean_soluteconcent_u_nf) + solutes_u[1,r,0] = -Vmax_per_area*mean_soluteconcent_u_af/(Km + mean_soluteconcent_u_af) + solutes_u[2,r,0] = -Vmax_per_area*mean_soluteconcent_u_d/(Km + mean_soluteconcent_u_d) + solutes_u[0,r,2:] = np.ones(NC) * mean_soluteconcent_u_nf + solutes_u[1,r,2:] = np.ones(NC) * mean_soluteconcent_u_af + solutes_u[2,r,2:] = np.ones(NC) * mean_soluteconcent_u_d + + + return watercontent_dumux, waterpotential_dumux, watercontent_sr, waterpotential_sr, solutes_dumux, solutes_sr_simp, solutes_ss, solutes_sr, solutes_d, solutes_u, solutes_TR + +if do_computation: + # save everything in the np arrays + for i in range(n_tests): + watercontent_dumux[i,:,:,:], waterpotential_dumux[i,:,:,:], watercontent_sr[i,:,:,:], waterpotential_sr[i,:,:,:], solutes_dumux[i,:,:,:], solutes_sr_simp[i,:,:,:], solutes_ss[i,:,:,:], solutes_sr[i,:,:,:], solutes_d[i,:,:,:], solutes_u[i,:,:,:], solutes_TR[i,:,:,:] = run_perirhizal_test(max_time, n_times, r_prhiz, r_root, NC, points, CC, volumes, length, n_scenarios, initial_waterpotential, initial_soluteconcentration) + np.savez("test_perirhizal.npz", + watercontent_dumux=watercontent_dumux, + waterpotential_dumux=waterpotential_dumux, + watercontent_sr=watercontent_sr, + waterpotential_sr=waterpotential_sr, + solutes_dumux=solutes_dumux, + solutes_sr_simp=solutes_sr_simp, + solutes_ss=solutes_ss, + solutes_sr=solutes_sr, + solutes_d=solutes_d, + solutes_u=solutes_u, + solutes_TR=solutes_TR) + +else: + simulation_results = np.load("test_perirhizal.npz") + watercontent_dumux = simulation_results["watercontent_dumux"] + waterpotential_dumux = simulation_results["waterpotential_dumux"] + watercontent_sr = simulation_results["watercontent_sr"] + waterpotential_sr = simulation_results["waterpotential_sr"] + solutes_dumux = simulation_results["solutes_dumux"] + solutes_sr_simp = simulation_results["solutes_sr_simp"] + solutes_ss = simulation_results["solutes_ss"] + solutes_sr = simulation_results["solutes_sr"] + solutes_d = simulation_results["solutes_d"] + solutes_u = simulation_results["solutes_u"] + solutes_TR = simulation_results["solutes_TR"] + + + + +# compare both for the differint means of water / solute content +run = 0 +timestep = np.array(np.linspace(1,9,num=5)) +for i in range(5): + timestep[i] = int(n_times * timestep[i] / 10) +timestep = timestep.astype(int) + + +linestyle_dumux = "solid" +linestyle_steadystate = "dotted" +linestyle_steadyrate = "dashed" +linestyle_special = "dashdot" + + +fig, ax1 = figure_style.subplots12(nrows=5, ncols=3) +# dumux(all 3) +# left: sr no flux +# middle: dirichlet BC water, advective flow solutes +# right: ss, sr, farfield for sr waterflow + + + +for i in range(5): + ax2_0 = ax1[i,0].twinx() + ax2_1 = ax1[i,1].twinx() + ax2_2 = ax1[i,2].twinx() + + #load data + water_dumux_nf = waterpotential_dumux[run, 0, timestep[i], 2:] + water_dumux_af = waterpotential_dumux[run, 1, timestep[i], 2:] + water_dumux_d = waterpotential_dumux[run, 2, timestep[i], 2:] + water_steadyrate_nf = waterpotential_sr[run, 0, timestep[i], 2:] + water_steadyrate_af = waterpotential_sr[run, 1, timestep[i], 2:] + water_steadyrate_d = waterpotential_sr[run, 2, timestep[i], 2:] + solute_dumux_nf = solutes_dumux[run, 0, timestep[i], 2:] + solute_dumux_af = solutes_dumux[run, 1, timestep[i], 2:] + solute_dumux_d = solutes_dumux[run, 2, timestep[i], 2:] + + solutes_sr_simp_nf = solutes_sr_simp[run, 0, timestep[i], 2:] + solutes_sr_simp_af = solutes_sr_simp[run, 1, timestep[i], 2:] + solutes_sr_simp_d = solutes_sr_simp[run, 2, timestep[i], 2:] + + solutes_ss_af = solutes_ss[run, 1, timestep[i], 2:] + solutes_ss_d = solutes_ss[run, 2, timestep[i], 2:] + + solutes_sr_nf = solutes_sr[run, 0, timestep[i], 2:] + solutes_sr_af = solutes_sr[run, 1, timestep[i], 2:] + solutes_sr_d = solutes_sr[run, 2, timestep[i], 2:] + + solutes_d_d = solutes_d[run, 2, timestep[i], 2:] + + solutes_u_nf = solutes_u[run, 0, timestep[i], 2:] + solutes_u_af = solutes_u[run, 1, timestep[i], 2:] + solutes_u_d = solutes_u[run, 2, timestep[i], 2:] + + + #left plot: no flux outer BC + ax1[i,0].plot(CC, water_dumux_nf, "b", linestyle = linestyle_dumux, label = "water_dumux") + ax1[i,0].plot(CC, water_steadyrate_nf, "b", linestyle = linestyle_steadyrate, label = "water_sr") + ax2_0.plot(CC, solute_dumux_nf, "m", linestyle = linestyle_dumux, label = "solute_dumux") + ax2_0.plot(CC, solutes_sr_simp_nf, "m", linestyle = linestyle_steadyrate, label = "solute_sr_simp") + ax2_0.plot(CC, solutes_sr_nf, "g", linestyle = linestyle_steadyrate, label = "solute_sr") + ax2_0.plot(CC, solutes_u_nf, "r", linestyle = linestyle_dumux, label = "solute_u") + + #middle plot: water dirichlet, advective flow solutes + ax1[i,1].plot(CC, water_dumux_af, "b", linestyle = linestyle_dumux, label = "water_dumux") + ax1[i,1].plot(CC, water_steadyrate_af, "b", linestyle = linestyle_steadyrate, label = "water_sr") + ax2_1.plot(CC, solute_dumux_af, "m", linestyle = linestyle_dumux, label = "solute_dumux") + ax2_1.plot(CC, solutes_sr_simp_af, "m", linestyle = linestyle_steadyrate, label = "solute_sr_simp") + ax2_1.plot(CC, solutes_sr_af, "g", linestyle = linestyle_steadyrate, label = "solute_sr") + ax2_1.plot(CC, solutes_ss_af, "y", linestyle = linestyle_steadystate, label = "solute_ss") + ax2_1.plot(CC, solutes_u_af, "r", linestyle = linestyle_dumux, label = "solute_u") + + #right plot: Dirichlet (initial conditions) outer BC + ax1[i,2].plot(CC, water_dumux_d, "b", linestyle = linestyle_dumux, label = "water_dumux") + ax1[i,2].plot(CC, water_steadyrate_d, "b", linestyle = linestyle_steadyrate, label = "water_sr") + ax2_2.plot(CC, solute_dumux_d, "m", linestyle = linestyle_dumux, label = "solute_dumux") + ax2_2.plot(CC, solutes_sr_simp_d, "m", linestyle = linestyle_steadyrate, label = "solute_sr_simp") + ax2_2.plot(CC, solutes_sr_d, "g", linestyle = linestyle_steadyrate, label = "solute_sr") + ax2_2.plot(CC, solutes_ss_d, "y", linestyle = linestyle_steadystate, label = "solute_ss") + ax2_2.plot(CC, solutes_d_d, "c", linestyle = linestyle_special, label = "solute_d") + ax2_2.plot(CC, solutes_u_d, "r", linestyle = linestyle_dumux, label = "solute_u") + +ax1[i,0].set_xlabel("distance root [cm]") +ax1[i,0].set_ylabel("water") +ax2_0.set_ylabel("nitrogen") +ax1[i,0].legend(["watercontent cm3/cm3"], loc="upper left") +ax2_0.legend(["nitrogen concentration mol/cm3"], loc="upper right") + +ax1[i,0].legend(loc="upper left") +ax2_0.legend(loc="upper right") + +ax1[i,1].set_xlabel("distance root [cm]") +ax1[i,1].set_ylabel("water") +ax2_1.set_ylabel("nitrogen") +ax1[i,1].legend(["watercontent cm3/cm3"], loc="upper left") +ax2_1.legend(["nitrogen concentration mol/cm3"], loc="upper right") + +ax1[i,1].legend(loc="upper left") +ax2_1.legend(loc="upper right") + +ax1[i,2].set_xlabel("distance root [cm]") +ax1[i,2].set_ylabel("water") +ax2_2.set_ylabel("nitrogen") +ax1[i,2].legend(["watercontent cm3/cm3"], loc="upper left") +ax2_2.legend(["nitrogen concentration mol/cm3"], loc="upper right") + +ax1[i,2].legend(loc="upper left") +ax2_2.legend(loc="upper right") + +#np.save("input/" + filename + "_fp", np.vstack((sim_times_, -np.array(t_act_), np.array(q_soil_)))) +figure = plt.gcf() # get current figure +figure.set_size_inches(8, 6) +plt.savefig("concentrations", dpi = 100) +plt.show() + + + +fig, ax1 = figure_style.subplots12(nrows=1, ncols=3) +# dumux(both) +# left: sr no flux, ss for the no flux +# right: ss, sr, farfield for sr waterflow + +#load solute uptake +suptake_dumux_nf = solutes_dumux[run, 0, 1:, 0] +suptake_dumux_af = solutes_dumux[run, 1, 1:, 0] +suptake_dumux_d = solutes_dumux[run, 2, 1:, 0] + +suptake_sr_simp_nf = solutes_sr_simp[run, 0, 1:, 0] +suptake_sr_simp_af = solutes_sr_simp[run, 1, 1:, 0] +suptake_sr_simp_d = solutes_sr_simp[run, 2, 1:, 0] + +suptake_ss_af = solutes_ss[run, 1, 1:, 0] +suptake_ss_d = solutes_ss[run, 2, 1:, 0] + +suptake_sr_nf = solutes_sr[run, 0, 1:, 0] +suptake_sr_af = solutes_sr[run, 1, 1:, 0] +suptake_sr_d = solutes_sr[run, 2, 1:, 0] + +suptake_d_d = solutes_d[run, 2, 1:, 0] + +suptake_u_nf = solutes_u[run, 0, 1:, 0] +suptake_u_af = solutes_u[run, 1, 1:, 0] +suptake_u_d = solutes_u[run, 2, 1:, 0] + +suptake_TR_nf = solutes_TR[run, 0, 1:, 0] +suptake_TR_af = solutes_TR[run, 1, 1:, 0] +suptake_TR_d = solutes_TR[run, 2, 1:, 0] + + +ax1[0].plot(suptake_dumux_nf, suptake_dumux_nf, "m", linestyle = linestyle_dumux, label = "dumux") +#ax1[0].scatter(suptake_dumux_nf, abs(suptake_dumux_nf), "m", marker = "*") +ax1[0].plot(suptake_dumux_nf, abs(suptake_sr_simp_nf), "b", linestyle = linestyle_dumux, label = "steady rate simp") +ax1[0].plot(suptake_dumux_nf, abs(suptake_sr_nf), "g", linestyle = linestyle_dumux, label = "steady rate") +ax1[0].plot(suptake_dumux_nf, abs(suptake_u_nf), "r", linestyle = linestyle_dumux, label = "uniform") +ax1[0].plot(suptake_dumux_nf, abs(suptake_TR_nf), "k", linestyle = linestyle_dumux, label = "TR") + +ax1[1].plot(suptake_dumux_af, suptake_dumux_af, "m", linestyle = linestyle_dumux, label = "dumux") +#ax1[1].scatter(suptake_dumux_af, abs(suptake_dumux_af), marker = "*") +ax1[1].plot(suptake_dumux_af, abs(suptake_sr_simp_af), "b", linestyle = linestyle_dumux, label = "steady rate simp") +ax1[1].plot(suptake_dumux_af, abs(suptake_sr_af), "g", linestyle = linestyle_dumux, label = "steady rate") +ax1[1].plot(suptake_dumux_af, abs(suptake_ss_af), "y", linestyle = linestyle_dumux, label = "steady state") +ax1[1].plot(suptake_dumux_af, abs(suptake_u_af), "r", linestyle = linestyle_dumux, label = "uniform") +ax1[1].plot(suptake_dumux_af, abs(suptake_TR_af), "k", linestyle = linestyle_dumux, label = "TR") + +ax1[2].plot(suptake_dumux_d, suptake_dumux_d, "m", linestyle = linestyle_dumux, label = "dumux") +#ax1[2].scatter(suptake_dumux_d, abs(suptake_dumux_d), marker = "*") +ax1[2].plot(suptake_dumux_d, abs(suptake_sr_simp_d), "b", linestyle = linestyle_dumux, label = "steady rate simp") +ax1[2].plot(suptake_dumux_d, abs(suptake_sr_d), "g", linestyle = linestyle_dumux, label = "steady rate") +ax1[2].plot(suptake_dumux_d, abs(suptake_ss_d), "y", linestyle = linestyle_dumux, label = "steady state") +ax1[2].plot(suptake_dumux_d, abs(suptake_d_d), "c", linestyle = linestyle_dumux, label = "steady state") +ax1[2].plot(suptake_dumux_d, abs(suptake_u_d), "r", linestyle = linestyle_dumux, label = "uniform") +ax1[2].plot(suptake_dumux_d, abs(suptake_TR_d), "k", linestyle = linestyle_dumux, label = "TR") + +ax1[0].set_xlabel("dumux solute uptake") +ax1[0].set_ylabel("analytical approximation") +ax1[0].legend(["solute utake mol/cm2d"], loc="upper left") + +ax1[1].set_xlabel("dumux solute uptake") +ax1[1].set_ylabel("analytical approximation") +ax1[1].legend(["solute utake mol/cm2d"], loc="upper left") + +ax1[2].set_xlabel("dumux solute uptake") +ax1[2].set_ylabel("analytical approximation") +ax1[2].legend(["solute utake mol/cm2d"], loc="upper left") + +#np.save("input/" + filename + "_fp", np.vstack((sim_times_, -np.array(t_act_), np.array(q_soil_)))) +figure = plt.gcf() # get current figure +figure.set_size_inches(8, 6) +plt.savefig("relative_uptake", dpi = 100) +plt.show() + + +fig, ax1 = figure_style.subplots12(nrows=1, ncols=3) +# dumux(both) +# left: sr no flux, ss for the no flux +# right: ss, sr, farfield for sr waterflow + +simtimes = np.linspace(0,max_time,n_times) + +#load solute uptake +suptake_dumux_nf = solutes_dumux[run, 0, 1:, 0] +#doublecheck if the uptake is correct +beginning_solutes = initial_soluteconcentration * vg.water_content(initial_waterpotential,sp) * (np.pi * (r_prhiz**2 - r_root**2) * length) +total_uptake_nf = [np.sum([solutes_dumux[0, 0, i, 2+j] * watercontent_dumux[0, 0, i, 2+j] * np.pi * (points[j+1]**2-points[j]**2) * length for j in range(NC)]) for i in range(n_times)] +total_uptake_nf = np.array([(beginning_solutes - total_uptake_nf[i])/(2*np.pi*r_root*length)/max_time*(n_times-1) for i in range(n_times)]) +suptake_dumux_af = solutes_dumux[run, 1, 1:, 0] +suptake_dumux_d = solutes_dumux[run, 2, 1:, 0] + +suptake_sr_simp_nf = solutes_sr_simp[run, 0, 1:, 0] +suptake_sr_simp_af = solutes_sr_simp[run, 1, 1:, 0] +suptake_sr_simp_d = solutes_sr_simp[run, 2, 1:, 0] + +suptake_ss_af = solutes_ss[run, 1, 1:, 0] +suptake_ss_d = solutes_ss[run, 2, 1:, 0] + +suptake_sr_nf = solutes_sr[run, 0, 1:, 0] +suptake_sr_af = solutes_sr[run, 1, 1:, 0] +suptake_sr_d = solutes_sr[run, 2, 1:, 0] + +suptake_d_d = solutes_d[run, 2, 1:, 0] + +suptake_u_nf = solutes_u[run, 0, 1:, 0] +suptake_u_af = solutes_u[run, 1, 1:, 0] +suptake_u_d = solutes_u[run, 2, 1:, 0] + +suptake_dumux_nf = np.array([ sum(suptake_dumux_nf[:i]) for i in range(n_times)]) +suptake_dumux_af = np.array([ sum(suptake_dumux_af[:i]) for i in range(n_times)]) +suptake_dumux_d = np.array([ sum(suptake_dumux_d[:i]) for i in range(n_times)]) + +suptake_sr_simp_nf = np.array([ sum(suptake_sr_simp_nf[:i]) for i in range(n_times)]) +suptake_sr_simp_af = np.array([ sum(suptake_sr_simp_af[:i]) for i in range(n_times)]) +suptake_sr_simp_d = np.array([ sum(suptake_sr_simp_d[:i]) for i in range(n_times)]) + +suptake_ss_af = np.array([ sum(suptake_ss_af[:i]) for i in range(n_times)]) +suptake_ss_d = np.array([ sum(suptake_ss_d[:i]) for i in range(n_times)]) + +suptake_sr_nf = np.array([ sum(suptake_sr_nf[:i]) for i in range(n_times)]) +suptake_sr_af = np.array([ sum(suptake_sr_af[:i]) for i in range(n_times)]) +suptake_sr_d = np.array([ sum(suptake_sr_d[:i]) for i in range(n_times)]) + +suptake_d_d = np.array([ sum(suptake_d_d[:i]) for i in range(n_times)]) + +suptake_u_nf = np.array([ sum(suptake_u_nf[:i]) for i in range(n_times)]) +suptake_u_af = np.array([ sum(suptake_u_af[:i]) for i in range(n_times)]) +suptake_u_d = np.array([ sum(suptake_u_d[:i]) for i in range(n_times)]) + +suptake_TR_nf = np.array([ sum(suptake_TR_nf[:i]) for i in range(n_times)]) +suptake_TR_af = np.array([ sum(suptake_TR_af[:i]) for i in range(n_times)]) +suptake_TR_d = np.array([ sum(suptake_TR_d[:i]) for i in range(n_times)]) + + + +ax1[0].plot(simtimes, abs(suptake_dumux_nf), "m", linestyle = linestyle_dumux, label = "dumux") +#ax1[0].plot(simtimes[1:], abs(total_uptake_nf[1:]), "g", linestyle = linestyle_dumux, label = "dumux 2 uptake") #should produce the same line +ax1[0].plot(simtimes, abs(suptake_sr_simp_nf), "b", linestyle = linestyle_dumux, label = "steady rate simp") +ax1[0].plot(simtimes, abs(suptake_sr_nf), "g", linestyle = linestyle_dumux, label = "steady rate") + + +ax1[1].plot(simtimes, abs(suptake_dumux_af), "m", linestyle = linestyle_dumux, label = "dumux") +ax1[1].plot(simtimes, abs(suptake_sr_simp_af), "b", linestyle = linestyle_dumux, label = "steady rate simp") +ax1[1].plot(simtimes, abs(suptake_sr_af), "g", linestyle = linestyle_dumux, label = "steady rate") +ax1[1].plot(simtimes, abs(suptake_ss_af), "y", linestyle = linestyle_dumux, label = "steady state") + +ax1[2].plot(simtimes, abs(suptake_dumux_d), "m", linestyle = linestyle_dumux, label = "dumux") +ax1[2].plot(simtimes, abs(suptake_sr_simp_d), "b", linestyle = linestyle_dumux, label = "steady rate simp") +ax1[2].plot(simtimes, abs(suptake_sr_d), "g", linestyle = linestyle_dumux, label = "steady rate nf BC") +ax1[2].plot(simtimes, abs(suptake_ss_d), "y", linestyle = linestyle_dumux, label = "steady state") +ax1[2].plot(simtimes, abs(suptake_d_d), "c", linestyle = linestyle_dumux, label = "steady rate D BC") + +if showuniform: + ax1[0].plot(simtimes, abs(suptake_u_nf), "r", linestyle = linestyle_dumux, label = "uniform") + ax1[1].plot(simtimes, abs(suptake_u_af), "r", linestyle = linestyle_dumux, label = "uniform") + ax1[2].plot(simtimes, abs(suptake_u_d), "r", linestyle = linestyle_dumux, label = "uniform") + +if showTiina: + ax1[0].plot(simtimes, abs(suptake_TR_nf), "k", linestyle = linestyle_dumux, label = "TR") + ax1[1].plot(simtimes, abs(suptake_TR_af), "k", linestyle = linestyle_dumux, label = "TR") + ax1[2].plot(simtimes, abs(suptake_TR_d), "k", linestyle = linestyle_dumux, label = "TR") +#ax1[0].set_xlabel("dumux solute uptake") +#ax1[0].set_ylabel("analytical approximation") +#ax1[0].legend(["solute utake mol/cm2d"], loc="upper left") + +#ax1[1].set_xlabel("dumux solute uptake") +#ax1[1].set_ylabel("analytical approximation") +#ax1[1].legend(["solute utake mol/cm2d"], loc="upper left") + +#ax1[2].set_xlabel("dumux solute uptake") +#ax1[2].set_ylabel("analytical approximation") +#ax1[2].legend(["solute utake mol/cm2d"], loc="upper left") + +ax1[0].legend(loc="upper left") +ax1[1].legend(loc="upper left") +ax1[2].legend(loc="upper left") + +#np.save("input/" + filename + "_fp", np.vstack((sim_times_, -np.array(t_act_), np.array(q_soil_)))) +figure = plt.gcf() # get current figure +figure.set_size_inches(8, 6) +plt.savefig("cumulative_uptake", dpi = 100) +plt.show() diff --git a/src/functional/Perirhizal.py b/src/functional/Perirhizal.py index 382b11c1a..a235bb36f 100644 --- a/src/functional/Perirhizal.py +++ b/src/functional/Perirhizal.py @@ -1,8 +1,16 @@ import numpy as np +from numpy import linalg as LA import scipy.linalg as la from scipy.interpolate import RegularGridInterpolator from scipy.optimize import fsolve, root_scalar from scipy.spatial import ConvexHull, Voronoi +from scipy.integrate import odeint +from scipy import integrate +from scipy.special.orthogonal import ts_roots + +import time #only used for a test + +import math import plantbox as pb import plantbox.functional.van_genuchten as vg @@ -17,14 +25,14 @@ class PerirhizalPython(Perirhizal): * calculates root soil interface potential using steady rate approximation (Schröder et al. 2008) * perirhizal_conductance_per_layer calculates the perirhizal conductance in a soil layer (or cell) (Vanderborght et al. 2023, Eqn [6]) - * support of 4D lookup tables (file type is a zipped archive) -> will be IMPROVED + * support of 2D lookup tables (file type is a zipped archive) -> will be IMPROVED * * analysis across soil grid, e.g. get_density, average, aggregate * calculates outer perirhizal radii (based on denisties or Voronoi) * TODO maybe seperate geometry related functions (e.g. get_density, get_outer_radii) from soil root interface potential related functions (e.g. soil_root_interface_potentials, create_lookup) into different classes - run script for different examples, uncomment (to create a 4D lookup table, usage lookup table, and voronoi outer radii) + run script for different examples, uncomment (to create a 2D lookup table, usage lookup table, and voronoi outer radii) """ def __init__(self, ms=None): @@ -34,24 +42,81 @@ def __init__(self, ms=None): else: super().__init__() - self.lookup_table = None # optional 4d look up table to find soil root interface potentials - self.sp = None # corresponding van gencuchten soil parameter + self.lookup_table = None # optional 2d look up table to find soil root interface potentials + self.global_lookup_table = None # optional 3d look up table to find soil root interface potentials + self.lookup_table_sr_solutes = None # optional 1d look up table for steady state solute flow + self.lookup_table_sr_solutes_simplified = None # optional 2d lookup tables for the steady rate solute flow (one every diffusion coefficient) + self.sp = None # corresponding van genuchten soil parameter + + self.alpha_0 = 0.3 # a constant that is used for numerically solving the perirhizal waterflow + self.h_wilt = -16000 + self.Ds0_ref = 1 # reference diffusion coefficient of solutes in water [cm2/d] + self.r0_ref = 1.0e-3 #reference lower bound on the radius[cm] + self.water_filename = "lookup_perirhizal_waterflow_global" #name and location for the global lookup table def set_soil(self, sp): """sets VG parameters, and no look up table (slow)""" vg.create_mfp_lookup(sp) self.sp = sp self.lookup_table = None - + self.global_lookup_table = None + self.lookup_table_sr_solutes = None + self.lookup_table_sr_solutes_simplified = None + + def open_lookup(self, filename): - """opens a look-up table from a file to quickly find soil root interface potentials""" + """opens a look-up table from a file to quickly find soil root interface potentials""" npzfile = np.load(filename + ".npz") interface = npzfile["interface"] - rx_, sx_, akr_, rho_ = npzfile["rx_"], npzfile["sx_"], npzfile["akr_"], npzfile["rho_"] + inner_kr_b, base_mfp = npzfile["inner_kr_b_"], npzfile["base_mfp_"] soil = npzfile["soil"] - self.lookup_table = RegularGridInterpolator((rx_, sx_, akr_, rho_), interface) # method = "nearest" fill_value = None , bounds_error=False + self.lookup_table = RegularGridInterpolator((inner_kr_b, base_mfp), interface) # method = "nearest" fill_value = None , bounds_error=False self.sp = vg.Parameters(soil) - + vg.create_mfp_lookup(self.sp) # I do not know why this has to be repeated here + + def open_global_lookup(self, filename): + """opens a global look-up table from a file to quickly find soil root interface potentials, this lookup table works for all van Genuchten parameter sets""" + npzfile = np.load(filename + ".npz") + interface, interface_vg = npzfile["interface"], npzfile["interface_vg"] + vg_m_, inner_kr_b_, base_mfp_, sx_ = npzfile["vg_m_"], npzfile["inner_kr_b_"], npzfile["base_mfp_"], npzfile["sx_"] + self.global_lookup_table = RegularGridInterpolator((vg_m_, inner_kr_b_, base_mfp_), interface) # method = "nearest" fill_value = None , bounds_error=False + self.lookup_global_mfp = RegularGridInterpolator((vg_m_, sx_), interface_vg) + #.sp = vg.Parameters(soil) + + def open_lookup_solutes_simplified(self, filename): + """opens a look-up table from a file to quickly find soil root solute concentrations in the steady state case""" + npzfile = np.load(filename + ".npz") + integral_AdvDiff_ = npzfile["integral_AdvDiff_"] + base_mfp_ = npzfile["base_mfp_"] + soil = npzfile["soil"] + self.lookup_table_sr_solutes_simplified = RegularGridInterpolator((base_mfp_,[0,1]) , integral_AdvDiff_) # method = "nearest" fill_value = None , bounds_error=False + self.sp = vg.Parameters(soil) + vg.create_mfp_lookup(self.sp) # does this have to be repeated here? + + def open_lookup_solutes(self, filename): + """opens an additional look-up table from a file to quickly find soil root solute concentrations in the steady rate case""" + npzfile_sr = np.load(filename + ".npz") + Phi_2_, Phi_1_, Phi_0_, r_eval_ = npzfile_sr["Phi_2_"], npzfile_sr["Phi_1_"], npzfile_sr["Phi_0_"], npzfile_sr["r_eval_"] + soil = npzfile_sr["soil"] + conc_rel_c, Uptake_rel_c, quadratic_rel_c = npzfile_sr["conc_rel_c"], npzfile_sr["Uptake_rel_c"], npzfile_sr["quadratic_rel_c"] + conc_mean_c, Uptake_mean_c, quatratic_mean_c = npzfile_sr["conc_mean_c"], npzfile_sr["Uptake_mean_c"], npzfile_sr["quatratic_mean_c"] + self.lookup_table_sr_solutes = { + "conc_rel_c" : [0], + "conc_mean_c" : [0], + "Uptake_rel_c" : [0], + "Uptake_mean_c" : [0], + "quadratic_rel_c" : [0], + "quatratic_mean_c" : [0] + } + self.lookup_table_sr_solutes["conc_rel_c"] = RegularGridInterpolator((Phi_2_, Phi_1_, Phi_0_, r_eval_) , conc_rel_c) # method = "nearest" fill_value = None , bounds_error=False + self.lookup_table_sr_solutes["conc_mean_c"] = RegularGridInterpolator((Phi_2_, Phi_1_, Phi_0_, r_eval_) , conc_mean_c) + self.lookup_table_sr_solutes["Uptake_rel_c"] = RegularGridInterpolator((Phi_2_, Phi_1_, Phi_0_, r_eval_) , Uptake_rel_c) + self.lookup_table_sr_solutes["Uptake_mean_c"] = RegularGridInterpolator((Phi_2_, Phi_1_, Phi_0_, r_eval_) , Uptake_mean_c) + self.lookup_table_sr_solutes["quadratic_rel_c"] = RegularGridInterpolator((Phi_2_, Phi_1_, Phi_0_, r_eval_) , quadratic_rel_c) + self.lookup_table_sr_solutes["quatratic_mean_c"] = RegularGridInterpolator((Phi_2_, Phi_1_, Phi_0_, r_eval_) , quatratic_mean_c) + self.sp = vg.Parameters(soil) + vg.create_mfp_lookup(self.sp) # does this have to be repeated here? + def soil_root_interface_potentials(self, rx, sx, inner_kr, rho): """ finds matric potentials at the soil root interface for as all segments @@ -63,10 +128,22 @@ def soil_root_interface_potentials(self, rx, sx, inner_kr, rho): rho geometry factor [1] (outer_radius / inner_radius) """ assert len(rx) == len(sx) == len(inner_kr) == len(rho), "rx, sx, inner_kr, and rho must have the same length" + + for i in range(len(rx)): + if rx[i] < -15989: + rx[i] = -15989 + if sx[i] > -1: + sx[i] = -1 + if self.lookup_table: rsx = self.soil_root_interface_potentials_table(rx, sx, inner_kr, rho) else: - rsx = np.array([PerirhizalPython.soil_root_interface_(rx[i], sx[i], inner_kr[i], rho[i], self.sp) for i in range(0, len(rx))]) + if self.global_lookup_table: + rsx = self.soil_root_interface_potentials_table_global(rx, sx, inner_kr, rho) + else: + rsx = np.array([PerirhizalPython.soil_root_interface_(rx[i], sx[i], inner_kr[i], rho[i], self.sp) for i in range(0, len(rx))]) + + return rsx @staticmethod @@ -80,17 +157,575 @@ def soil_root_interface_(rx, sx, inner_kr, rho, sp): rho geometry factor [1] (outer_radius / inner_radius) sp soil parameter: van Genuchten parameter set (type vg.Parameters) """ + #vg.create_mfp_lookup(sp) + #print(vg.fast_mfp[sp](-16000)) if inner_kr < 1.0e-7: return sx - k_soilfun = lambda hsoil, hint: (vg.fast_mfp[sp](hsoil) - vg.fast_mfp[sp](hint)) / (hsoil - hint) # Vanderborgth et al. 2023, Eqn [7] + k_soilfun = lambda hsoil, hint: (vg.fast_mfp[sp](hsoil) - vg.fast_mfp[sp](hint)) * (hsoil - hint) / ((hsoil - hint)**2 + 0.001) # Vanderborgth et al. 2023, Eqn [7] rho2 = np.square(rho) # rho squared b = 2 * (rho2 - 1) / (1 - 0.53 * 0.53 * rho2 + 2 * rho2 * (np.log(rho) + np.log(0.53))) # Vanderborgth et al. 2023, Eqn [8] fun = lambda x: (inner_kr * rx + b * sx * k_soilfun(sx, x)) / (b * k_soilfun(sx, x) + inner_kr) - x + if rx == sx: # degenerate bracket: no gradient, interface equals bulk potential - return sx - # rsx = root_scalar(fun, method="brentq", bracket=[min(rx, sx), max(rx, sx)]) - rsx = fsolve(fun, (rx + sx) / 2) - return rsx + return sx + rsx = root_scalar(fun, method="brentq", bracket=[min(rx, sx), max(rx, sx)]) + return rsx.root + + @staticmethod + def soil_root_interface_simp_(self, inner_kr_b, base_mfp, sp): + """ + finds matric potential at the soil root interface for as single segment + gives (numerically) the same results as "soil_root_interface_", but is much simpler to create a lookup table for + + rx xylem matric potential [cm] + sx bulk soil matric potential [cm] + inner_kr root radius times hydraulic conductivity [cm/day] + rho geometry factor [1] (outer_radius / inner_radius) + sp soil parameter: van Genuchten parameter set (type vg.Parameters) + + inner_kr_b: inner_kr/b + base_mfp: -(inner_kr/b*rx+vg.fast_mfp[sp](sx)) + """ + fun = lambda x: vg.fast_mfp[sp](x) + inner_kr_b * (x - 0*self.h_wilt) +base_mfp + x_int = [-15999.0,-1.0] # interval for search + if fun(x_int[0])> 0: + return x_int[0] + if fun(x_int[1])< 0: + return x_int[1] + rsx = root_scalar(fun, method="brentq", bracket=[x_int[0], x_int[1]]) + return rsx.root + + @staticmethod + def soil_root_interface_global_(self, inner_kr_b, base_mfp, vg_m): + """ + finds matric potential at the soil root interface for as single segment + + rx xylem matric potential [cm] + sx bulk soil matric potential [cm] + inner_kr root radius times hydraulic conductivity [cm/day] + rho geometry factor [1] (outer_radius / inner_radius) + sp soil parameter: van Genuchten parameter set (type vg.Parameters) + + inner_kr_b: (inner_kr * alpha)/(alpha_0 * vg_Ks * b) + base_mfp: -((inner_kr * alpha)/(alpha_0 * Ks * b) * rx + vg.fast_mfp[sp_base*(alpha_0/alpha)](sx)) + vg_m: van genuchten parameter m + + (sp_base is the van genuchten parameter set of Ks = 1 and alpha = alpha_0, m is an input, theta_r and theta_s are not important for the steady rate model computaiton right now) + """ + + + + fun = lambda x: self.lookup_global_mfp((vg_m, x)) + inner_kr_b * x +base_mfp + x_int = [-15990.0,-0.1] # interval for search, get closer to 0 as x = h_sr * alpha / alpha_0, alpha / alpha_0 <= 1 + if fun(x_int[0])> 0: + return x_int[0] + if fun(x_int[1])< 0: + return x_int[1] + rsx = root_scalar(fun, method="brentq", bracket=[x_int[0], x_int[1]]) + return rsx.root + + def solutesuptake_convdiff_(self, watercontent, c_bulk, Vmax, Km, Ds, waterflow, r_root, E, t, sp): + """ + finds solute concentration at the soil root interface for all segments following T. Roose and G. Kirk 2009 doi:10.1007/s11104-008-9777-z + + It assumes a constant watercontent with a constant wateruptake by the root + + watercontent watercontent of the perirhizal zone [cm3/cm3] + c_bulk starting solute concentration at the bulk soil [mol/cm3] + Vmax maximum solute uptake rate, Michaelis Menten Kinetics [mol/(cm2d)] + Km half saturation constant Michaelis Menten Kinetics [mol/cm3] + Ds Diffusion constant in water [cm2/d] + waterflow steady state waterflow (entire root circumference) [cm2/d] + r_root root radius [cm] + E minimum net influx into the plant [mol/cm2d] + t root age [d] + sp van Genuchten parameter set + """ + assert len(watercontent) == len(c_bulk) == len(Vmax) == len(Km) == len(waterflow) == len(r_root) == len(E) == len(t), "error in Perirhizal.py, solutesuptake_convdiff_: watercontent, c_bulk, Vmax, Km, Ds, waterflow, r_root, E and t must have the same length" + + n_segments = len(c_bulk) + segLength = 1 #reference segment length, should not impact the results [cm] + + rsc = np.zeros(n_segments) # solute concentration at the root soil interface + F = np.zeros(n_segments) # uptake of solutes mol/(cm2d) + gamma = 0.577 # Eulers constant + l_func = lambda time : 1/2*np.log(4*np.exp(-gamma)*time+1) + + + for i in range(n_segments): + #compute diffusion coefficient according to Millington and Quirk + D = Ds * math.pow(watercontent[i], 10/3) / (sp.theta_S**2) + + #unit conversions at the end + unitconversion_F = Km[i] * D * r_root[i] # to [mol/cm2d] #TODO: look at this again, it doesn't seem right + unitconversion_c = Km[i] # to [mol/cm3] + + #express the inputs without units + c_inf = c_bulk[i] / Km[i] + Pe = waterflow[i]/D #Peclet number + lamb = segLength*r_root[i]*Vmax[i]/(D*Km[i]) + epsilon = segLength*r_root[i]*E[i]/(D*Km[i]) + + #compute the unitless uptake + if Pe<=(lamb/(1+c_inf)): + F[i] = 2*lamb*(c_inf+epsilon*l_func(t[i])) / (1+c_inf+l_func(t[i])*(lamb + epsilon) +math.sqrt(4*(c_inf+epsilon * l_func(t[i]))+(1-c_inf+(lamb-epsilon)*l_func(t[i]))**2)) #Eqn. [9] + else: + F[i] = 2*lamb*(c_inf+epsilon*l_func(t[i])) / (1+c_inf+l_func(t[i])*(lamb - Pe + epsilon) +math.sqrt(4*(c_inf+epsilon * l_func(t[i]))*(1-Pe*l_func(t[i]))+(1-c_inf+l_func(t[i])*(lamb-Pe-epsilon))**2)) #Eqm. [10] + #compute the unitless solute concentration next to the root + #F(t) = lamb * c / (1+c) - epsilon shortly after Eqn. [10] + #(1+c)*(F(t)+epsilon) = lamb * c + #F(t)+epsilon = (-F(t)-epsilon+lamb) * c + rsc[i] = (F[i]+epsilon)/(lamb-F[i]-epsilon) + + #add units + F[i] = F[i] * unitconversion_F + rsc[i] = rsc[i] * unitconversion_c + + return rsc + + def soil_root_solutes_sr(self, Phi_out, rootwateruptake, waterinflow, r_root, r_prhiz, c_soil, c_outer, Vmax, Km, Ds, sp, mode = "ss"): + """ + steady rate assumption of solute uptake by roots TODO: insert citation + + Phi_out outer matrix flux potential [cm2/d] + rootwateruptake radial root water uptake [cm2/d] + waterinflow radial water inflow at r_prhiz [cm2/d] + r_root root radius [cm] + r_prhiz outer radius [cm] + c_soil solute concentration of the cylinder [mol/cm3] + c_outer solute concentration outside of the cylinder [mol/cm3], only used in the general steady rate case + Vmax maximum radial solute uptake rate, Michaelis Menten Kinetics [mol/(cmd)] + Km half saturation constant Michaelis Menten Kinetics [mol/cm3] + Ds Diffusion constant in water [cm2/d] + sp van Genuchten parameter set + mode what mode is used for the solute uptake? + "ss" for steady state solute uptake + "sr" steady rate solute uptake with a no flux outer BC + "ff" for steady rate uptake in combination with the far field approximation to give an outer Cauchy BC #TODO: cite paper for far field approximation + "dirichlet" for Dirichlet outer BC + + output: + rsc solute concentration next to the root [mol/cm3] + Uptake solute uptake of the root [mol/(cm d)] + + """ + assert len(Phi_out) == len(rootwateruptake) == len(waterinflow) == len(r_root) == len(r_prhiz) == len(c_soil) == len(c_outer) == len(Vmax) == len(Km) == len(Ds), "Phi_soil, rootwateruptake, waterinflow, r_root, r_prhiz, c_soil, c_outer, Vmax, Km and Ds must have the same length" + + n_segments = len(c_soil) + + rsc = np.zeros(n_segments) #concentration next to the root + Uptake = np.zeros(n_segments) #total solute uptake of the root + quadratic_flow = np.zeros(n_segments) #quadratic flow into the parirhizal zone + noflux = np.zeros(n_segments) #concentration where diffusion and advection cancel out + + r_prhiz = r_prhiz.copy() # make sure that the original does not get scaled + r_root = r_root.copy() # make sure that the original does not get scaled + + #solve quadratic eqation for root uptake # TODO: Link to publication + for i in range(n_segments): + #waterflow + rho = r_prhiz[i] / r_root[i] + Phi_0 = Phi_out[i] - 1/2 *(rootwateruptake[i] - waterinflow[i]) / (2*np.pi) * (rho**2) / (1 - rho**2) + Phi_1 = waterinflow[i] / (2*np.pi) - (rootwateruptake[i] - waterinflow[i]) / (2*np.pi) * rho**2 / (1 - rho**2) + Phi_2 = (rootwateruptake[i] - waterinflow[i]) / (2*np.pi) * (rho**2) / (1 - rho**2) + + Ds0 = self.Ds0_ref + scaling = math.sqrt(Ds[i] / Ds0) # Ds_i > Ds0 means scaling >1, the region is modeled to be smaller than it actually is + #scaling = Ds[i] / Ds0 # Ds_i > Ds0 means scaling >1, the region is modeled to be smaller than it actually is + Phi_2 = Phi_2 * (scaling**2) + Phi_1 = Phi_1 + Phi_0 = Phi_0 + Phi_1 * np.log(scaling) + #Phi = lambda r : Phi_2 * r**2/2 + Phi_1 * np.log(r) + Phi_0 + #waterpotential_func = lambda r : vg.fast_imfp[sp](Phi(r)) + #watercontent_func = lambda r : vg.water_content(waterpotential_func(r),sp) + #radial_waterflow = lambda r : 2 * np.pi * (Phi_2 * r**2 + Phi_1) + #Ds_func = lambda r : Ds[i] * math.pow(watercontent_func(r),10/3) / (sp.theta_S**2) # Millington and Quirk + + r_eval = np.logspace(np.log10(r_root[i]), np.log10(r_prhiz[i]), num = 30) + r_eval = [r / scaling for r in r_eval] + r_prhiz[i] = r_prhiz[i] / scaling + r_root[i] = r_root[i] / scaling + + if self.lookup_table_sr_solutes: + for i in range(n): + conc_rel_c[i] = self.lookup_table_sr_solutes["conc_rel_c"]((Phi_2, Phi_1, Phi_0, r_eval[i])) + conc_mean_c[i] = self.lookup_table_sr_solutes["conc_mean_c"]((Phi_2, Phi_1, Phi_0, r_eval[i])) + Uptake_rel_c[i] = self.lookup_table_sr_solutes["Uptake_rel_c"]((Phi_2, Phi_1, Phi_0, r_eval[i])) + Uptake_mean_c[i] = self.lookup_table_sr_solutes["Uptake_mean_c"]((Phi_2, Phi_1, Phi_0, r_eval[i])) + quadratic_rel_c[i] = self.lookup_table_sr_solutes["quadratic_rel_c"]((Phi_2, Phi_1, Phi_0, r_eval[i])) + quatratic_mean_c[i] = self.lookup_table_sr_solutes["quatratic_mean_c"]((Phi_2, Phi_1, Phi_0, r_eval[i])) + else: + conc_rel_c, conc_mean_c, Uptake_rel_c, Uptake_mean_c, quadratic_rel_c, quadratic_mean_c = self.solute_linearequation_sr_(Phi_2, Phi_1, Phi_0, r_eval, sp) + + + #default prefactors for the steady state case + pre_c = 0 + pre_Uptake = 0 + pre_quadratic = 1 + absolute = 0 + + #pre_c * c(r_prhiz) + pre_Uptake * Uptake + pre_quadratic * quadratic = absolute + match mode: + case "ss": #steady state + pre_c, pre_Uptake, pre_quadratic, absolute = 0, 0, 1, 0 + case "sr": #steady rate uptake, no influx + pre_c, pre_Uptake, pre_quadratic, absolute = 0, r_prhiz[i]**2, -1, 0 + case "ff": #far field approximation # TODO insert link to Tiina Roose publication, and current application + B_abs = c_outer[i] * (r_prhiz[i]**2)/ (math.exp(-2*r_prhiz[i]**2)-math.exp(-r_prhiz[i]**2)) + B_cprhiz = -1 * (r_prhiz[i]**2) / (math.exp(-2*r_prhiz[i]**2)-math.exp(-r_prhiz[i]**2)) + + pre_Uptake = 1 + pre_quadratic = r_prhiz[i]**2 + pre_c = Ds_func(r_prhiz[i]) * B_cprhiz * math.exp(-r_prhiz[i]**2) / r_prhiz[i]**2 + waterinflow[i] / (2 * np.pi * r_prhiz[i]) + absolute = Ds_func(r_prhiz[i]) * B_abs * math.exp(-r_prhiz[i]**2) / r_prhiz[i]**2 + case "dirichlet": + pre_c = 1 + pre_Uptake = 0 + pre_quadratic = 0 + absolute = c_outer[i] + case _: #default + pre_c, pre_Uptake, pre_quadratic, absolute = 0, 0, 1, 0 + print("The mode of a steady rate solute flow approximation could not be identified, use steady state.") + + #variables to eliminate: Uptake, quadratic, c_prhiz + # variable for the final equations: rsc + A = np.zeros((4,4)) + b = np.zeros((4,2)) #right hand side of the linear equation, one for absolute values, one for the rsc value + + #linear equations for c_prhiz, Uptake, quadratic, rsc (3 linear, 1 quadratic) + + #c11 * Uptake + c12 * quadratic + c13 * c0 = c_mean + A[0,0] = (Uptake_mean_c[-1]*(r_prhiz[i]**2)- Uptake_mean_c[0]*(r_root[i]**2))/(r_prhiz[i]**2-r_root[i]**2) + A[0,1] = (quadratic_mean_c[-1]*(r_prhiz[i]**2)- quadratic_mean_c[0]*(r_root[i]**2))/(r_prhiz[i]**2-r_root[i]**2) + A[0,2] = (conc_mean_c[-1]*(r_prhiz[i]**2)- conc_mean_c[0]*(r_root[i]**2))/(r_prhiz[i]**2-r_root[i]**2) + b[0,0] = c_soil[i] + + #c21 * Uptake + c22 * quadratic + c23 * c0 = c_prhiz + A[1,0] = Uptake_rel_c[-1] + A[1,1] = quadratic_rel_c[-1] + A[1,2] = conc_rel_c[-1] + A[1,3] = -1 + + #pre_uptake * Uptake + pre_quadratic * quadratic + pre_c * c_prhiz = absolute + A[2,0] = pre_Uptake + A[2,1] = pre_quadratic + A[2,3] = pre_c + b[2,0] = absolute + + # c41 * Uptake + c42 * quadratic + c43 * c0 = rsc + A[3,0] = Uptake_rel_c[0] + A[3,1] = quadratic_rel_c[0] + A[3,2] = conc_rel_c[0] + b[3,1] = 1 + + #m*rsc+n = [ss uptake, sr uptake, c(rprhiz), Uptake] + n = la.solve(A,b[:,0]) + m = la.solve(A,b[:,1]) + + #solve quadratic equation + #Uptake means radial uptake + #Uptake = Vmax * rsc / (Km + rsc) + #Uptake = (m[0]+r_root**2 * m[1])*rsc+(n[0]+r_root**2 * n[1]) =: w * rsc + v + #(Km + rsc) * (w*rsc+v) = Vmax * rsc + # w * rsc**2 + (Km[i]*w+v-Vmax[i]) * rsc + Km[i]*v = 0 + w = (m[0] + (r_root[i]**2) * m[1]) #*(2*np.pi*r_root[i]) + v = (n[0] + (r_root[i]**2) * n[1]) #*(2*np.pi*r_root[i]) + tol = 1e-9 #arbitrary for now, TODO think of a reasonable tolerance + + if abs(w) Ds0_ref means scaling >1, the region is modeled to be smaller than it actually is + #scaling = Ds0 / self.Ds0_ref # Ds0 > Ds0_ref means scaling >1, the region is modeled to be smaller than it actually is + r_eval = [r / scaling for r in r_eval] + Phi_2 = Phi_2 * (scaling**2) + Phi_1 = Phi_1 + Phi_0 = Phi_0 + Phi_1 * np.log(scaling) + + #values necessary for plotting + n = len(r_eval) + conc_rel_c, conc_mean_c, Uptake_rel_c, Uptake_mean_c, quadratic_rel_c, quadratic_mean_c = np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n), np.zeros(n) + + if self.lookup_table_sr_solutes: + for i in range(n): + conc_rel_c[i] = self.lookup_table_sr_solutes["conc_rel_c"]((Phi_2, Phi_1, Phi_0, r_eval[i])) + conc_mean_c[i] = self.lookup_table_sr_solutes["conc_mean_c"]((Phi_2, Phi_1, Phi_0, r_eval[i])) + Uptake_rel_c[i] = self.lookup_table_sr_solutes["Uptake_rel_c"]((Phi_2, Phi_1, Phi_0, r_eval[i])) + Uptake_mean_c[i] = self.lookup_table_sr_solutes["Uptake_mean_c"]((Phi_2, Phi_1, Phi_0, r_eval[i])) + quadratic_rel_c[i] = self.lookup_table_sr_solutes["quadratic_rel_c"]((Phi_2, Phi_1, Phi_0, r_eval[i])) + quatratic_mean_c[i] = self.lookup_table_sr_solutes["quatratic_mean_c"]((Phi_2, Phi_1, Phi_0, r_eval[i])) + else: + conc_rel_c, conc_mean_c, Uptake_rel_c, Uptake_mean_c, quadratic_rel_c, quadratic_mean_c = self.solute_linearequation_sr_(Phi_2, Phi_1, Phi_0, r_eval, sp) + + soluteconcentration = noflux * conc_rel_c + Uptake * Uptake_rel_c + quadratic_flow * quadratic_rel_c + soluteconcentration_mean = noflux * conc_mean_c + Uptake * Uptake_mean_c + quadratic_flow * quadratic_mean_c + + return watercontent, waterpotential, soluteconcentration, soluteconcentration_mean + + def solute_linearequation_sr_(self, Phi_2, Phi_1, Phi_0, r_eval, sp): + """ + computes a linear equation from the diffusion advection ODE on the perirhizal solute flow #TODO: insert equation number + the radius (and steady rate wateruptake) are scaled by the diffusion coefficient of the solute already included in the inputs + + Phi1 matrix flux potential at \tilde{r} = 1 [cm2/d] + waterflow waterflow at r = 1 [cm/d] + wateruptake steady rate wateruptake throughout the perirhizal zone, scaled to \tilde{r} [cm3/(cm3d)] + r_eval the equation is solved from r_eval[0] to r_eval[-1], it is a 1d array [cm] + + outputs: + conc_rel_c ratio of c(r_rel[i]) to c(r_rel[0]) for d_r c(1)=0 + conc_mean_c mean solute concentrations for the steady state solute uptake + inflow_rel_c ratio of c(r_rel[i]) to d_r c(1) for c(1)=0 + inflow_mean_c mean solute concentration for c(1)=0 + Uptake_rel_c ratio of c(r_rel[i]) to rel_soluteUptake for c(1)=0, d_r c(1)=0 (steady rate uptake case) + Uptake_mean_c mean solute concentrations relative to rel_soluteUptake + """ + + #reference diffusion coefficient to which r was scaled to [cm2/d] + Ds0 = self.Ds0_ref + + #begin at a predefined + n = len(r_eval) + r0 = self.r0_ref + insertion = False #was the radius inserted in the beginning? + if r_eval[0]==r0: + pass + else: + r_eval = np.insert(r_eval, 0, r0, axis=0) + insertion = True + n = n+1 + + conc_rel_c = np.ones(n) + conc_mean_c = np.ones(n) + inflow_rel_c = np.ones(n) + inflow_mean_c = np.ones(n) + Uptake_rel_c = np.ones(n) + Uptake_mean_c = np.ones(n) + + #used for weighting + watercontent_times_r = np.ones(n) + + #waterflow + Phi = lambda r : Phi_2 * r**2/2 + Phi_1 * np.log(r) + Phi_0 + radial_waterflow = lambda r : - 2 * np.pi * (Phi_2 * r**2 + Phi_1) + watercontent = lambda r : vg.water_content(vg.fast_imfp[sp](Phi(r)), sp) + Ds = lambda r : Ds0 * math.pow(watercontent(r),10/3) / (sp.theta_S**2) + + f_homogen = lambda c, r : ( radial_waterflow(r) * c) / (2 * np.pi * Ds(r) * r * watercontent(r)) #TODO: add reference + f_steadystate = lambda c, r : (1 + radial_waterflow(r) * c) / (2 * np.pi * Ds(r) * r * watercontent(r)) #TODO: add reference + f_quadratic = lambda c, r : (r**2 + radial_waterflow(r) * c) / (2 * np.pi * Ds(r) * r * watercontent(r)) #TODO: add reference + + #numerically solve the ODE c' = f(r_rel,c) starting at r_rel=1 + conc_rel_c = odeint(f_homogen, y0 = 1, t = r_eval) + conc_rel_c = np.array([element[0] for element in conc_rel_c]) + inflow_rel_c = odeint(f_steadystate, y0 = 0, t = r_eval) + inflow_rel_c = np.array([element[0] for element in inflow_rel_c]) + Uptake_rel_c = odeint(f_quadratic, y0 = 0, t = r_eval) + Uptake_rel_c = np.array([element[0] for element in Uptake_rel_c]) + + #compute the means + conc_mean_c[0] = conc_rel_c[0] + inflow_mean_c[0] = inflow_rel_c[0] + Uptake_mean_c[0] = Uptake_rel_c[0] + watercontent_times_r[0] = 0 + watercontent_times_r[1:] = np.array([watercontent(r_eval[i])*(r_eval[i]**2-r_eval[i-1]**2) for i in range(1,n)]) + conc_mean_c[1:] = np.array([np.average(conc_rel_c[0:(j+1)], weights=watercontent_times_r[0:(j+1)]) for j in range(1,n)]) + inflow_mean_c[1:] = np.array([np.average(inflow_rel_c[0:(j+1)], weights=watercontent_times_r[0:(j+1)]) for j in range(1,n)]) + Uptake_mean_c[1:] = np.array([np.average(Uptake_rel_c[0:(j+1)], weights=watercontent_times_r[0:(j+1)]) for j in range(1,n)]) + + #remove the inserted start at self.r0_ref + if insertion: + return conc_rel_c[1:], conc_mean_c[1:], inflow_rel_c[1:], inflow_mean_c[1:], Uptake_rel_c[1:], Uptake_mean_c[1:] + else: + return conc_rel_c, conc_mean_c, inflow_rel_c, inflow_mean_c, Uptake_rel_c, Uptake_mean_c + + + def soil_root_solutes_steadyrate_simplified_(self, Phi_root, Phi_soil, r_root, r_prhiz, c_bulk, Vmax, Km, Ds, waterflow, sp, n_approx = 1): + """ + finds solute concentration at the soil root interface for all segments assuming that the waterflow and soluteflow are proportional throughout the perirhizal zone + in this simplification the mode of soluteuptake (steady state, steady rate with no influx, other steady rate) is prescribed by the water uptake + uses a look up tables if present (see create_lookup, and open_lookup) #TODO: add links to the exact lookup table + + Phi_root matrix flux potential at the root-soil-interface [cm2/d] + Phi_soil matrix flux potential at the bulk soil [cm2/d] + r_root root radius [cm] + r_prhiz perirhizal radius [cm] + c_bulk solute concentration at the bulk soil [mol/cm3] + Vmax maximum solute uptake rate, Michaelis Menten Kinetics [mol/(cm2d)] + Km half saturation constant Michaelis Menten Kinetics [mol/cm3] + Ds Diffusion constant in water [cm2/d] + waterflow steady state waterflow [cm/d] + sp van Genuchten parameter set + n_approx how many points are used to approximate the perirhizal solute concentration? + n_approx = 1 means only approximation at the perirhizal radius + """ + assert len(Phi_root) == len(Phi_soil) == len(r_root) == len(r_prhiz) == len(c_bulk) == len(Vmax) == len(Km) == len(waterflow), "Phi_root, Phi_soil, r_root, r_prhiz, c_bulk, Vmax, Km and waterflow must have the same length" + + n_segments = len(c_bulk) + + #the radii at which the solute concentration will be tested. They are given as relative values between r_root and r_prhiz. r=1 means r_prhiz. + # approx(int[fun,0,1]) = sum(weights[i] * fun(x(i))) + if n_approx ==1: + x = [1] + weights = [1] + else: + [x,weights] = ts_roots(n_approx) # returns the nodes and weights for the chebyshev numerical integration + + + rho = [r_prhiz[i]/r_root[i] for i in range(n_segments)] + c_sol_mean2root = np.zeros(n_segments) # ratio of the mean solute concentration to solute concentration next to the root + mean_watercontent = np.zeros(n_segments) #watercontent next to the root + rsc = np.zeros(n_segments) #solute concentration next to the root (output) + F = np.zeros(n_segments) #F is a helper value for computing the ratio between mean solute concentraiton and solute concentration next to the root + F0 = np.zeros(n_segments) #helper value for the ratio of solute next to the root + + #compute a prefactor + D_tilde = 1/Ds/math.pow(sp.theta_S-sp.theta_R,13/3)*(sp.theta_S*sp.theta_S) + + #determine mfp depending on r + Phi_A = np.zeros(n_segments) + Phi_C = np.zeros(n_segments) + for i in range(n_segments): + Phi_A[i], Phi_C[i] = self.determine_mfp_function(Phi_root[i], Phi_soil[i], rho[i]) #Phi(r/r_prhiz)= A(s^2-ln(s^2))+C + + + + #determine F0 + if self.lookup_table_sr_solutes_simplified: + F0=[self.lookup_table_sr_solutes_simplified((Phi_root[i],0)) for i in range(n_segments)] + else: + F0=[self.integral_AdvectionDiffusion_(Phi_root[i],self.sp) for i in range(n_segments)] + + #compute the ratio of (c_mean-C_SW) to (c_root-C_SW) (weighted by water content and cylinder volume) + for j in range(n_approx): + test_radii = np.array([math.sqrt(x[j] * (r_prhiz[i]**2-r_root[i]**2)+r_root[i]**2) for i in range(n_segments)]) #scale the chebyshev nodes according to volume + s = [test_radii[i]/r_prhiz[i] for i in range(n_segments)] + Phi_current = [Phi_A[i]* (s[i]**2-2*np.log(s[i])) + Phi_C[i] for i in range(n_segments)] + current_watercontent = [vg.water_content(vg.fast_imfp[self.sp](Phi_current[i]),self.sp) for i in range(n_segments)] + if self.lookup_table_sr_solutes_simplified: + F=[(self.lookup_table_sr_solutes_simplified((Phi_current[i],0))-F0[i]) for i in range(n_segments)] + else: + F=[(self.integral_AdvectionDiffusion_(Phi_current[i],self.sp)-F0[i]) for i in range(n_segments)] + for i in range(n_segments): + c_sol_mean2root[i] += weights[j] * math.exp(max(min(D_tilde*F[i],10),-10)) * current_watercontent[i] + mean_watercontent[i] += weights[j] * current_watercontent[i] + for i in range(n_segments): + c_sol_mean2root[i] = c_sol_mean2root[i] / mean_watercontent[i] + + #solve quadratic eqation # TODO: Link to publication + for i in range(n_segments): + a1=c_bulk[i]/c_sol_mean2root[i] + a2=(1-1/c_sol_mean2root[i])/(waterflow[i]) + p=Km[i]-a2*Vmax[i]-a1 + q=-Km[i]*a1 + r1=-p/2-math.sqrt(pow(p/2,2)-q) + r2=-p/2+math.sqrt(pow(p/2,2)-q) + if r1<0: + rsc[i]=r2 + else: + rsc[i]=r1 + return rsc + + @staticmethod + def integral_AdvectionDiffusion_(Phi_input, sp): + """ + Integral of (D_s(theta_s-theta_r)^(13/3))/(theta*D(theta)) from Phi = 0.001 to Phi_input + + Phi matrix flux potential [cm2/d] + sp soil parameter: van Genuchten parameter set (type vg.Parameters) + """ + if Phi_input <=0: + return 0 + theta_rel = sp.theta_R/(sp.theta_S-sp.theta_R) + integral_fun = lambda Phi: pow(theta_rel+vg.effective_saturation(vg.fast_imfp[sp](Phi),sp),-13/3) + integral_AdvDiff, _ = integrate.quad(integral_fun, 1.0e-3, Phi_input) + + return integral_AdvDiff + + + def determine_mfp_function(self, Phi_root, Phi_soil, rho): + """ + determine the spatial function of the mfp + Phi_root matrix flux potential at the root-soil-interface [cm2/d] + Phi_soil matrix flux potential at the bulk soil [cm2/d] + rho geometry factor (outer_radius / inner_radius) [1] + + Phi(r/r_prhiz)= A(s^2-ln(s^2))+C + (steady rate kinetics with a no flux outer BC at s = 1 + + A (0.53^2-ln(0.53^2))+C=Phi_soil + A (0.01^2-ln(0.01^2))+C=Phi_root + + A = (Phi_soil-Phi_root)/(0.53^2-ln(0.53^2)-0.01^2+ln(0.01^2) + C = ((0.53^2-ln(0.53^2)) * Phi_root - (0.01^2-ln(0.01^2)) * Phi_soil)/(0.53^2-ln(0.53^2)-0.01^2+ln(0.01^2) + """ + + det = 0.53**2-2*math.log(0.53)-(1/rho)**2+2*math.log(1/rho) + a = 0.53**2-2*math.log(0.53) + c = (1/rho)**2-2*math.log(1/rho) + + Phi_A = min((Phi_soil-Phi_root) / det, 0) + Phi_C = max((a*Phi_root-c*Phi_soil) / det + c * Phi_A, 1.0e-6) - c * Phi_A + + return Phi_A, Phi_C + def perirhizal_conductance_per_layer(self, h_bs, h_sr, sp): """ @@ -113,8 +748,11 @@ def perirhizal_conductance_per_layer(self, h_bs, h_sr, sp): l_root = rld * dz # [cm-1] return 2 * np.pi * np.multiply(l_root, np.multiply(b, k_prhiz(h_bs, h_sr))) # [day-1], see Vanderborght et al. (2023), Eqn [6] + # depricated, leave the function in, in case we ever want to compute big lookup tables again def create_lookup_mpi(self, filename, sp): """ + depricated + Precomputes all soil root interface potentials for a specific soil type and saves results into a 4D lookup table @@ -122,6 +760,10 @@ def create_lookup_mpi(self, filename, sp): sp van genuchten soil parameters, , call vg.create_mfp_lookup(sp) before """ + + if rank == 0: + print("The function *create_lookup_mpi* is deprecated, use *create_lookup* without mpi instead.") + from mpi4py import MPI comm = MPI.COMM_WORLD @@ -139,7 +781,7 @@ def create_lookup_mpi(self, filename, sp): akrn = 100 akr_ = np.logspace(np.log10(1.0e-7), np.log10(1.0e-4), akrn) rhon = 30 - rho_ = np.logspace(np.log10(1.0), np.log10(200.0), rhon) + rho_ = np.logspace(np.log10(2.0), np.log10(200.0), rhon) if rank == 0: print(filename, "calculating", rxn * sxn * rhon * akrn, "supporting points on", size, "thread(s)") @@ -200,35 +842,213 @@ def create_lookup_mpi(self, filename, sp): def create_lookup(self, filename, sp): """ Precomputes all soil root interface potentials for a specific soil type - and saves results into a 4D lookup table + and saves results into a 2D lookup table filename three files are written (filename, filename_, and filename_soil) sp van genuchten soil parameters, , call vg.create_mfp_lookup(sp) before """ - rxn = 150 - rx_ = -np.logspace(np.log10(1.0), np.log10(16000), rxn) - rx_ = rx_ + np.ones((rxn,)) - rx_ = rx_[::-1] - sxn = 150 - sx_ = -np.logspace(np.log10(1.0), np.log10(16000), sxn) - sx_ = sx_ + np.ones((sxn,)) - sx_ = sx_[::-1] - akrn = 100 - akr_ = np.logspace(np.log10(1.0e-7), np.log10(1.0e-4), akrn) - rhon = 30 - rho_ = np.logspace(np.log10(1.0), np.log10(200.0), rhon) - interface = np.zeros((rxn, sxn, akrn, rhon)) - for i, rx in enumerate(rx_): - print(i) - for j, sx in enumerate(sx_): - for k, akr in enumerate(akr_): - for l, rho in enumerate(rho_): - interface[i, j, k, l] = PerirhizalPython.soil_root_interface_(rx, sx, akr, rho, sp) - np.savez(filename, interface=interface, rx_=rx_, sx_=sx_, akr_=akr_, rho_=rho_, soil=list(sp)) - self.lookup_table = RegularGridInterpolator((rx_, sx_, akr_, rho_), interface) + + + # intervals for all inputs + rx_int_abs = [1.0, 15999.0] #absolute values for logarithmic scaling + sx_int_abs = [1.0, 15999.0] + akr_int = [1.0e-7,1.0e-4] + rho_int = [2.0,200.0] + + #the lookup table only needs 2 inputs from combinations of the 4 general inputs + b_func = lambda rho: 2 * (rho*rho - 1) / (1 - 0.53 * 0.53 * rho*rho + 2 * rho*rho * (np.log(rho) + np.log(0.53))) # Vanderborgth et al. 2023, Eqn [8] + b_0 = b_func(rho_int[0]) + b_1 = b_func(rho_int[1]) + b_int = [min(b_0,b_1),max(b_0,b_1)] + + #compute the intervals of these 2 possible inputs for the lookup table + + inner_kr_bn = 100 + inner_kr_b_int = [akr_int[0]/b_int[1],akr_int[1]/b_int[0]] + inner_kr_b_= np.logspace(np.log10(inner_kr_b_int[0]), np.log10(inner_kr_b_int[1]), inner_kr_bn) + + min_int, max_int = vg.fast_mfp[sp](-sx_int_abs[1]), vg.fast_mfp[sp](-sx_int_abs[0]) + base_mfp_n_1 = 500 + base_mfp_n_2 = 500 + base_mfp_n = base_mfp_n_1 + base_mfp_n_2 + base_mfp_int = [inner_kr_b_int[0]*(0*self.h_wilt-(-rx_int_abs[0])) - vg.fast_mfp[sp](-sx_int_abs[0]), inner_kr_b_int[1]*(0*self.h_wilt-(-rx_int_abs[1])) - vg.fast_mfp[sp](-sx_int_abs[1])] + #base_mfp_int are negative and positive + base_mfp_interval1= -np.logspace(np.log10(-base_mfp_int[0]), np.log10(1.0e-9), base_mfp_n_1) + base_mfp_interval2= np.logspace(np.log10(1.0e-9), np.log10(base_mfp_int[1]), base_mfp_n_2) + base_mfp_ = np.concatenate((base_mfp_interval1,base_mfp_interval2)) + print("boundaries lookup", base_mfp_[0], base_mfp_[1]) + interface = np.zeros((inner_kr_bn,base_mfp_n)) + + print("Creating a lookup table for the hydraulic perirhizal resistance model") + + for i, inner_kr_b in enumerate(inner_kr_b_): + print(i, " / ", inner_kr_bn) + for j, base_mfp in enumerate(base_mfp_): + interface[i, j] = PerirhizalPython.soil_root_interface_simp_(self, inner_kr_b, base_mfp, sp) + np.savez(filename, interface=interface, inner_kr_b_=inner_kr_b_, base_mfp_=base_mfp_, soil=list(sp)) + self.lookup_table = RegularGridInterpolator((inner_kr_b_, base_mfp_), interface) self.sp = sp + + print("Done with creating a lookup table for the hydraulic perirhizal resistance model") + return 0 + def create_lookup_global(self, filename, dummy_sp = vg.Parameters([0.078, 0.43, 0.036, 1.56, 24.96])): + """ + Precomputes all soil root interface potentials for a specific soil type + and saves results into a 3D lookup table + + filename three files are written (filename, filename_, and filename_soil) + dummy_sp van genuchten soil parameters set for the perirhizal zone, not needed in the actual lookup table + has been set to hydrus loam as a default + """ + + + # intervals for all inputs + rx_int_abs = [0.1, 16000.0] #absolute values for logarithmic scaling + sx_int_abs = [0.09, 16000.0] + akr_int = [1.0e-7,1.0e-4] + rho_int = [2.0,200.0] + vg_m_int = [0.1,1.5] + vg_Ks_int = [1.0,100.0] + vg_alpha_int = [0.01,0.3] + + #the lookup table only needs 2 inputs from combinations of the 4 general inputs + b_func = lambda rho: 2 * (rho*rho - 1) / (1 - 0.53 * 0.53 * rho*rho + 2 * rho*rho * (np.log(rho) + np.log(0.53))) # Vanderborgth et al. 2023, Eqn [8] + b_0 = b_func(rho_int[0]) + b_1 = b_func(rho_int[1]) + b_int = [min(b_0,b_1),max(b_0,b_1)] + + vg_m_n = 100 + vg_m_ = np.logspace(np.log10(vg_m_int[0]), np.log10(vg_m_int[1]), vg_m_n) + + sx_n = 200 + sx_ = -np.logspace(np.log10(sx_int_abs[0]), np.log10(sx_int_abs[1]), sx_n) + + interface_vg = np.zeros((vg_m_n,sx_n)) + + #construct a smaller lookup table for the simplified matrix flux potential + print("Creating a general lookup table for the matrix flux potential") + for i, vg_m in enumerate(vg_m_): + print(i, " / ", vg_m_n) + sp_dummy = vg.Parameters([0.1,0.4,self.alpha_0,1./(1.-vg_m),1.0]) + vg.create_mfp_lookup(sp_dummy, verbose=False) + for j, sx in enumerate(sx_): + interface_vg[i, j] = vg.matric_flux_potential(sx, sp_dummy) + self.lookup_global_mfp = RegularGridInterpolator((vg_m_, sx_), interface_vg) + + + inner_kr_bn = 200 + inner_kr_b_abs_int = [akr_int[0]/(b_int[1]*vg_Ks_int[1]),akr_int[1]/(b_int[0]*vg_Ks_int[0])] + inner_kr_b_= np.logspace(np.log10(inner_kr_b_abs_int[0]), np.log10(inner_kr_b_abs_int[1]), inner_kr_bn) + + #base_mfp_n = 200 + base_mfp_n_1 = 100 + base_mfp_n_2 = 100 + base_mfp_n = base_mfp_n_1 + base_mfp_n_2 + #base_mfp_int are negative and positive + base_mfp_int = [inner_kr_b_abs_int[0]*(-(-rx_int_abs[0]) * vg_alpha_int[0]/self.alpha_0) - self.lookup_global_mfp((vg_m_int[1],min(-sx_int_abs[0] * vg_alpha_int[0]/self.alpha_0,-1))), inner_kr_b_abs_int[1]*(-(-rx_int_abs[1]) *vg_alpha_int[1]/self.alpha_0) -self.lookup_global_mfp((vg_m_int[0],-sx_int_abs[1] * vg_alpha_int[1]/self.alpha_0))] + base_mfp_interval1= -np.logspace(np.log10(-base_mfp_int[0]), np.log10(1.0e-9), base_mfp_n_1) + base_mfp_interval2= np.logspace(np.log10(1.0e-9), np.log10(base_mfp_int[1]), base_mfp_n_2) + base_mfp_ = np.concatenate((base_mfp_interval1,base_mfp_interval2)) + + interface = np.zeros((vg_m_n,inner_kr_bn,base_mfp_n)) + + tol = 0.001 #tolerance for extreme values of the matrix flux potential at the soil -> output h_rs = h_soil or h_x + + print("Creating a global lookup table for the hydraulic perirhizal resistance model") + for i, vg_m in enumerate(vg_m_): + print(i, " / ", vg_m_n) + for j, inner_kr_b in enumerate(inner_kr_b_): + for k, base_mfp in enumerate(base_mfp_): + interface[i, j, k] = PerirhizalPython.soil_root_interface_global_(self, inner_kr_b, base_mfp, vg_m) + np.savez(filename, interface=interface, interface_vg=interface_vg, vg_m_=vg_m_, inner_kr_b_=inner_kr_b_, base_mfp_=base_mfp_, sx_=sx_, soil=list(sp)) + self.global_lookup_table = RegularGridInterpolator((vg_m_, inner_kr_b_, base_mfp_), interface) + self.sp = dummy_sp + + print("Done with creating a general lookup table for the matrix flux potential") + return 0 + + def create_integralAdvectionDiffusion_lookup(self, filename, sp): + """ + Precomputes all integrals for the steady state solute flow + + Phi upper matrix flux potential (bottom is set to 0) [cm2/d] + sp van genuchten soil parameters, , call + vg.create_mfp_lookup(sp) before + """ + Phin = 300 + Phi_ = np.logspace(np.log10(1.0e-9), np.log10(500), Phin) + base_mfp_=Phi_ + integral_AdvDiff_ = np.zeros((Phin, 2)) + print("Creating a lookup table for the steady state solute flow") + for i, Phi in enumerate(Phi_): + print(i, " / ", Phin) + integral_AdvDiff_[i,0] = PerirhizalPython.integral_AdvectionDiffusion_(Phi, sp) + integral_AdvDiff_[i,1] = integral_AdvDiff_[i,0] + np.savez(filename, integral_AdvDiff_ = integral_AdvDiff_, base_mfp_ = base_mfp_, soil = list(sp)) + self.lookup_table_sr_solutes_simplified = RegularGridInterpolator((base_mfp_,[0,1]) , integral_AdvDiff_) + self.sp = sp + return integral_AdvDiff_, base_mfp_ + + def create_solute_sr_lookup(self, filename, sp): + """ + Precomputes for the steady rate solute flow + + Phi_2, Phi_1, Phi_0 Phi = lambda r : Phi_2 * r**2/2 + Phi_1 * np.log(r) + Phi_0 + r_eval radii at which the conenctrations are avaluated + sp van genuchten soil parameters, , call + vg.create_mfp_lookup(sp) before + """ + + #repeat steady state case + + n_Ds = len(Ds_) + + Phi_2_n = 100 + Phi_1_n = 100 + Phi_0_n = 100 + r_eval_n = 100 + r_eval_ = np.logspace(np.log(self.r0_ref), np.log(10)) + + conc_rel_c = np.zeros((Phi_2_n, Phi_1_n, Phi_0_n, r_eval_n)) + conc_mean_c = np.zeros((Phi_2_n, Phi_1_n, Phi_0_n, r_eval_n)) + Uptake_rel_c = np.zeros((Phi_2_n, Phi_1_n, Phi_0_n, r_eval_n)) + Uptake_mean_c = np.zeros((Phi_2_n, Phi_1_n, Phi_0_n, r_eval_n)) + quadratic_rel_c = np.zeros((Phi_2_n, Phi_1_n, Phi_0_n, r_eval_n)) + quatratic_mean_c = np.zeros((Phi_2_n, Phi_1_n, Phi_0_n, r_eval_n)) + + + Phi_2_ = -np.logspace(np.log(1.0e-5), np.log(1.0e-9), Phi_2_n) + print("Creating a big lookup table for the steady rate solute flow") + for i, Phi_2 in enumerate(Phi_2_): + print("Done with ", str(i+1), " out of ", str(Phi_2_n)) + Phi_1_ = np.logspace(np.log(1.0e-5), np.log(1.0e-1), Phi_1_n) + for j, Phi_1 in enumerate(Phi_1_): + Phi_0_ = np.logspace(np.log(1.0e-9), np.log(1.0e1), Phi_0_n) + for k, Phi_0 in enumerate(Phi_0_): + conc_rel_c[i,j,k,:], conc_mean_c[i,j,k,:], Uptake_rel_c[i,j,k,:], Uptake_mean_c[i,j,k,:], quadratic_rel_c[i,j,k,:], quatratic_mean_c[i,j,k,:] = self.solute_linearequation_sr_(Phi_2, Phi_1, Phi_0, r_eval_, sp) + + np.savez(filename, Phi_2_n=Phi_2_n, Phi_1_n=Phi_1_n, Phi_0_n = Phi_0_n, r_eval_n = r_eval_n, conc_rel_c = conc_rel_c, conc_mean_c = conc_mean_c, Uptake_rel_c = Uptake_rel_c, Uptake_mean_c = Uptake_mean_c, quatratic_rel_c = quatratic_rel_c, quatratic_mean_c = quatratic_mean_c, soil = list(sp)) + #like 3 GB? + self.lookup_table_sr_solutes = { + "conc_rel_c" : conc_rel_c, + "conc_mean_c" : conc_mean_c, + "Uptake_rel_c" : Uptake_rel_c, + "Uptake_mean_c" : Uptake_mean_c, + "quadratic_rel_c" : quadratic_rel_c, + "quatratic_mean_c" : quatratic_mean_c + } + self.lookup_table_sr_solutes["conc_rel_c"] = RegularGridInterpolator((Phi_2_, Phi_1_, Phi_0_, r_eval_) , conc_rel_c) # method = "nearest" fill_value = None , bounds_error=False + self.lookup_table_sr_solutes["conc_mean_c"] = RegularGridInterpolator((Phi_2_, Phi_1_, Phi_0_, r_eval_) , conc_mean_c) + self.lookup_table_sr_solutes["Uptake_rel_c"] = RegularGridInterpolator((Phi_2_, Phi_1_, Phi_0_, r_eval_) , Uptake_rel_c) + self.lookup_table_sr_solutes["Uptake_mean_c"] = RegularGridInterpolator((Phi_2_, Phi_1_, Phi_0_, r_eval_) , Uptake_mean_c) + self.lookup_table_sr_solutes["quadratic_rel_c"] = RegularGridInterpolator((Phi_2_, Phi_1_, Phi_0_, r_eval_) , quadratic_rel_c) + self.lookup_table_sr_solutes["quatratic_mean_c"] = RegularGridInterpolator((Phi_2_, Phi_1_, Phi_0_, r_eval_) , quatratic_mean_c) + self.sp = vg.Parameters(soil) + vg.create_mfp_lookup(self.sp) # does this have to be repeated here? + + def soil_root_interface_potentials_table(self, rx, sx, inner_kr_, rho_): """ finds potential at the soil root interface using a lookup table @@ -239,11 +1059,34 @@ def soil_root_interface_potentials_table(self, rx, sx, inner_kr_, rho_): rho geometry factor [1] f function to look up the potentials """ + + """ + the original equations were + vg.fast_mfp[sp](x) = int_(h_wilting)^(x)K(S(h))dh + k_soilfun(sx, h_sr)= (vg.fast_mfp[sp](h_sr)-vg.fast_mfp[sp](sx))/(h_sr-sx) Vanderborgth et al. 2023, Eqn [7] + (inner_kr * rx + b * sx * k_soilfun(sx, x)) / (b * k_soilfun(sx, x) + inner_kr) - x = 0 + + equivalent: + vg.fast_mfp[sp](x) + (inner_kr / b) * x - (inner_kr / b * rx + vg.fast_mfp[sp](sx)) = 0 + vg.fast_mfp[sp](x) + inner_kr_b * x + base_mfp = 0 + x only depends on 2 variables + + """ + try: sx = np.array(sx) mask = inner_kr_ == 0 inner_kr_[mask] = 1.0e-7 - rsx = self.lookup_table((rx, sx, inner_kr_, rho_)) + + #compute two inputs for the lookup table + rho = np.array(rho_) + rho[rho > 199] = 199#maximal for rho + rho2 = np.multiply(rho,rho) + b = 2 * (rho2 - 1) / (1 - 0.53 * 0.53 * rho2 + 2 * rho2 * (np.log(rho) + np.log(0.53))) # Vanderborgth et al. 2023, Eqn [8] + inner_kr_b = np.divide(inner_kr_,b) + base_mfp = - (inner_kr_b * rx + vg.fast_mfp[self.sp](sx)) + + rsx = self.lookup_table((inner_kr_b, np.array([max(min(base_mfp[i],-51.01),-53.59) for i in range(len(base_mfp))]))) rsx[mask] = sx[mask] # if inner_kr is zero, there is no flow, and the interface potential is the same as the soil potential except: if np.max(rx) > 0: @@ -281,6 +1124,90 @@ def soil_root_interface_potentials_table(self, rx, sx, inner_kr_, rho_): raise return rsx + + def soil_root_interface_potentials_table_global(self, rx, sx, inner_kr_, rho_): + """ + finds potential at the soil root interface using a globla lookup table, the lookup table works for all van Genuchten parameter sets + + rx xylem matric potential [cm] + sx bulk soil matric potential [cm] + inner_kr root radius times hydraulic conductivity [cm/day] + rho geometry factor [1] + f function to look up the potentials + """ + + """ + the original equations were + vg.fast_mfp[sp](x) = int_(h_wilting)^(x)K(S(h))dh + k_soilfun(sx, h_sr)= (vg.fast_mfp[sp](h_sr)-vg.fast_mfp[sp](sx))/(h_sr-sx) Vanderborgth et al. 2023, Eqn [7] + (inner_kr * rx + b * sx * k_soilfun(sx, x)) / (b * k_soilfun(sx, x) + inner_kr) - x = 0 + + equivalent: + vg.fast_mfp[sp](x) = alpha_0/alpha*vg_Ks * int_(h_wilting*(alpha/alpha_0))^(x*(alpha/alpha_0))K(S(alpha_0/alpha*h))/vg_Ks dh + vg.fast_mfp[sp](x) = alpha_0/alpha*vg_Ks * (vg.fast_mfp[sp_base](x*(alpha/alpha_0))-vg.fast_mfp[sp_base](h_wilt*(alpha/alpha_0))) + vg.fast_mfp[sp](x)-vg.fast_mfp[sp](sx) = alpha_0/alpha*vg.Ks * (vg.fast_mfp[sp_base](x*(alpha/alpha_0))-vg.fast_mfp[sp_base](sx*(alpha/alpha_0))) + sp_base is the van genuchten parameter set of Ks = 1 and alpha = alpha_0. Since sx*(alpha/alpha_0) should ideally be above the wilting point, choose alpha_0 large + + vg.fast_mfp[sp_base](x) + (inner_kr * alpha)/(alpha_0 * vg_Ks * b) * x - ((inner_kr * alpha)/(alpha_0 * Ks * b) * rx + vg.fast_mfp[sp_base](sx)) = 0 + vg.fast_mfp[sp_base](x) + inner_kr_b * x + base_mfp = 0 + + x only depends on 3 variables: inner_kr_b, base_mfp, vg_m (vg.fast_mfp[sp_base](x) depends only on x and this van Genuchten parameter) + + note: this function might not be particular fast as the creation of the sp_base parameter set might take a while, + + """ + + try: + sx = np.array(sx) + mask = inner_kr_ == 0 + inner_kr_[mask] = 1.0e-7 + + + #compute two inputs for the lookup table + rho = np.array(rho_) + rho2 = np.multiply(rho,rho) + b = 2 * (rho2 - 1) / (1 - 0.53 * 0.53 * rho2 + 2 * rho2 * (np.log(rho) + np.log(0.53))) # Vanderborgth et al. 2023, Eqn [8] + inner_kr_b = np.divide(inner_kr_,self.sp.Ksat*b) + base_mfp = - ( inner_kr_b * rx * self.sp.alpha / self.alpha_0 + self.lookup_global_mfp((self.sp.m,sx*self.sp.alpha/self.alpha_0))) + #lookup table + rsx = np.array([self.global_lookup_table((self.sp.m, inner_kr_b[i], base_mfp[i])) for i in range(len(inner_kr_b))])*self.alpha_0/self.sp.alpha + rsx[mask] = sx[mask] # if inner_kr is zero, there is no flow, and the interface potential is the same as the soil potential + except: + if np.max(rx) > 0: + print("xylem matric potential positive", np.max(rx), "at", np.argmax(rx)) + if np.min(rx) < -16000: + print("xylem matric potential under -16000 cm", np.min(rx), "at", np.argmin(rx)) + if np.max(sx) > 0: + print("soil matric potential positive", np.max(sx), "at", np.argmax(sx)) + if np.min(sx) < -16000: + print("soil matric potential under -16000 cm", np.min(sx), "at", np.argmin(sx)) + if np.min(inner_kr_) < 1.0e-7: + print("radius times radial conductivity below 1.e-7", np.min(inner_kr_), "at", np.argmin(inner_kr_)) + if np.max(inner_kr_) > 1.0e-4: + print("radius times radial conductivity above 1.e-4", np.max(inner_kr_), "at", np.argmax(inner_kr_)) + if np.min(rho_) < 1: + print("geometry factor below 1", np.min(rho_), "at", np.argmin(rho_)) + if np.max(rho_) > 200: + print("geometry factor above 200", np.max(rho_), "at", np.argmax(rho_)) + + print("PerirhizalPython.soil_root_interface_potentials_table(): table look up failed, value exceeds table") + print( + "\trx", + np.min(rx), + np.max(rx), + "sx", + np.min(sx), + np.max(sx), + "inner_kr", + np.min(inner_kr_), + np.max(inner_kr_), + "rho", + np.min(rho_), + np.max(rho_), + ) + raise + + return rsx def get_density(self, type: str, volumes=None): """retrieves length, surface or volume density [cm/cm3, cm2/cm3, or cm3/cm3] per soil cells @@ -683,11 +1610,23 @@ def get_default_volumes_(self): # sand = [0.045, 0.43, 0.15, 3, 1000] # loam = [0.08, 0.43, 0.04, 1.6, 50] # clay = [0.1, 0.4, 0.01, 1.1, 10] - - # hydrus_loam = [0.078, 0.43, 0.036, 1.56, 24.96] - # hydrus_clay = [0.068, 0.38, 0.008, 1.09, 4.8] - # hydrus_sand = [0.045, 0.43, 0.145, 2.68, 712.8] - # hydrus_sandyloam = [0.065, 0.41, 0.075, 1.89, 106.1] + hydrus_loam = [0.078, 0.43, 0.036, 1.56, 24.96] + filename = "hydrus_loam" + sp = vg.Parameters(hydrus_loam) + vg.create_mfp_lookup(sp) + peri = PerirhizalPython() + #peri.create_lookup_mpi(filename, sp) # takes some hours; mpiexec -n 4 python Perirhizal.py + start = time.time() + peri.create_lookup(filename, sp) # should be better + end = time.time() + print("Creating the simplified lookup table for hydrus loam took ", end - start, " seconds.") + start = time.time() + peri.create_lookup_global(peri.water_filename) # this is global + end = time.time() + print("Creating the general lookup table for all van Genuchten sets took ", end - start, " seconds.") + + # generate some tests + # filename = "hydrus_loam" # sp = vg.Parameters(hydrus_loam) diff --git a/src/functional/van_genuchten.py b/src/functional/van_genuchten.py index 61a193278..34b60f675 100644 --- a/src/functional/van_genuchten.py +++ b/src/functional/van_genuchten.py @@ -154,9 +154,10 @@ def matric_potential_mfp(mfp, sp): call create_mfp_lookup first, once for each soil parameter @param sp""" -def create_mfp_lookup(sp, wilting_point = -16000, n = 16001): +def create_mfp_lookup(sp, wilting_point = -16000, n = 16001, verbose = True): """ initializes the look up tables for soil parameter to use fast_mfp, and fast_imfp """ - print("initializing look up tables") + if verbose: + print("initializing mfp look up tables") global fast_mfp global fast_imfp @@ -176,7 +177,8 @@ def create_mfp_lookup(sp, wilting_point = -16000, n = 16001): imfp[i] = h_[i] fast_imfp[sp] = interpolate.interp1d(mfp, imfp, bounds_error = False, fill_value = (imfp[0], imfp[-1])) # - print("done") + if verbose: + print("done") # fast_specific_moisture_storage = {} # fast_water_content = {} diff --git a/tutorial/chapter7_coupled/example7_3_lookup.py b/tutorial/chapter7_coupled/example7_3_lookup.py index 7e575d680..dbe3c4e6c 100644 --- a/tutorial/chapter7_coupled/example7_3_lookup.py +++ b/tutorial/chapter7_coupled/example7_3_lookup.py @@ -9,4 +9,4 @@ filename = "hydrus_loam" sp = vg.Parameters(hydrus_loam) # |\label{l73l:soil_end}| vg.create_mfp_lookup(sp) # |\label{l73l:mfp}| -peri.create_lookup_mpi("results/" + filename, sp) # |\label{l73l:lookup}| +peri.create_lookup("results/" + filename, sp) # |\label{l73l:lookup}| diff --git a/tutorial/chapter7_coupled/results/.gitignore b/tutorial/chapter7_coupled/results/.gitignore index ca49de51a..33498976f 100644 --- a/tutorial/chapter7_coupled/results/.gitignore +++ b/tutorial/chapter7_coupled/results/.gitignore @@ -1,2 +1,4 @@ *.vtk *.vtu +*.npy +*.npz \ No newline at end of file