From 9c872946a1f02341925a240906ed2ac80ee60a91 Mon Sep 17 00:00:00 2001 From: ftomei Date: Thu, 9 Jan 2025 20:19:46 +0100 Subject: [PATCH] update 3d fluxes --- crop/root.cpp | 11 ++- soilFluxes3D/balance.cpp | 106 +++++++++++--------- soilFluxes3D/boundary.cpp | 6 +- soilFluxes3D/header/parameters.h | 5 +- soilFluxes3D/header/soilPhysics.h | 20 ++-- soilFluxes3D/soilFluxes3D.cpp | 2 +- soilFluxes3D/soilPhysics.cpp | 159 +++++++++++++++--------------- soilFluxes3D/solver.cpp | 26 ++--- soilFluxes3D/water.cpp | 133 ++++++++++++++----------- 9 files changed, 249 insertions(+), 219 deletions(-) diff --git a/crop/root.cpp b/crop/root.cpp index 17f9f1dd3..003045986 100644 --- a/crop/root.cpp +++ b/crop/root.cpp @@ -261,12 +261,13 @@ namespace root lunetteDensity.resize(nrLayersWithRoot*2); double sinAlfa, cosAlfa, alfa; + double halfPI = PI / 2.; for (i = 0; i < nrLayersWithRoot; i++) { sinAlfa = 1 - double(i+1) / double(nrLayersWithRoot); - cosAlfa = MAXVALUE(sqrt(1 - pow(sinAlfa,2)), 0.0001); + cosAlfa = MAXVALUE(sqrt(1 - pow(sinAlfa, 2)), 0.0001); alfa = atan(sinAlfa/cosAlfa); - lunette[i] = ((PI/2) - alfa - sinAlfa*cosAlfa) / PI; + lunette[i] = (halfPI - alfa - sinAlfa*cosAlfa) / PI; } lunetteDensity[0] = lunette[0]; @@ -282,8 +283,7 @@ namespace root LiMin = -log(0.2) / nrLayersWithRoot; Limax = -log(0.05) / nrLayersWithRoot; - // TODO verify - k = LiMin + (Limax - LiMin) * (shapeFactor-1); + k = LiMin + (Limax - LiMin) * (shapeFactor-1); // TODO verify rootDensitySum = 0 ; for (i = 0; i < (2*nrLayersWithRoot); i++) @@ -291,10 +291,13 @@ namespace root lunetteDensity[i] *= exp(-k*(i+0.5)); rootDensitySum += lunetteDensity[i]; } + + // normalize for (i = 0; i < (2*nrLayersWithRoot); i++) { lunetteDensity[i] /= rootDensitySum; } + for (i = 0; i < totalLayers; i++) { densityThinLayers[i] = 0; diff --git a/soilFluxes3D/balance.cpp b/soilFluxes3D/balance.cpp index 439718fe6..bf265344b 100644 --- a/soilFluxes3D/balance.cpp +++ b/soilFluxes3D/balance.cpp @@ -77,48 +77,54 @@ void InitializeBalanceWater() balanceWholePeriod.waterMBE = 0.; balanceWholePeriod.waterMBR = 0.; - /*! initialize link flow */ + // initialize link flow for (long n = 0; n < myStructure.nrNodes; n++) - { + { nodeList[n].up.sumFlow = 0.; nodeList[n].down.sumFlow = 0.; for (short i = 0; i < myStructure.nrLateralLinks; i++) + { nodeList[n].lateral[i].sumFlow = 0.; } + } - /*! initialize boundary flow */ + // initialize boundary flow for (long n = 0; n < myStructure.nrNodes; n++) + { if (nodeList[n].boundary != nullptr) nodeList[n].boundary->sumBoundaryWaterFlow = 0.; + } } /*! - * \brief computes total water content [m^3] - * \return result + * \brief computes the total water content + * \return total water content [m3] */ double computeTotalWaterContent() { - double theta, sum = 0.0; + double sum = 0.0; for (unsigned long i = 0; i < unsigned(myStructure.nrNodes); i++) - if (nodeList[i].isSurface) + { + if (nodeList[i].isSurface) { sum += (nodeList[i].H - double(nodeList[i].z)) * nodeList[i].volume_area; } else { - theta = theta_from_Se(i); - sum += theta * nodeList[i].volume_area; + sum += theta_from_Se(i) * nodeList[i].volume_area; } - return(sum); + } + + return sum; } /*! - * \brief computes sum of water sink/source [m^3] - * \param deltaT - * \return result + * \brief computes sum of water sink/source + * \param deltaT [s] + * \return sum of water sink/source [m3] */ double sumWaterFlow(double deltaT) { @@ -128,69 +134,77 @@ double sumWaterFlow(double deltaT) if (nodeList[n].Qw != 0.) sum += nodeList[n].Qw * deltaT; } - return (sum); + return sum; } - +/*! + * \brief computes the mass balance error in balanceCurrentTimeStep + * \param deltaT [s] + */ void computeMassBalance(double deltaT) { - balanceCurrentTimeStep.storageWater = computeTotalWaterContent(); + balanceCurrentTimeStep.storageWater = computeTotalWaterContent(); // [m3] - double dStorage = balanceCurrentTimeStep.storageWater - balancePreviousTimeStep.storageWater; + double dStorage = balanceCurrentTimeStep.storageWater - balancePreviousTimeStep.storageWater; // [m3] - balanceCurrentTimeStep.sinkSourceWater = sumWaterFlow(deltaT); + balanceCurrentTimeStep.sinkSourceWater = sumWaterFlow(deltaT); // [m3] - balanceCurrentTimeStep.waterMBE = dStorage - balanceCurrentTimeStep.sinkSourceWater; + balanceCurrentTimeStep.waterMBE = dStorage - balanceCurrentTimeStep.sinkSourceWater; // [m3] - /*! reference water: sumWaterFlow or 0.1% of storage */ - double denominator = MAXVALUE(fabs(balanceCurrentTimeStep.sinkSourceWater), balanceCurrentTimeStep.storageWater * 1e-3); + // minimum reference water storage: 0.1% of current storage, minimum 1 liter + double minRefWaterStorage = MAXVALUE(balanceCurrentTimeStep.storageWater * 1e-3, 0.001); // [m3] - /*! no water - minimum 1 liter */ - denominator = MAXVALUE(denominator, 0.001); + // reference water for computation of mass balance ratio + double referenceWater = MAXVALUE(fabs(balanceCurrentTimeStep.sinkSourceWater), minRefWaterStorage); // [m3] - balanceCurrentTimeStep.waterMBR = balanceCurrentTimeStep.waterMBE / denominator; + balanceCurrentTimeStep.waterMBR = balanceCurrentTimeStep.waterMBE / referenceWater; } double getMatrixValue(long i, TlinkedNode *link) { if (link != nullptr) - { + { int j = 1; - while ((j < myStructure.maxNrColumns) && (A[i][j].index != NOLINK) && (A[i][j].index != (*link).index)) j++; + while ((j < myStructure.maxNrColumns) && (A[i][j].index != NOLINK) && (A[i][j].index != (*link).index)) + j++; /*! Rebuild the A elements (previously normalized) */ if (A[i][j].index == (*link).index) + { return (A[i][j].val * A[i][0].val); } + } + return double(INDEX_ERROR); } /*! - * \brief updating of in- and out-flows [m^3] + * \brief updates in and out flows [m3] * \param index - * \param link TlinkedNode pointer - * \param delta_t + * \param link TlinkedNode pointer + * \param delta_t [s] */ void update_flux(long index, TlinkedNode *link, double delta_t) { if (link->index != NOLINK) - (*link).sumFlow += float(getWaterExchange(index, link, delta_t)); + { + (*link).sumFlow += float(getWaterExchange(index, link, delta_t)); // [m3] + } } - void saveBestStep() { for (long n = 0; n < myStructure.nrNodes; n++) + { nodeList[n].bestH = nodeList[n].H; + } } - - void restoreBestStep(double deltaT) { for (unsigned long n = 0; n < unsigned(myStructure.nrNodes); n++) @@ -198,8 +212,10 @@ void restoreBestStep(double deltaT) nodeList[n].H = nodeList[n].bestH; /*! compute new soil moisture (only sub-surface nodes) */ - if (!nodeList[n].isSurface) - nodeList[n].Se = computeSe(n); + if (! nodeList[n].isSurface) + { + nodeList[n].Se = computeSe(n); + } } computeMassBalance(deltaT); @@ -215,18 +231,20 @@ void acceptStep(double deltaT) /*! update sum of flow */ for (long i = 0; i < myStructure.nrNodes; i++) - { + { update_flux(i, &(nodeList[i].up), deltaT); update_flux(i, &(nodeList[i].down), deltaT); for (short j = 0; j < myStructure.nrLateralLinks; j++) + { update_flux(i, &(nodeList[i].lateral[j]), deltaT); + } if (nodeList[i].boundary != nullptr) nodeList[i].boundary->sumBoundaryWaterFlow += nodeList[i].boundary->waterFlow * deltaT; - } - + } } + bool waterBalance(double deltaT, int approxNr) { computeMassBalance(deltaT); @@ -243,15 +261,15 @@ bool waterBalance(double deltaT, int approxNr) /*! best case */ if (MBRerror < myParameters.MBRThreshold) - { + { acceptStep(deltaT); - if ((approxNr < 2) && (Courant < 0.5) && (MBRerror < (myParameters.MBRThreshold * 0.5))) - { + if ((approxNr <= 2) && (Courant < 0.5) && (MBRerror < (myParameters.MBRThreshold * 0.5))) + { /*! system is stable: double time step */ doubleTimeStep(); - } - return (true); } + return true; + } /*! worst case: error high or last approximation */ if ((MBRerror > (bestMBRerror * 2.0)) @@ -272,7 +290,7 @@ bool waterBalance(double deltaT, int approxNr) } // default - return (false); + return false; } diff --git a/soilFluxes3D/boundary.cpp b/soilFluxes3D/boundary.cpp index e39c7f2d4..d571d4d35 100644 --- a/soilFluxes3D/boundary.cpp +++ b/soilFluxes3D/boundary.cpp @@ -285,9 +285,9 @@ void updateBoundaryWater (double deltaT) else if (nodeList[i].boundary->type == BOUNDARY_PRESCRIBEDTOTALPOTENTIAL) { // water table - double L = 1.0; // [m] - double boundaryZ = nodeList[i].z - L; // [m] - double boundaryK; // [m s-1] + double L = 1.0; // [m] + double boundaryZ = nodeList[i].z - L; // [m] + double boundaryK; // [m s-1] if (nodeList[i].boundary->prescribedTotalPotential >= boundaryZ) { diff --git a/soilFluxes3D/header/parameters.h b/soilFluxes3D/header/parameters.h index 37de879f9..17fa0a307 100644 --- a/soilFluxes3D/header/parameters.h +++ b/soilFluxes3D/header/parameters.h @@ -13,8 +13,7 @@ double delta_t_min; double delta_t_max; double current_delta_t; - int iterazioni_min; - int iterazioni_max; + int maxIterationsNumber; int maxApproximationsNumber; int waterRetentionCurve; int meanType; @@ -27,7 +26,7 @@ delta_t_min = 1; delta_t_max = 600; current_delta_t = delta_t_max; - iterazioni_max = 200; + maxIterationsNumber = 200; maxApproximationsNumber = 10; MBRThreshold = 1E-6; ResidualTolerance = 1E-10; diff --git a/soilFluxes3D/header/soilPhysics.h b/soilFluxes3D/header/soilPhysics.h index 1cc56c25d..d48b1ddc8 100644 --- a/soilFluxes3D/header/soilPhysics.h +++ b/soilFluxes3D/header/soilPhysics.h @@ -4,16 +4,16 @@ struct Tsoil; double computeWaterConductivity(double Se, Tsoil *mySoil); - double computeSefromPsi_unsat(double myPsi, Tsoil *mySoil); - double theta_from_Se(unsigned long myIndex); - double theta_from_Se (double Se, unsigned long myIndex); - double theta_from_sign_Psi (double myPsi, unsigned long myIndex); - double Se_from_theta (unsigned long myIndex, double myTheta); - double psi_from_Se(unsigned long myIndex); - double computeSe(unsigned long myIndex); - double dTheta_dH(unsigned long myIndex); - double dThetav_dH(unsigned long myIndex, double temperature, double dTheta_dH); - double computeK(unsigned long myIndex); + double computeSefromPsi_unsat(double psi, Tsoil *mySoil); + double theta_from_Se(unsigned long index); + double theta_from_Se (double Se, unsigned long index); + double theta_from_sign_Psi (double myPsi, unsigned long index); + double Se_from_theta (unsigned long index, double myTheta); + double psi_from_Se(unsigned long index); + double computeSe(unsigned long index); + double dTheta_dH(unsigned long index); + double dThetav_dH(unsigned long index, double temperature, double dTheta_dH); + double computeK(unsigned long index); double compute_K_Mualem(double Ksat, double Se, double VG_Sc, double VG_m, double Mualem_L); double getThetaMean(long i); double getTheta(long i, double H); diff --git a/soilFluxes3D/soilFluxes3D.cpp b/soilFluxes3D/soilFluxes3D.cpp index b3a9a9790..3783c2908 100644 --- a/soilFluxes3D/soilFluxes3D.cpp +++ b/soilFluxes3D/soilFluxes3D.cpp @@ -153,7 +153,7 @@ int DLL_EXPORT __STDCALL setNumericalParameters(float minDeltaT, float maxDeltaT if (maxIterationNumber < 10) maxIterationNumber = 10; if (maxIterationNumber > MAX_NUMBER_ITERATIONS) maxIterationNumber = MAX_NUMBER_ITERATIONS; - myParameters.iterazioni_max = maxIterationNumber; + myParameters.maxIterationsNumber = maxIterationNumber; if (maxApproximationsNumber < 1) maxApproximationsNumber = 1; diff --git a/soilFluxes3D/soilPhysics.cpp b/soilFluxes3D/soilPhysics.cpp index 291d64e55..90f675caf 100644 --- a/soilFluxes3D/soilPhysics.cpp +++ b/soilFluxes3D/soilPhysics.cpp @@ -37,29 +37,26 @@ /*! * \brief Computes volumetric water content from current degree of saturation - * \param myIndex * \return result */ - double theta_from_Se (unsigned long myIndex) + double theta_from_Se (unsigned long index) { - return ((nodeList[myIndex].Se * (nodeList[myIndex].Soil->Theta_s - nodeList[myIndex].Soil->Theta_r)) + nodeList[myIndex].Soil->Theta_r); + return ((nodeList[index].Se * (nodeList[index].Soil->Theta_s - nodeList[index].Soil->Theta_r)) + nodeList[index].Soil->Theta_r); } /*! * \brief Computes volumetric water content from degree of saturation - * \param Se - * \param myIndex + * \param Se degree of saturation [-] * \return result */ - double theta_from_Se (double Se, unsigned long myIndex) + double theta_from_Se (double Se, unsigned long index) { - return ((Se * (nodeList[myIndex].Soil->Theta_s - nodeList[myIndex].Soil->Theta_r)) + nodeList[myIndex].Soil->Theta_r); + return ((Se * (nodeList[index].Soil->Theta_s - nodeList[index].Soil->Theta_r)) + nodeList[index].Soil->Theta_r); } /*! * \brief Computes volumetric water content from water potential (with sign) * \param signPsi water potential with sign [m] - * \param index * \return volumetric water content [m3 m-3] */ double theta_from_sign_Psi (double signPsi, unsigned long index) @@ -81,44 +78,43 @@ /*! * \brief Computes degree of saturation from volumetric water content - * \param myIndex - * \param theta + * \param theta volumetric water content [m3 m-3] * \return result */ - double Se_from_theta (unsigned long myIndex, double theta) + double Se_from_theta (unsigned long index, double theta) { /*! check range */ - if (theta >= nodeList[myIndex].Soil->Theta_s) return(1.); - else if (theta <= nodeList[myIndex].Soil->Theta_r) return(0.); - else return ((theta - nodeList[myIndex].Soil->Theta_r) / (nodeList[myIndex].Soil->Theta_s - nodeList[myIndex].Soil->Theta_r)); + if (theta >= nodeList[index].Soil->Theta_s) return(1.); + else if (theta <= nodeList[index].Soil->Theta_r) return(0.); + else return ((theta - nodeList[index].Soil->Theta_r) / (nodeList[index].Soil->Theta_s - nodeList[index].Soil->Theta_r)); } /*! - * \brief Computes degree of saturation from matric potential (Van Genutchen) Se = [1+(alfa*|h|)^n]^-m - * valid only for unsaturated soil - * \param myPsi water potential (absolute value) [m] + * \brief Computes degree of saturation from matric potential Se = [1+(alfa*|h|)^n]^-m + * \param psi [m] water potential (absolute value) * \param mySoil * \return degree of saturation [-] */ - double computeSefromPsi_unsat(double myPsi, Tsoil *mySoil) + double computeSefromPsi_unsat(double psi, Tsoil *mySoil) { double Se = NODATA; if (myParameters.waterRetentionCurve == MODIFIEDVANGENUCHTEN) { - if (myPsi <= mySoil->VG_he) + if (psi <= mySoil->VG_he) { + // saturated Se = 1.; } else { - Se = pow(1. + pow(mySoil->VG_alpha * myPsi, mySoil->VG_n), - mySoil->VG_m); + Se = pow(1. + pow(mySoil->VG_alpha * psi, mySoil->VG_n), - mySoil->VG_m); Se *= (1. / mySoil->VG_Sc); } } else if (myParameters.waterRetentionCurve == VANGENUCHTEN) { - Se = pow(1. + pow(mySoil->VG_alpha * myPsi, mySoil->VG_n), - mySoil->VG_m); + Se = pow(1. + pow(mySoil->VG_alpha * psi, mySoil->VG_n), - mySoil->VG_m); } return Se; @@ -126,12 +122,11 @@ /*! * \brief Computes current degree of saturation - * \param myIndex - * \return result + * \return degree of saturation [-] */ - double computeSe(unsigned long myIndex) + double computeSe(unsigned long index) { - if (nodeList[myIndex].H >= nodeList[myIndex].z) + if (nodeList[index].H >= nodeList[index].z) { // saturated return 1.; @@ -139,84 +134,83 @@ else { // unsaturated - double psi = fabs(nodeList[myIndex].H - nodeList[myIndex].z); /*!< [m] */ - return computeSefromPsi_unsat(psi, nodeList[myIndex].Soil); + double psi = fabs(nodeList[index].H - nodeList[index].z); /*!< [m] */ + return computeSefromPsi_unsat(psi, nodeList[index].Soil); } } /*! - * \brief Computes hydraulic conductivity [m/sec] (Mualem) + * \brief Computes hydraulic conductivity, passing a soil structure + * Mualem equation: * K(Se) = Ksat * Se^(L) * {1-[1-Se^(1/m)]^m}^2 * WARNING: very low values are possible (es: 10^12) - * \param Se - * \param mySoil Tsoil pointer - * \return result + * \param Se degree of saturation [-] + * \param mySoil Tsoil pointer + * \return hydraulic conductivity [m/sec] */ double computeWaterConductivity(double Se, Tsoil *mySoil) { if (Se >= 1.) return(mySoil->K_sat ); - double myTmp = NODATA; + double tmp = NODATA; if (myParameters.waterRetentionCurve == MODIFIEDVANGENUCHTEN) { double myNumerator = 1. - pow(1. - pow(Se*mySoil->VG_Sc, 1./mySoil->VG_m), mySoil->VG_m); - myTmp = myNumerator / (1. - pow(1. - pow(mySoil->VG_Sc, 1./mySoil->VG_m), mySoil->VG_m)); + tmp = myNumerator / (1. - pow(1. - pow(mySoil->VG_Sc, 1./mySoil->VG_m), mySoil->VG_m)); } else if (myParameters.waterRetentionCurve == VANGENUCHTEN) - myTmp = 1. - pow(1. - pow(Se, 1./mySoil->VG_m), mySoil->VG_m); + { + tmp = 1. - pow(1. - pow(Se, 1./mySoil->VG_m), mySoil->VG_m); + } - return (mySoil->K_sat * pow(Se, mySoil->Mualem_L) * pow(myTmp , 2.)); + return (mySoil->K_sat * pow(Se, mySoil->Mualem_L) * pow(tmp, 2.)); } /*! - * \brief Computes hydraulic conductivity [m/sec] (Mualem) + * \brief Computes hydraulic conductivity, passing soil parameters + * Mualem equation: * K(Se) = Ksat * Se^(L) * {1-[1-Se^(1/m)]^m}^2 * WARNING: very low values are possible (es: 10^12) - * \param Ksat - * \param Se - * \param VG_Sc - * \param VG_n - * \param VG_m - * \param Mualem_L - * \return result + * \param Se degree of saturation [-] + * \return hydraulic conductivity [m/sec] */ double compute_K_Mualem(double Ksat, double Se, double VG_Sc, double VG_m, double Mualem_L) { if (Se >= 1.) return(Ksat); - double temp= NODATA; + double tmp= NODATA; if (myParameters.waterRetentionCurve == MODIFIEDVANGENUCHTEN) { double num = 1. - pow(1. - pow(Se*VG_Sc, 1./VG_m), VG_m); - temp = num / (1. - pow(1. - pow(VG_Sc, 1./VG_m), VG_m)); + tmp = num / (1. - pow(1. - pow(VG_Sc, 1./VG_m), VG_m)); } else if (myParameters.waterRetentionCurve == VANGENUCHTEN) { - temp = 1. - pow(1. - pow(Se, 1./VG_m), VG_m); + tmp = 1. - pow(1. - pow(Se, 1./VG_m), VG_m); } - return (Ksat * pow(Se, Mualem_L) * pow(temp , 2.)); + return (Ksat * pow(Se, Mualem_L) * pow(tmp , 2.)); } /*! * \brief Computes current soil water total (liquid + vapor) conductivity [m sec^-1] - * \param myIndex + * \param index * \return result */ - double computeK(unsigned long myIndex) + double computeK(unsigned long index) { - double k = compute_K_Mualem(nodeList[myIndex].Soil->K_sat, nodeList[myIndex].Se, - nodeList[myIndex].Soil->VG_Sc, nodeList[myIndex].Soil->VG_m, - nodeList[myIndex].Soil->Mualem_L); + double k = compute_K_Mualem(nodeList[index].Soil->K_sat, nodeList[index].Se, + nodeList[index].Soil->VG_Sc, nodeList[index].Soil->VG_m, + nodeList[index].Soil->Mualem_L); // vapor isothermal flow if (myStructure.computeHeat && myStructure.computeHeatVapor) { - double avgT = getTMean(myIndex); - double kv = IsothermalVaporConductivity(myIndex, nodeList[myIndex].H - nodeList[myIndex].z, avgT); + double avgT = getTMean(index); + double kv = IsothermalVaporConductivity(index, nodeList[index].H - nodeList[index].z, avgT); // from kg s m-3 to m s-1 kv *= (GRAVITY / WATER_DENSITY); @@ -228,26 +222,26 @@ /*! - * \brief Computes Water Potential from degree of saturation - * \param myIndex - * \return result + * \brief Computes current water potential from degree of saturation + * \param index + * \return water potential [m] */ - double psi_from_Se(unsigned long myIndex) + double psi_from_Se(unsigned long index) { - double m = nodeList[myIndex].Soil->VG_m; + double m = nodeList[index].Soil->VG_m; double temp = NODATA; if (myParameters.waterRetentionCurve == MODIFIEDVANGENUCHTEN) - temp = pow(1./ (nodeList[myIndex].Se * nodeList[myIndex].Soil->VG_Sc) , 1./ m ) - 1.; + temp = pow(1./ (nodeList[index].Se * nodeList[index].Soil->VG_Sc) , 1./ m ) - 1.; else if (myParameters.waterRetentionCurve == VANGENUCHTEN) - temp = pow(1./ nodeList[myIndex].Se, 1./ m ) - 1.; + temp = pow(1./ nodeList[index].Se, 1./ m ) - 1.; - return((1./ nodeList[myIndex].Soil->VG_alpha) * pow(temp, 1./ nodeList[myIndex].Soil->VG_n)); + return((1./ nodeList[index].Soil->VG_alpha) * pow(temp, 1./ nodeList[index].Soil->VG_n)); } /*! * \brief [m-1] dThetaV/dH - * \param myIndex + * \param index * \return derivative of vapor volumetric content with respect to H */ double dThetav_dH(unsigned long i, double temperature, double dTheta_dH) @@ -267,48 +261,49 @@ * \brief [m-1] dTheta/dH (Van Genutchen) * dTheta/dH = dSe/dH (Theta_s-Theta_r) * dSe/dH = -sgn(H-z) alfa n m [1+(alfa|(H-z)|)^n]^(-m-1) (alfa|(H-z)|)^n-1 - * \param myIndex * \return derivative of water volumetric content with respect to H */ - double dTheta_dH(unsigned long myIndex) + double dTheta_dH(unsigned long index) { - double alfa = nodeList[myIndex].Soil->VG_alpha; - double n = nodeList[myIndex].Soil->VG_n; - double m = nodeList[myIndex].Soil->VG_m; + double alfa = nodeList[index].Soil->VG_alpha; + double n = nodeList[index].Soil->VG_n; + double m = nodeList[index].Soil->VG_m; - double psi_abs = fabs(MINVALUE(nodeList[myIndex].H - nodeList[myIndex].z, 0.)); - double psiPrevious_abs = fabs(MINVALUE(nodeList[myIndex].oldH - nodeList[myIndex].z, 0.)); + double currentPsi = fabs(MINVALUE(nodeList[index].H - nodeList[index].z, 0.)); + double previousPsi = fabs(MINVALUE(nodeList[index].oldH - nodeList[index].z, 0.)); - if (myParameters.waterRetentionCurve == MODIFIEDVANGENUCHTEN) + if (myParameters.waterRetentionCurve == VANGENUCHTEN) { // saturated - if ((psi_abs <= nodeList[myIndex].Soil->VG_he) && (psiPrevious_abs <= nodeList[myIndex].Soil->VG_he)) return 0.; + if ((currentPsi == 0.) && (previousPsi == 0.)) + return 0.; } - if (myParameters.waterRetentionCurve == VANGENUCHTEN) + if (myParameters.waterRetentionCurve == MODIFIEDVANGENUCHTEN) { - if ((psi_abs == 0.) && (psiPrevious_abs == 0.)) return 0.; + // saturated + if ((currentPsi <= nodeList[index].Soil->VG_he) && (previousPsi <= nodeList[index].Soil->VG_he)) + return 0.; } double dSe_dH; - if (psi_abs == psiPrevious_abs) + if (currentPsi == previousPsi) { - dSe_dH = alfa * n * m * pow(1. + pow(alfa * psi_abs, n), -(m + 1.)) * pow(alfa * psi_abs, n - 1.); + dSe_dH = alfa * n * m * pow(1. + pow(alfa * currentPsi, n), -(m + 1.)) * pow(alfa * currentPsi, n - 1.); if (myParameters.waterRetentionCurve == MODIFIEDVANGENUCHTEN) { - dSe_dH *= (1. / nodeList[myIndex].Soil->VG_Sc); + dSe_dH *= (1. / nodeList[index].Soil->VG_Sc); } } else { - double theta = computeSefromPsi_unsat(psi_abs, nodeList[myIndex].Soil); - double thetaPrevious = computeSefromPsi_unsat(psiPrevious_abs, nodeList[myIndex].Soil); - double delta_H = nodeList[myIndex].H - nodeList[myIndex].oldH; - dSe_dH = fabs((theta - thetaPrevious) / delta_H); + double theta = computeSefromPsi_unsat(currentPsi, nodeList[index].Soil); + double previousTheta = computeSefromPsi_unsat(previousPsi, nodeList[index].Soil); + dSe_dH = fabs((theta - previousTheta) / (nodeList[index].H - nodeList[index].oldH)); } - return dSe_dH * (nodeList[myIndex].Soil->Theta_s - nodeList[myIndex].Soil->Theta_r); + return dSe_dH * (nodeList[index].Soil->Theta_s - nodeList[index].Soil->Theta_r); } diff --git a/soilFluxes3D/solver.cpp b/soilFluxes3D/solver.cpp index edb53148c..dbcb97363 100644 --- a/soilFluxes3D/solver.cpp +++ b/soilFluxes3D/solver.cpp @@ -100,11 +100,11 @@ TlinkedNode* getLink(long i, long j) } -int calcola_iterazioni_max(int num_approssimazione) +int getMaxIterationsNr(int approximationNr) { - float max_iterazioni = float(myParameters.iterazioni_max) - / float(myParameters.maxApproximationsNumber) * float(num_approssimazione + 1); - return MAXVALUE(20, int(max_iterazioni)); + int maxIterationsNr = int((approximationNr + 1) * (float(myParameters.maxIterationsNumber) + / float(myParameters.maxApproximationsNumber))); + return MAXVALUE(20, maxIterationsNr); } @@ -139,8 +139,10 @@ double GaussSeidelIterationWater(short direction) // surface check (H cannot go below z) if (nodeList[i].isSurface) + { if (newX < double(nodeList[i].z)) newX = double(nodeList[i].z); + } // water potential [m] double psi = fabs(newX - double(nodeList[i].z)); @@ -191,16 +193,15 @@ double GaussSeidelIterationHeat() } -bool GaussSeidelRelaxation (int approximation, double residualTolerance, int process) +bool GaussSeidelRelaxation (int approximation, double toleranceThreshold, int process) { - const double MAX_NORM = 1.0; - double currentNorm = MAX_NORM; - double bestNorm = MAX_NORM; - int iteration = 0; + double currentNorm = 1.0; + double bestNorm = currentNorm; - int maxIterationsNr = calcola_iterazioni_max(approximation); + int maxIterationsNr = getMaxIterationsNr(approximation); + int iteration = 0; - while ((currentNorm > residualTolerance) && (iteration < maxIterationsNr)) + while ((currentNorm > toleranceThreshold) && (iteration < maxIterationsNr)) { if (process == PROCESS_HEAT) currentNorm = GaussSeidelIterationHeat(); @@ -218,7 +219,8 @@ bool GaussSeidelRelaxation (int approximation, double residualTolerance, int pro if (currentNorm > (bestNorm * 10.0)) { - return false; // not converging + // non-convergent system + return false; } else if (currentNorm < bestNorm) { diff --git a/soilFluxes3D/water.cpp b/soilFluxes3D/water.cpp index efc46af41..13bb92d22 100644 --- a/soilFluxes3D/water.cpp +++ b/soilFluxes3D/water.cpp @@ -31,6 +31,7 @@ #include "physics.h" #include "commonConstants.h" +#include "basicMath.h" #include "header/types.h" #include "header/water.h" #include "header/soilPhysics.h" @@ -114,21 +115,28 @@ double runoff(long i, long j, TlinkedNode *link, double deltaT, unsigned long ap } - double infiltration(long sup, long inf, TlinkedNode *link, double deltaT) { - double cellDistance = (nodeList[sup].z - nodeList[inf].z) * 2.0; + double cellDistance = nodeList[sup].z - nodeList[inf].z; - /*! unsaturated */ - if (nodeList[inf].H < nodeList[sup].z) - { + /*! unsaturated */ + if (nodeList[inf].H < nodeList[sup].z) + { /*! surface water content [m] */ - // double surfaceH = (nodeList[sup].H + nodeList[sup].oldH) * 0.5; - double surfaceH = nodeList[sup].H; + //double surfaceH = (nodeList[sup].H + nodeList[sup].oldH) * 0.5; + + double initialSurfaceWater = MAXVALUE(nodeList[sup].oldH - nodeList[sup].z, 0); // [m] + + // precipitation: positive - evaporation: negative + double precOrEvapRate = nodeList[sup].Qw / nodeList[sup].volume_area; // [m s-1] /*! maximum water infiltration rate [m/s] */ - double maxInfiltrationRate = (surfaceH - nodeList[sup].z) / deltaT; - if (maxInfiltrationRate <= 0.0) return(0.0); + double maxInfiltrationRate = initialSurfaceWater / deltaT + precOrEvapRate; // [m s-1] + if (maxInfiltrationRate <= 0) + return 0.; + + double dH = nodeList[sup].H - nodeList[inf].H; + double maxK = maxInfiltrationRate * (cellDistance / dH); // [m s-1] /*! first soil layer: mean between current k and k_sat */ double meanK = computeMean(nodeList[inf].Soil->K_sat, nodeList[inf].k); @@ -146,17 +154,12 @@ double infiltration(long sup, long inf, TlinkedNode *link, double deltaT) } } - double dH = nodeList[sup].H - nodeList[inf].H; - double maxK = maxInfiltrationRate * (cellDistance / dH); - - double k = MINVALUE(meanK , maxK); - + double k = MINVALUE(meanK, maxK); return (k * link->area) / cellDistance; - } - - /*! saturated */ - else + } + else { + /*! saturated */ if (nodeList[inf].boundary != nullptr) { if (nodeList[inf].boundary->type == BOUNDARY_URBAN) @@ -172,7 +175,6 @@ double infiltration(long sup, long inf, TlinkedNode *link, double deltaT) return(nodeList[inf].Soil->K_sat * link->area) / cellDistance; } - } @@ -246,10 +248,7 @@ bool computeFlux(long i, int matrixIndex, TlinkedNode *link, double deltaT, unsi bool waterFlowComputation(double deltaT) { - bool isValidStep; - long i; - double dThetadH, dthetavdh; - double avgTemperature; + bool isValidStep = false; int approximationNr = 0; do @@ -257,46 +256,53 @@ bool waterFlowComputation(double deltaT) Courant = 0.0; if (approximationNr == 0) { - for (i = 0; i < myStructure.nrNodes; i++) + // diagonal indexes + for (int i = 0; i < myStructure.nrNodes; i++) { A[i][0].index = i; } } /*! hydraulic conductivity and theta derivative */ - for (i = 0; i < myStructure.nrNodes; i++) + for (unsigned i = 0; i < unsigned(myStructure.nrNodes); i++) { invariantFlux[i] = 0.; - if (!nodeList[i].isSurface) + if (! nodeList[i].isSurface) { - nodeList[i].k = computeK(unsigned(i)); - dThetadH = dTheta_dH(unsigned(i)); - C[i] = nodeList[i].volume_area * dThetadH; - - // vapor capacity term - if (myStructure.computeHeat && myStructure.computeHeatVapor) - { - avgTemperature = getTMean(i); - dthetavdh = dThetav_dH(unsigned(i), avgTemperature, dThetadH); - C[i] += nodeList[i].volume_area * dthetavdh; - } + nodeList[i].k = computeK(i); + double dThetadH = dTheta_dH(i); + C[i] = nodeList[i].volume_area * dThetadH; + + // vapor capacity term + if (myStructure.computeHeat && myStructure.computeHeatVapor) + { + double avgTemperature = getTMean(i); + double dthetavdh = dThetav_dH(i, avgTemperature, dThetadH); + C[i] += nodeList[i].volume_area * dthetavdh; + } } } // update boundary conditions - // updateBoundaryWater(deltaT); + //updateBoundaryWater(deltaT); /*! computes the matrix elements */ - for (i = 0; i < myStructure.nrNodes; i++) + for (int i = 0; i < myStructure.nrNodes; i++) { short j = 1; - if (computeFlux(i, j, &(nodeList[i].up), deltaT, approximationNr, UP)) j++; + if (computeFlux(i, j, &(nodeList[i].up), deltaT, approximationNr, UP)) + j++; for (short l = 0; l < myStructure.nrLateralLinks; l++) - if (computeFlux(i, j, &(nodeList[i].lateral[l]), deltaT, approximationNr, LATERAL)) j++; - if (computeFlux(i, j, &(nodeList[i].down), deltaT, approximationNr, DOWN)) j++; + { + if (computeFlux(i, j, &(nodeList[i].lateral[l]), deltaT, approximationNr, LATERAL)) + j++; + } + if (computeFlux(i, j, &(nodeList[i].down), deltaT, approximationNr, DOWN)) + j++; /*! closure */ - while (j < myStructure.maxNrColumns) A[i][j++].index = NOLINK; + while (j < myStructure.maxNrColumns) + A[i][j++].index = NOLINK; j = 1; double sum = 0.; @@ -307,47 +313,54 @@ bool waterFlowComputation(double deltaT) j++; } - /*! sum of the diagonal elements */ + /*! diagonal */ A[i][0].val = C[i]/deltaT + sum; - /*! b vector(vector of constant terms) */ + /*! b vector (vector of constant terms) */ b[i] = ((C[i] / deltaT) * nodeList[i].oldH) + nodeList[i].Qw + invariantFlux[i]; /*! preconditioning */ j = 1; - while ((j < myStructure.maxNrColumns) && (A[i][j].index != NOLINK)) - A[i][j++].val /= A[i][0].val; + while ((A[i][j].index != NOLINK) && (j < myStructure.maxNrColumns)) + { + A[i][j++].val /= A[i][0].val; + } b[i] /= A[i][0].val; } - if (Courant > 1.0) - if (deltaT > myParameters.delta_t_min) - { - halveTimeStep(); - setForcedHalvedTime(true); - return false; - } + if (Courant > 1.0 && deltaT > myParameters.delta_t_min) + { + halveTimeStep(); + setForcedHalvedTime(true); + return false; + } if (! GaussSeidelRelaxation(approximationNr, myParameters.ResidualTolerance, PROCESS_WATER)) + { if (deltaT > myParameters.delta_t_min) { halveTimeStep(); setForcedHalvedTime(true); - return (false); + return false; } + } /*! set new potential - compute new degree of saturation */ - for (i = 0; i < myStructure.nrNodes; i++) + for (int i = 0; i < myStructure.nrNodes; i++) { nodeList[i].H = X[i]; - if (!nodeList[i].isSurface) + if (! nodeList[i].isSurface) + { nodeList[i].Se = computeSe(unsigned(i)); + } } - /*! water balance */ + // check water balance isValidStep = waterBalance(deltaT, approximationNr); - if (getForcedHalvedTime()) return (false); - } + + if (getForcedHalvedTime()) + return false; + } while ((! isValidStep) && (++approximationNr < myParameters.maxApproximationsNumber)); return isValidStep;