numerical solving in C

emf-creaf · Oct 23, 2024 · 9373462 · 9373462
1 parent 92184a0
commit 9373462
Show file tree

Hide file tree

Showing 4 changed files with 23 additions and 13 deletions.
diff --git a/src/hydrology.cpp b/src/hydrology.cpp
@@ -1053,8 +1053,9 @@ NumericVector soilWaterBalance(DataFrame soil, String soilFunctions,
               d[l] = (lambda[l]*C_step_m[l]*dZ_m[l]/halftsubstep)*Psi_step_m[l] + drain_above[l] - drain_below[l] + capill_below[l] + finalSourceSinks_m3s[l];
             }
           }
-          NumericVector Psi_step_m_t05 = tridiagonalSolving(a,b,c,d, nlayers);
-          NumericVector Psi_step_t05 = Psi_step_m_t05*mTOMPa; // m to MPa
+          double* Psi_step_m_t05 = tridiagonalSolving(a,b,c,d, nlayers);
+          double* Psi_step_t05 = new double[nlayers];
+          for(int l=0;l<nlayers;l++) Psi_step_t05[l] = Psi_step_m_t05[l]*mTOMPa; // m to MPa
           //Calculate K and C at t05
           for(int l=0;l<nlayers;l++) {
             if(soilDomains=="single") {
@@ -1120,8 +1121,10 @@ NumericVector soilWaterBalance(DataFrame soil, String soilFunctions,
               d[l] = (lambda[l]*C_step_m05[l]*dZ_m[l]/tsubstep)*Psi_step_m[l] + a[l]*(Psi_step_m[l - 1] - Psi_step_m[l]) + capill_below[l] + drain_above[l] - drain_below[l] + finalSourceSinks_m3s[l];
             }
           }
-          NumericVector Psi_step_m_t1 = tridiagonalSolving(a,b,c,d, nlayers);
-          NumericVector Psi_step_t1 = Psi_step_m_t1*mTOMPa; // m to MPa
+          double* Psi_step_m_t1 = tridiagonalSolving(a,b,c,d, nlayers);
+          double* Psi_step_t1 = new double[nlayers];
+          for(int l=0;l<nlayers;l++) Psi_step_t1[l] = Psi_step_m_t1[l]*mTOMPa; // m to MPa
+
           //calculate drainage (m3)
           if(freeDrainage) {
             drainage_matrix_step_m3 += drain_below[nlayers -1]*tsubstep;

diff --git a/src/numerical_solving.cpp b/src/numerical_solving.cpp
@@ -3,10 +3,10 @@
 #include <Rcpp.h>
 using namespace Rcpp;
 
-NumericVector tridiagonalSolving(double* a, double* b, double* c, double* d, int n) {
+double* tridiagonalSolving(double* a, double* b, double* c, double* d, int n) {
   double* e = new double[n];
   double* f = new double[n];
-  NumericVector u(n);
+  double* u = new double[n];
 
   //Forward steps
   double e_prev = 0.0;

diff --git a/src/numerical_solving.h b/src/numerical_solving.h
@@ -5,4 +5,4 @@
 #endif
 using namespace Rcpp;
 
-NumericVector tridiagonalSolving(double* a, double* b, double* c, double* d, int n);
+double* tridiagonalSolving(double* a, double* b, double* c, double* d, int n);
diff --git a/src/soil_thermodynamics.cpp b/src/soil_thermodynamics.cpp
@@ -185,11 +185,15 @@ NumericVector temperatureChange(NumericVector widths, NumericVector Temp,
                                           Temp);
   int nlayers = Temp.length();
 
-  NumericVector dZ_m = widths*0.001; //mm to m
-
   //Estimate layer interfaces
-  NumericVector dZUp(nlayers), dZDown(nlayers), Zcent(nlayers), Zup(nlayers), Zdown(nlayers);
+  double* dZ_m = new double[nlayers];
+  double* dZUp = new double[nlayers];
+  double* dZDown = new double[nlayers];
+  double* Zcent = new double[nlayers];
+  double* Zup = new double[nlayers];
+  double* Zdown = new double[nlayers];
   for(int l=0;l<nlayers;l++) {
+    dZ_m[l] = widths[l]*0.001; //mm to m
     if(l==0) { //first layer
       dZUp[l] = dZ_m[0]/2.0; //Distance from ground to mid-layer
       Zcent[l] = -1.0*dZ_m[0]/2.0; //Center of the layer (negative downwards)
@@ -207,7 +211,8 @@ NumericVector temperatureChange(NumericVector widths, NumericVector Temp,
       dZDown[l] = dZ_m[l]/2.0;
     }
   }
-  NumericVector k_up(nlayers), k_down(nlayers);
+  double* k_up = new double[nlayers];
+  double* k_down = new double[nlayers];
   double* a = new double[nlayers];
   double* b = new double[nlayers];
   double* c = new double[nlayers];
@@ -233,8 +238,10 @@ NumericVector temperatureChange(NumericVector widths, NumericVector Temp,
       d[l] = (k_up[l]/dZUp[l])*(Temp[l-1] - Temp[l]);
     }
   }
-  NumericVector tempch = tridiagonalSolving(a,b,c,d, nlayers);
-  return(tempch);
+  double* tempch = tridiagonalSolving(a,b,c,d, nlayers);
+  NumericVector out(nlayers);
+  for(int l=0;l<nlayers;l++) out[l] = tempch[l];
+  return(out);
 }
 // NumericVector temperatureChange(NumericVector widths, NumericVector Temp,
 //                                 NumericVector sand, NumericVector clay,