This should be (one of) the last commit(s) before merging with main. …

…It includes removing the typo; any mainpulation of the porosity diffusion coefficient should be done through modifying b. Canonical values of phi0 and phi00 are now 0.8 and 0.8, respectively. These differ from Scenario A, but we are looking for oscillations. We will investigating larger grids and a higher grid resolution, but we will set for 5m and 2.5cm as default values. The default ad-hoc bottom boundary conditions for aragonite and calcite were chosen to align with rhythmite. This means that the spatial derivatives of the deepest and second deepest node are equal. If two adjacent spatial derivatives are equal, this means that the second spatial derivative is equal to 0. This is reflected by {'curvature': 0}.

marlpde/LHeureux_model.py

@@ -16,8 +16,8 @@ def __init__(self, AragoniteSurface, CalciteSurface, CaSurface,
         self.CaSurface = CaSurface
         self.CO3Surface = CO3Surface
         self.PorSurface = PorSurface
-        self.bc_CA = [{"value": CA0}, {"curvature" : 0}]
-        self.bc_CC = [{"value": CC0}, {"curvature": 0}]
+        self.bc_CA = [{"value": CA0}, {"curvature": 0}]
+        self.bc_CC = [{"value": CC0}, {"curvature": 0}]        
         self.bc_cCa = [{"value": cCa0}, {"derivative": 0}]
         self.bc_cCO3 = [{"value": cCO30}, {"derivative": 0}]
         self.bc_Phi = [{"value": Phi0}, {"derivative": 0}]
@@ -79,7 +79,7 @@ def __init__(self, AragoniteSurface, CalciteSurface, CaSurface,
         self.PhiIni = PhiIni
         self.F_fixed = 1 - np.exp(10 - 10 / self.PhiIni)
         self.dPhi_fixed = self.auxcon * self.F_fixed *\
-                          self.PhiIni * 3 / (1 - self.PhiIni) 
+                          self.PhiIni ** 3 / (1 - self.PhiIni) 
 
 
     def get_state(self, AragoniteSurface, CalciteSurface, CaSurface, CO3Surface, 

marlpde/parameters.py

@@ -39,15 +39,15 @@ class Scenario:
     S: quantity      = 0.1 * u.cm / u.a
     # cAthy: quantity  = 0.1 * u.dimensionless
     phiinf: quantity = 0.01 * u.dimensionless
-    phi0: quantity   = 0.9 * u.dimensionless
+    phi0: quantity   = 0.8 * u.dimensionless
     ca0: quantity    = 0.326e-3 * u.M
     co30: quantity   = 0.326e-3 * u.M
     ccal0: quantity  = 0.3 * u.dimensionless
     cara0: quantity  = 0.6 * u.dimensionless
     xdis: quantity   = 50.0 * u.cm       # x_d   (start of dissolution zone)
-    length: quantity = 1000.0 * u.cm
+    length: quantity = 500.0 * u.cm
     Th: quantity     = 100.0 * u.cm      # h_d   (height of dissolution zone)
-    phi00: quantity  = 0.9 * u.dimensionless
+    phi00: quantity  = 0.8 * u.dimensionless
     ca00: quantity   = 0.326e-3 * u.M    # sqrt(Kc) / 2
     co300: quantity  = 0.326e-3 * u.M    # sqrt(Kc) / 2
     ccal00: quantity = 0.3 * u.dimensionless
@@ -157,7 +157,7 @@ class Solver:
     # T* is more suitable as a default value
     # than the original value (100_000 years).
     tmax: int     = 250_000
-    N: int        = 400
+    N: int        = 200
     solver: str   = "scipy"
     scheme: str   = "euler"
     backend: str = "numba"