Skip to content

Commit

Permalink
update 3d fluxes
Browse files Browse the repository at this point in the history
  • Loading branch information
@ftomei
ftomei committed Jan 9, 2025
1 parent 359855d commit 9199078
Show file tree
Hide file tree
Showing 11 changed files with 252 additions and 222 deletions.
2 changes: 1 addition & 1 deletion DATA/PROJECT/Troy/SETTINGS/parameters.ini
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,4 +86,4 @@ conductivityHorizVertRatio=10
freeCatchmentRunoff=true
freeBottomDrainage=true
freeLateralDrainage=true
modelAccuracy=5
modelAccuracy=4
11 changes: 7 additions & 4 deletions agrolib/crop/root.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand All @@ -282,19 +283,21 @@ 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++)
{
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;
Expand Down
106 changes: 62 additions & 44 deletions agrolib/soilFluxes3D/balance.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
Expand All @@ -128,78 +134,88 @@ 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++)
{
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);
Expand All @@ -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);
Expand All @@ -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))
Expand All @@ -272,7 +290,7 @@ bool waterBalance(double deltaT, int approxNr)
}

// default
return (false);
return false;
}


Expand Down
6 changes: 3 additions & 3 deletions agrolib/soilFluxes3D/boundary.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
Expand Down
5 changes: 2 additions & 3 deletions agrolib/soilFluxes3D/header/parameters.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down
20 changes: 10 additions & 10 deletions agrolib/soilFluxes3D/header/soilPhysics.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion agrolib/soilFluxes3D/soilFluxes3D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
Loading

0 comments on commit 9199078

Please sign in to comment.