Oscillations like fortran

.gitignore

@@ -3,5 +3,5 @@ __pycache__
 *.npy
 .*.swp
 Results
-.vscode
+.vscode/*.json
 *.lock

README.md

@@ -52,8 +52,6 @@ python marlpde/Evolve_scenario.py
 ```
 Results in the form of an .hdf5 file will be stored in a subdirectory of a `Results` directory, which will be in the root folder of the cloned repo.
 
-After two minutes you should see plots similar to figure 3e from [L'Heureux (2018)](https://www.hindawi.com/journals/geofluids/2018/4968315/).
-
 ### Alternative: poetry
 If you prefer [`poetry`](https://python-poetry.org/) over `pipenv`, you may install all the dependencies and activate the environment using the command `poetry install`. Next, either:
 

marlpde/Evolve_scenario.py

@@ -16,7 +16,7 @@ def integrate_equations(**kwargs):
     '''
     This function retrieves the parameters of the Scenario to be simulated and 
     the solution parameters for the integration. It then integrates the five
-    partial differential equations form L'Heureux, stores and returns the 
+    partial differential equations from L'Heureux, stores and returns the 
     solution, to be used for plotting.
     '''
 
@@ -93,16 +93,27 @@ def integrate_equations(**kwargs):
     else:
         data_tracker = None
 
-    sol, info = eq.solve(state, t_range=End_time, dt=dt, \
-                         solver=kwargs["solver"], scheme=kwargs["scheme"],\
-                         tracker=["progress", \
-                         storage.tracker(kwargs["progress_tracker_interval"]),\
-                         live_plots, data_tracker], \
-                         backend=kwargs["backend"], ret_info=kwargs["retinfo"],\
-                         adaptive=kwargs["adaptive"])
+    if kwargs["solver"]=='scipy':
+        sol, info = eq.solve(state, t_range=End_time, dt=dt, 
+                             solver='scipy', method=kwargs["method"], 
+                             lband=kwargs["lband"], uband=kwargs["uband"], 
+                             tracker=["progress", storage.tracker(
+                                 kwargs["progress_tracker_interval"]),
+                                 live_plots, data_tracker], 
+                             backend=kwargs["backend"], 
+                             ret_info=kwargs["retinfo"])
+    else:
+        sol, info = eq.solve(state, t_range=End_time, dt=dt, 
+                             solver=kwargs["solver"], scheme=kwargs["scheme"],
+                             tracker=["progress", storage.tracker(
+                                 kwargs["progress_tracker_interval"]),
+                                 live_plots, data_tracker], 
+                             backend=kwargs["backend"], 
+                             ret_info=kwargs["retinfo"],\
+                             adaptive=kwargs["adaptive"])
 
     print()
-    print(f"Meta-information about the solution : {info}")        
+    print(f"Meta-information about the solution : {info}\n")        
 
     covered_time_span = Tstar * info["controller"]["t_final"]
 

marlpde/LHeureux_model.py

@@ -16,15 +16,8 @@ def __init__(self, AragoniteSurface, CalciteSurface, CaSurface,
         self.CaSurface = CaSurface
         self.CO3Surface = CO3Surface
         self.PorSurface = PorSurface
-        # There are no bottom boundary conditions for CA and CC in equations 
-        # 35 from L'Heureux, but py-pde demands left and right boundary 
-        # conditions. This means that '"curvature" : 0' is a dummy boundary
-        # condition. We will make sure that it does not leak into our 
-        # integration, by using only backward differencing for the spatial
-        # derivatives in the right-hand sides of the time derivative equations
-        # for CA and CC.
-        self.bc_CA = [{"value": CA0}, {"value" : 1e99}]
-        self.bc_CC = [{"value": CC0}, {"value": 1e99}]
+        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}]
@@ -87,6 +80,7 @@ def __init__(self, AragoniteSurface, CalciteSurface, CaSurface,
         self.F_fixed = 1 - np.exp(10 - 10 / self.PhiIni)
         self.dPhi_fixed = self.auxcon * self.F_fixed *\
                           self.PhiIni ** 3 / (1 - self.PhiIni) 
+
 
     def get_state(self, AragoniteSurface, CalciteSurface, CaSurface, CO3Surface, 
                   PorSurface):
@@ -167,13 +161,18 @@ def evolution_rate(self, state, t=0):
         one_minus_Phi = 1 - Phi
         U = self.presum + self.rhorat * Phi ** 3 * F /one_minus_Phi
 
-        # Enforce no bottom boundary condition for CA and CC by using
-        # backwards differencing only.
+        # Choose either forward or backward differencing for CA and CC
+        # depending on the sign of U
+
         CA_grad_back = CA.apply_operator("grad_back", self.bc_CA)
-        CA_grad = CA_grad_back
+        CA_grad_forw = CA.apply_operator("grad_forw", self.bc_CA)
+        CA_grad = ScalarField(state.grid, np.where(U.data>0, CA_grad_back.data, \
+            CA_grad_forw.data))
 
         CC_grad_back = CC.apply_operator("grad_back", self.bc_CC)
-        CC_grad = CC_grad_back
+        CC_grad_forw = CC.apply_operator("grad_forw", self.bc_CC)
+        CC_grad = ScalarField(state.grid, np.where(U.data>0, CC_grad_back.data, \
+            CC_grad_forw.data))
 
         W = self.presum - self.rhorat * Phi ** 2 * F
 
@@ -269,11 +268,10 @@ def _make_pde_rhs_numba(self, state):
         Peclet_min = self.Peclet_min
         Peclet_max = self.Peclet_max
         delta_x = state.grid._axes_coords[0][1] - state.grid._axes_coords[0][0]
-
-        # Enforce no bottom boundary condition for CA and CC by using
-        # backwards differencing only.
         grad_back_CA = state.grid.make_operator("grad_back", bc = self.bc_CA)
+        grad_forw_CA = state.grid.make_operator("grad_forw", bc = self.bc_CA)
         grad_back_CC = state.grid.make_operator("grad_back", bc = self.bc_CC)
+        grad_forw_CC = state.grid.make_operator("grad_forw", bc = self.bc_CC)
         grad_back_cCa = state.grid.make_operator("grad_back", bc = self.bc_cCa)
         grad_forw_cCa = state.grid.make_operator("grad_forw", bc = self.bc_cCa)
         laplace_cCa = state.grid.make_operator("laplace", bc = self.bc_cCa)
@@ -292,9 +290,11 @@ def pde_rhs(state_data, t=0):
             # either one of them, based on the sign of U.
             CA = state_data[0]
             CA_grad_back = grad_back_CA(CA)
+            CA_grad_forw = grad_forw_CA(CA)
 
             CC = state_data[1]
             CC_grad_back = grad_back_CC(CC)
+            CC_grad_forw = grad_forw_CC(CC)
 
             cCa = state_data[2]
             cCa_grad_back = grad_back_cCa(cCa)
@@ -348,10 +348,12 @@ def pde_rhs(state_data, t=0):
 
                 U[i] = presum + rhorat * Phi[i] ** 3 * F[i]/ (1 - Phi[i])
 
-                # Enforce no bottom boundary condition for CA and CC by using
-                # backwards differencing only.
-                CA_grad[i] = CA_grad_back[i]
-                CC_grad[i] = CC_grad_back[i]
+                if U[i] > 0:
+                    CA_grad[i] = CA_grad_back[i]
+                    CC_grad[i] = CC_grad_back[i]
+                else:
+                    CA_grad[i] = CA_grad_forw[i]
+                    CC_grad[i] = CC_grad_forw[i]
 
                 W[i] = presum - rhorat * Phi[i] ** 2 * F[i]
 

marlpde/parameters.py

@@ -1,6 +1,3 @@
-# ~\~ language=Python filename=marlpde/marlpde.py
-# ~\~ begin <<lit/python-interface.md|marlpde/marlpde.py>>[init]
-# import numpy
 import configparser
 from pathlib import Path
 from dataclasses import (dataclass, asdict, make_dataclass, fields)
@@ -39,17 +36,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.6 * 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)
-    xcem: quantity   = -100.0 * u.cm
-    xcemf: 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.5 * 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
@@ -155,24 +150,32 @@ class Solver:
     So parameters like time interval, time step and tolerance.
     '''
     dt: float     = 1.e-6
-    # T* is more suitable as a default value
-    # than the original value (100_000 years).
-    tmax: float     = Map_Scenario().Tstar
-    # timesteps in between writing.
+    # tmax is the integration time in years.
+    tmax: int     = Map_Scenario().Tstar
     N: int        = 200
-    solver: str   = "explicit"
-    scheme: str   = "rk"
+    solver: str   = "scipy"
+    # Beware that "scheme" and "adaptive" will only be propagated if you have 
+    # chosen py-pde's native "explicit" solver above.
+    scheme: str   = "euler"
+    adaptive: bool = True
+    # solve_ivp from scipy offers six methods. They can be set here.
+    method: str = "LSODA"
+    # Setting lband and uband for method="LSODA" leads to tremendous performance
+    # increase. See Scipy's solve_ivp documentation for background. Consider it
+    # equivalent to providing a sparsity matrix for the "Radau" and "BDF"
+    # implicit methods.
+    lband: int = 1
+    uband: int = 1
     backend: str = "numba"
     retinfo: bool = True
-    adaptive: bool = True
 
 @dataclass
 class Tracker:
     '''
     Initialises all the tracking parameters, such as tracker interval.
     Also indicates the quantities to be tracked, as boolean values.
     '''
-    progress_tracker_interval: float = 0.01
+    progress_tracker_interval: float = Solver().tmax/ 100 / Map_Scenario().Tstar
     live_plotting: bool = False
     plotting_interval: str = '0:05'
     data_tracker_interval: float = 0.01

tests/Regression_test/test_regression.py

@@ -38,10 +38,10 @@ def test_integration_Scenario_A():
     # dict containing the solver parameters (such as required tolerance).
     # The ground truth data were generated with PhiNR=Phi0. In this codebase
     # that was later changed to PhiNR=PhiIni, to comply with the L'Heureux
-    # paper, so we need to reapply the old condition PhiNR=Phi0 to arrive at the
-    # ground truth data.
+    # paper, so we need to reapply the old condition PhiNR=Phi0 to arrive at 
+    # the ground truth data.
     all_kwargs = asdict(Map_Scenario()) | asdict(Solver()) | \
-                 asdict(Tracker()) | {"PhiNR": 0.6}
+                 asdict(Tracker()) | {"Phi0": 0.6, "PhiIni": 0.5, "PhiNR": 0.6}
     # integrate_equations returns four variables, we only need the first one.
     solution, _, _, _, _ = \
         integrate_equations(**all_kwargs)