Merge pull request #42 from astro-turing/Oscillations_like_Fortran

Oscillations like fortran

.gitignore

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

Pipfile

@@ -18,6 +18,7 @@ h5py = "*"
 pint = "*"
 ipython = "*"
 pandas = "*"
+jupyter = "*"
 
 [dev-packages]
 pylint = "*"

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

@@ -12,29 +12,26 @@
 from parameters import Map_Scenario, Solver, Tracker
 from LHeureux_model import LMAHeureuxPorosityDiff
 
-def integrate_equations(**kwargs):
+def integrate_equations(solver_parms, tracker_parms, pde_parms):
     '''
     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.
     '''
 
-    Xstar = kwargs["Xstar"]
-    Tstar = kwargs["Tstar"]
-    max_depth = kwargs["max_depth"]
-    ShallowLimit = kwargs["ShallowLimit"]
-    DeepLimit = kwargs["DeepLimit"]
-    CAIni = kwargs["CAIni"]
-    CCIni = kwargs["CCIni"]
-    cCaIni = kwargs["cCaIni"]
-    cCO3Ini = kwargs["cCO3Ini"]
-    PhiIni = kwargs["PhiIni"]
+    Xstar = pde_parms["Xstar"]
+    Tstar = pde_parms["Tstar"]
+    max_depth = pde_parms["max_depth"]
+    ShallowLimit = pde_parms["ShallowLimit"]
+    DeepLimit = pde_parms["DeepLimit"]
+    CAIni = pde_parms["CAIni"]
+    CCIni = pde_parms["CCIni"]
+    cCaIni = pde_parms["cCaIni"]
+    cCO3Ini = pde_parms["cCO3Ini"]
+    PhiIni = pde_parms["PhiIni"]
 
-    Number_of_depths = kwargs["N"]
-    # End_time is in units of Tstar.
-    End_time = kwargs["tmax"]/Tstar
-    dt = kwargs["dt"]
+    Number_of_depths = pde_parms["N"]
 
     depths = CartesianGrid([[0, max_depth/Xstar]], [Number_of_depths], periodic=False)
     # We will be needing forward and backward differencing for
@@ -61,8 +58,8 @@ def integrate_equations(**kwargs):
     # Not all keys from kwargs are LMAHeureuxPorosityDiff arguments.
     # Taken from https://stackoverflow.com/questions/334655/passing-a-\
     # dictionary-to-a-function-as-keyword-parameters
-    filtered_kwargs = {k: v for k, v in kwargs.items() if k in [p.name for p in 
-                      inspect.signature(LMAHeureuxPorosityDiff).parameters.\
+    filtered_kwargs = {k: v for k, v in pde_parms.items() if k in [p.name for p
+                       in inspect.signature(LMAHeureuxPorosityDiff).parameters.\
                         values()]}
 
     eq = LMAHeureuxPorosityDiff(AragoniteSurface, CalciteSurface, CaSurface, 
@@ -77,36 +74,36 @@ def integrate_equations(**kwargs):
                    datetime.now().strftime("%d_%m_%Y_%H_%M_%S" + "/")
     os.makedirs(store_folder)
     stored_results = store_folder + "LMAHeureuxPorosityDiff.hdf5"
-    storage = FileStorage(stored_results, info=kwargs)
+    # Keep a record of all parameters.
+    storage_parms = solver_parms | tracker_parms | pde_parms
+    storage = FileStorage(stored_results, info=storage_parms)
 
-    if kwargs["live_plotting"]:
-        live_plots = LivePlotTracker(interval=kwargs["plotting_interval"], \
+    if tracker_parms["live_plotting"]:
+        live_plots = LivePlotTracker(interval=\
+                                     tracker_parms["plotting_interval"], \
                                      title="Integration results",
                                      show=True, max_fps=1, \
                                      plot_args ={"ax_style": {"ylim": (0, 1.5)}})
     else:
         live_plots = None
 
-    if kwargs["track_U_at_bottom"]:
+    if tracker_parms["track_U_at_bottom"]:
         data_tracker = DataTracker(eq.track_U_at_bottom, \
-                               interval = kwargs["data_tracker_interval"])
+                               interval = tracker_parms["data_tracker_interval"])
     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"])
+
+    sol, info = eq.solve(state, **solver_parms, tracker=["progress", \
+                         storage.tracker(\
+                             tracker_parms["progress_tracker_interval"]),\
+                         live_plots, data_tracker])
 
     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"]
 
-    if kwargs["track_U_at_bottom"]:
+    if tracker_parms["track_U_at_bottom"]:
         with h5py.File(stored_results, 'a') as hf:
             U_grp = hf.create_group("U")
             U_grp.create_dataset("U_at_bottom", \
@@ -135,10 +132,9 @@ def Plot_results(sol, covered_time, depths, Xstar, store_folder):
     fig.savefig(store_folder + 'Final_distributions.pdf', bbox_inches="tight")
 
 if __name__ == '__main__':
-    # Concatenate the dict containing the Scenario parameters with the
-    # dict containing the solver parameters (such as required tolerance) and
-    # with the dict containing the tracker parameters.
-    all_kwargs = asdict(Map_Scenario()) | asdict(Solver()) | asdict(Tracker())
+    pde_parms = asdict(Map_Scenario()) 
+    solver_parms = asdict(Solver()) 
+    tracker_parms = asdict(Tracker())
     solution, covered_time, depths, Xstar, store_folder = \
-        integrate_equations(**all_kwargs)
+        integrate_equations(solver_parms, tracker_parms, pde_parms)
     Plot_results(solution, covered_time, depths, Xstar, store_folder)

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]