diff --git a/Pipfile b/Pipfile index 27905f1..3469c3d 100644 --- a/Pipfile +++ b/Pipfile @@ -15,8 +15,6 @@ py-pde = "*" tqdm = "*" h5py = "*" pint = "*" -ipython = "*" -pandas = "*" [dev-packages] pylint = "*" diff --git a/README.md b/README.md index c0ac52f..3597600 100644 --- a/README.md +++ b/README.md @@ -1,17 +1,17 @@ # Integrating diagenetic equations using Python -This repo was created as an attempt to reproduce the plots shown at the kickoff of the AstroTOM ("Turing or Milankovitch") project by Niklas Hohmann, from his Matlab scripts (available at [github.com/MindTheGap-ERC/LMA-Matlab](https://github.com/MindTheGap-ERC/LMA-Matlab)). +This repo was created in an attempt to reproduce the plots shown at the kickoff of the AstroTOM ("Turing or Milankovitch") project by Niklas Hohmann, from his Matlab scripts (available at [github.com/MindTheGap-ERC/LMA-Matlab](https://github.com/MindTheGap-ERC/LMA-Matlab)). AstroTOM is an OpenSSI 2021b project from the Netherlands eScience Center and Utrecht University (UU). Dr. Emilia Jarochowska (UU) is the lead applicant of this project. After replacing central differencing for the gradients in the five diagenetic equations 40-43 from [L'Heureux (2018)](https://www.hindawi.com/journals/geofluids/2018/4968315/) by forward and backward differencing depending on the sign of U and W as a first step and a Fiadeiro-Veronis spatial difference scheme as a second step, it turns out that these equations can be integrated for more than 13.190 years (the full T*) with an implicit or explicit (in time) solver, but not with a simple Eulerian scheme. A Runge-Kutta solver, with an adaptive timestep will, however, suffice. -After correcting the value of b (5-->5e-4) it turned out that a stable integration is also possible without a Fiadeiro-Veronis scheme. The `main` branch makes use of a constant porosity diffusion coefficient. +After correcting the value of b (5-->5e-4) it turned out that a stable integration is also possible without a Fiadeiro-Veronis scheme. We currently make use of a constant porosity diffusion coefficient. -Implicit (in time) solvers with use of Jacobians (in functional forms, so without numerical approximations) are available in the `Use_solve_ivp_without_py-pde_wrapper` branch. +The implicit (in time) solvers that `solve_ivp` offers can be deployed with its numerically approximated Jacobians and a Jacobian sparsity matrix. -Wide use is made of the [py-pde](https://py-pde.readthedocs.io/en/latest/) package, especially in the `main` branch. +Wide use is made of the [py-pde](https://py-pde.readthedocs.io/en/latest/) package, especially `CartesianGrid` and `ScalarField`. ## Installing and using To run this code, you need `git` and `conda` or `pip` to install . diff --git a/marlpde/Evolve_scenario.py b/marlpde/Evolve_scenario.py index a8ac670..d576246 100644 --- a/marlpde/Evolve_scenario.py +++ b/marlpde/Evolve_scenario.py @@ -1,16 +1,20 @@ #!/usr/bin/env python -from datetime import datetime -import os from dataclasses import asdict import inspect -import matplotlib.pyplot as plt +import os +import time +from datetime import datetime import h5py -from pde import CartesianGrid, ScalarField, FileStorage, LivePlotTracker -from pde import DataTracker -from pde.grids.operators.cartesian import _make_derivative -from parameters import Map_Scenario, Solver, Tracker from LHeureux_model import LMAHeureuxPorosityDiff +from parameters import Map_Scenario, Solver, Tracker +from pde import CartesianGrid, ScalarField +from pde.grids.operators.cartesian import _make_derivative +from scipy.integrate import solve_ivp +from tqdm import tqdm +import matplotlib +import matplotlib.pyplot as plt +matplotlib.use("AGG") def integrate_equations(solver_parms, tracker_parms, pde_parms): ''' @@ -41,12 +45,6 @@ def integrate_equations(solver_parms, tracker_parms, pde_parms): depths.register_operator("grad_forw", \ lambda grid: _make_derivative(grid, method="forward")) - AragoniteSurface = ScalarField(depths, CAIni) - CalciteSurface = ScalarField(depths, CCIni) - CaSurface = ScalarField(depths, cCaIni) - CO3Surface = ScalarField(depths, cCO3Ini) - PorSurface = ScalarField(depths, PhiIni) - # I need those two fields for computing coA, which is rather involved. # There may be a simpler way of selecting these depths, but I haven't # found one yet. For now these two Heaviside step functions. @@ -55,77 +53,152 @@ def integrate_equations(solver_parms, tracker_parms, pde_parms): not_too_deep = ScalarField.from_expression(depths, f"heaviside({DeepLimit/Xstar}-x, 0)") - # Not all keys from kwargs are LMAHeureuxPorosityDiff arguments. + # Not all keys from pde_parms 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 pde_parms.items() if k in [p.name for p - in inspect.signature(LMAHeureuxPorosityDiff).parameters.\ - values()]} + filtered_pde_parms = {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, - CO3Surface, PorSurface, not_too_shallow, - not_too_deep, **filtered_kwargs) + slices_all_fields = [slice(i * Number_of_depths, (i+1) * Number_of_depths) \ + for i in range(5)] + + eq = LMAHeureuxPorosityDiff(depths, slices_all_fields, not_too_shallow, + not_too_deep, **filtered_pde_parms) + + # Setting initial values for all five fields in three steps. + # This is perhaps more involved than necessary, since the + # end result is a 1D Numpy array with a length of five times + # the number of depths. The number five comes from the five + # fields we are tracking, by integrating five pdes. + # First step. + AragoniteSurface = ScalarField(depths, CAIni) + CalciteSurface = ScalarField(depths, CCIni) + CaSurface = ScalarField(depths, cCaIni) + CO3Surface = ScalarField(depths, cCO3Ini) + PorSurface = ScalarField(depths, PhiIni) + # Second step. state = eq.get_state(AragoniteSurface, CalciteSurface, CaSurface, CO3Surface, PorSurface) + # Third step. + y0 = state.data.ravel() + no_progress_updates = tracker_parms["no_progress_updates"] + + start_computing = time.time() + with tqdm(total=no_progress_updates) as pbar: + t0 = solver_parms["t_span"][0] + end_time = solver_parms["t_span"][1] + progress_bar_args = [pbar, (end_time - t0) / no_progress_updates, t0] + + # The backend parameter determines which function should be called to + # determine the right-hand sides of the five pdes. It can be either + # a Numpy-based function or a Numba-based function. The latter is + # faster. + backend = solver_parms["backend"] + # "backend" is not an argument for solve_ivp, so remove it now. + del solver_parms["backend"] + + sol = solve_ivp(eq.fun if backend=="numpy" else eq.fun_numba, + y0=y0, **solver_parms, + t_eval= tracker_parms["t_eval"], + events = [eq.zeros, eq.zeros_CA, eq.zeros_CC, \ + eq.ones_CA_plus_CC, eq.ones_Phi, eq.zeros_U, \ + eq.zeros_W], args=progress_bar_args) + end_computing = time.time() + + print() + print(f"Number of rhs evaluations = {sol.nfev} \n") + print(f"Number of Jacobian evaluations = {sol.njev} \n") + print(f"Number of LU decompositions = {sol.nlu} \n") + print(f"Status = {sol.status} \n") + print(f"Success = {sol.success} \n") + f = sol.t_events[0] + print(("Times, in years, at which any field at any depth crossed zero: "\ + +', '.join(['%.2f']*len(f))+"") % tuple([Tstar * time for time in f])) + print() + g = sol.t_events[1] + print(("Times, in years, at which CA at any depth crossed zero: "\ + +', '.join(['%.2f']*len(g))+"") % tuple([Tstar * time for time in g])) + print() + h = sol.t_events[2] + print(("Times, in years, at which CC at any depth crossed zero: "\ + +', '.join(['%.2f']*len(h))+"") % tuple([Tstar * time for time in h])) + print() + k = sol.t_events[3] + print(("Times, in years, at which CA + CC at any depth crossed one: "\ + +', '.join(['%.2f']*len(k))+"") % tuple([Tstar * time for time in k])) + print() + l = sol.t_events[4] + print(("Times, in years, at which the porosity at any depth crossed one: "\ + +', '.join(['%.2f']*len(l))+"") % tuple([Tstar * time for time in l])) + print() + m = sol.t_events[5] + print(("Times, in years, at which U at any depth crossed zero: "\ + +', '.join(['%.2f']*len(m))+"") % tuple([Tstar * time for time in m])) + print() + n = sol.t_events[6] + print(("Times, in years, at which W at any depth crossed zero: "\ + +', '.join(['%.2f']*len(n))+"") % tuple([Tstar * time for time in n])) + print() + + print(f"Message from solve_ivp = {sol.message} \n") + print(("Time taken for solve_ivp is " + f"{end_computing - start_computing:.2e}s. \n")) + + if sol.status == 0: + covered_time = Tstar * end_time + else: + covered_time = pbar.n * Tstar * end_time / no_progress_updates + # Store your results somewhere in a subdirectory of a parent directory. store_folder = "../Results/" + \ datetime.now().strftime("%d_%m_%Y_%H_%M_%S" + "/") os.makedirs(store_folder) stored_results = store_folder + "LMAHeureuxPorosityDiff.hdf5" # Keep a record of all parameters. - storage_parms = solver_parms | tracker_parms | pde_parms - storage = FileStorage(stored_results, info=storage_parms) - - 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 + stored_parms = solver_parms | tracker_parms | pde_parms + # Remove items which will raise a problem when storing as metadata in an + # hdf5 file + if "jac_sparsity" in stored_parms: + del stored_parms["jac_sparsity"] - if tracker_parms["track_U_at_bottom"]: - data_tracker = DataTracker(eq.track_U_at_bottom, \ - interval = tracker_parms["data_tracker_interval"]) - else: - data_tracker = None + # Reshape the solutions such that you get one row per field. + # The third axis, i.e.the time axis, can remain unchanged. + field_solutions = sol.y.reshape(5, Number_of_depths, sol.y.shape[-1]) - sol, info = eq.solve(state, **solver_parms, tracker=["progress", \ - storage.tracker(\ - tracker_parms["progress_tracker_interval"]),\ - live_plots, data_tracker]) + with h5py.File(stored_results, "w") as stored: + stored.create_dataset("solutions", data=field_solutions) + stored.create_dataset("times", data=sol.t) + for event_index, _ in enumerate(sol.t_events): + stored.create_dataset("event_" + str(event_index), + data=sol.t_events[event_index]) + stored.attrs.update(stored_parms) - print() - print(f"Meta-information about the solution : {info} \n") - - covered_time_span = Tstar * info["controller"]["t_final"] - - 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", \ - data=data_tracker.dataframe.to_numpy()) + # We will be plotting only the distributions corresponding to the last time. + # Thus no point in returning all data. Moreover, all data have been saved. - return sol, covered_time_span, depths, Xstar, store_folder + return field_solutions[:, :, -1], covered_time, depths, Xstar, store_folder -def Plot_results(sol, covered_time, depths, Xstar, store_folder): +def Plot_results(last_field_sol, covered_time, depths, Xstar, store_folder): ''' Plot the five fields at the end of the integration interval as a function of depth. ''' fig, ax = plt.subplots() - fig.suptitle(f"Distributions after {covered_time:.2f} years") + fig.suptitle(f"Distributions after {covered_time:.2e} years") # Marker size ms = 5 plotting_depths = ScalarField.from_expression(depths, "x").data * Xstar - ax.plot(plotting_depths, sol.data[0], "v", ms = ms, label = "CA") - ax.plot(plotting_depths, sol.data[1], "^", ms = ms, label = "CC") - ax.plot(plotting_depths, sol.data[2], ">", ms = ms, label = "cCa") - ax.plot(plotting_depths, sol.data[3], "<", ms = ms, label = "cCO3") - ax.plot(plotting_depths, sol.data[4], "o", ms = ms, label = "Phi") + + ax.plot(plotting_depths, last_field_sol[0], "v", ms = ms, label = "CA") + ax.plot(plotting_depths, last_field_sol[1], "^", ms = ms, label = "CC") + ax.plot(plotting_depths, last_field_sol[2], ">", ms = ms, label = "cCa") + ax.plot(plotting_depths, last_field_sol[3], "<", ms = ms, label = "cCO3") + ax.plot(plotting_depths, last_field_sol[4], "o", ms = ms, label = "Phi") + ax.set_xlabel("Depth (cm)") ax.set_ylabel("Compositions and concentrations (dimensionless)") ax.legend(loc='upper right') @@ -135,6 +208,4 @@ def Plot_results(sol, covered_time, depths, Xstar, store_folder): pde_parms = asdict(Map_Scenario()) solver_parms = asdict(Solver()) tracker_parms = asdict(Tracker()) - solution, covered_time, depths, Xstar, store_folder = \ - integrate_equations(solver_parms, tracker_parms, pde_parms) - Plot_results(solution, covered_time, depths, Xstar, store_folder) + Plot_results(*integrate_equations(solver_parms, tracker_parms, pde_parms)) \ No newline at end of file diff --git a/marlpde/LHeureux_model.py b/marlpde/LHeureux_model.py index b5422ac..80b6f4d 100644 --- a/marlpde/LHeureux_model.py +++ b/marlpde/LHeureux_model.py @@ -1,26 +1,33 @@ +from pde.grids.operators.cartesian import _make_derivative +from pde import FieldCollection, ScalarField import numpy as np -from pde import FieldCollection, PDEBase, ScalarField from numba import njit -np.seterr(divide="raise", over="raise", under="raise", invalid="raise") - -class LMAHeureuxPorosityDiff(PDEBase): - - def __init__(self, AragoniteSurface, CalciteSurface, CaSurface, - CO3Surface, PorSurface, not_too_shallow, not_too_deep, CA0, CC0, - cCa0, cCO30, Phi0, sedimentationrate, Xstar, Tstar, k1, k2, k3, - k4, m1, m2, n1, n2, b, beta, rhos, rhow, rhos0, KA, KC, muA, - D0Ca, PhiNR, PhiInfty, PhiIni, DCa, DCO3): - - self.AragoniteSurface = AragoniteSurface - self.CalciteSurface = CalciteSurface - self.CaSurface = CaSurface - self.CO3Surface = CO3Surface - self.PorSurface = PorSurface - self.bc_CA = [{"value": CA0}, {"curvature": 0}] - self.bc_CC = [{"value": CC0}, {"curvature": 0}] + +np.seterr(divide="raise", over="raise", under="warn", invalid="raise") + +class LMAHeureuxPorosityDiff(): + ''' + This class defines all the parameters for the diagenetic model. + ''' + def __init__(self, Depths, slices_all_fields, not_too_shallow, not_too_deep, + CA0, CC0, cCa0, cCO30, Phi0, sedimentationrate, Xstar, Tstar, + k1, k2, k3, k4, m1, m2, n1, n2, b, beta, rhos, rhow, rhos0, + KA, KC, muA, D0Ca, PhiNR, PhiInfty, PhiIni, DCa, DCO3, + FV_switch): + self.no_fields = 5 + self.Depths = Depths + self.Depths.register_operator("grad_back", \ + lambda grid: _make_derivative(grid, method="backward")) + self.Depths.register_operator("grad_forw", \ + lambda grid: _make_derivative(grid, method="forward")) + self.delta_x = self.Depths._axes_coords[0][1] - \ + self.Depths._axes_coords[0][0] + self.slices_all_fields = slices_all_fields + 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}] + self.bc_Phi = [{"value": Phi0}, {"derivative": 0}] self.sedimentationrate = sedimentationrate self.Xstar = Xstar self.Tstar = Tstar @@ -49,9 +56,6 @@ def __init__(self, AragoniteSurface, CalciteSurface, CaSurface, self.Phi0 = Phi0 self.DCa = DCa self.DCO3 = DCO3 - self.not_too_shallow = not_too_shallow - self.not_too_deep = not_too_deep - self.g = 100 * 9.81 self.dCa = self.DCa / self.D0Ca self.dCO3 = self.DCO3 / self.D0Ca @@ -65,27 +69,72 @@ def __init__(self, AragoniteSurface, CalciteSurface, CaSurface, self.rhorat = (self.rhos / self.rhow - 1) * self.beta / \ self.sedimentationrate self.presum = 1 - self.rhorat0 * self.Phi0 ** 3 * \ - (1 - np.exp(10 - 10 / self.Phi0)) / (1 - self.Phi0) - - # Fiadeiro-Veronis differentiation involves a coth and a reciprocal, which can + (1 - np.exp(10 - 10 / self.Phi0)) / (1 - self.Phi0) + self.not_too_shallow = not_too_shallow.data + self.not_too_deep = not_too_deep.data + self.FV_switch = FV_switch + + self.CA_sl= self.slices_all_fields[0] + self.CC_sl = self.slices_all_fields[1] + self.cCa_sl = self.slices_all_fields[2] + self.cCO3_sl = self.slices_all_fields[3] + self.Phi_sl = self.slices_all_fields[4] + + # Fiadeiro-Veronis integration involves a coth and a reciprocal, which can # easily lead to FloatingPointError: overflow encountered in double_scalars. # To avoid this, better revert to either backwards or central differencing # the Peclet number is very large or very small. self.Peclet_min = 1e-2 self.Peclet_max = 1/self.Peclet_min - # Need this number for Fiadeiro-Veronis differentiation. - self.delta_x = self.AragoniteSurface.grid._axes_coords[0][1] - \ - self.AragoniteSurface.grid._axes_coords[0][0] + + self.last_t = 0. # Helper variable for progress bar. + + # Make operators for use in Numba compiled (i.e. jit compiled) functions. + # Instead of the default central differenced gradient from py-pde + # construct forward and backward differenced gradients and apply + # either one of them, based on the sign of U. + self.CA_grad_back_op = self.Depths.make_operator("grad_back", self.bc_CA) + self.CA_grad_forw_op = self.Depths.make_operator("grad_forw", self.bc_CA) + self.CC_grad_back_op = self.Depths.make_operator("grad_back", self.bc_CC) + self.CC_grad_forw_op = self.Depths.make_operator("grad_forw", self.bc_CC) + # Instead of the default central differenced gradient from py-pde + # construct forward and backward differenced gradients and give them + # appropriate weights according to a Fiadeiro-Veronis scheme. + self.cCa_grad_back_op = self.Depths.make_operator("grad_back", self.bc_cCa) + self.cCa_grad_forw_op = self.Depths.make_operator("grad_forw", self.bc_cCa) + self.cCa_laplace_op = self.Depths.make_operator("laplace", self.bc_cCa) + self.cCO3_grad_back_op = self.Depths.make_operator("grad_back", self.bc_cCO3) + self.cCO3_grad_forw_op = self.Depths.make_operator("grad_forw", self.bc_cCO3) + self.cCO3_laplace_op = self.Depths.make_operator("laplace", self.bc_cCO3) + self.Phi_grad_back_op = self.Depths.make_operator("grad_back", self.bc_Phi) + self.Phi_grad_forw_op = self.Depths.make_operator("grad_forw", self.bc_Phi) + self.Phi_laplace_op = self.Depths.make_operator("laplace", self.bc_Phi) + + # Make sure integration stops when we field values become less than zero + # or more than one, in some cases. Or just register that this is happening + # and continue integration, which corresponds to "False". + setattr(self.zeros.__func__, "terminal", False) + setattr(self.zeros_CA.__func__, "terminal", False) + setattr(self.zeros_CC.__func__, "terminal", False) + setattr(self.ones_CA_plus_CC.__func__, "terminal", False) + setattr(self.ones_Phi.__func__, "terminal", False) + setattr(self.zeros_U.__func__, "terminal", False) + setattr(self.zeros_W.__func__, "terminal", False) + + # Accommodate for fixed porosity diffusion coefficient. + # Beware that the functional Jacobian has not been aligned with this, + # i.e. it was derived using a porosity diffusion coefficient depending + # on a time-varying porosity. This will all be fine if you do not set + # 'jac=eq.jac' in solve_ivp, but leave it empty such that a numerically + # approximated Jacobian will be used. 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, PorSurface): - # Return initial state and register forward and backward difference - # operators. + """generate a suitable initial state""" AragoniteSurface.label = "ARA" CalciteSurface.label = "CAL" CaSurface.label = "Ca" @@ -95,39 +144,10 @@ def get_state(self, AragoniteSurface, CalciteSurface, CaSurface, CO3Surface, return FieldCollection([AragoniteSurface, CalciteSurface, CaSurface, CO3Surface, PorSurface]) - - def track_U_at_bottom(self, state, t): - # First, extract the porosity at the bottom of the system. - # The derived quantities will then also be at the bottom of the system. - Phi = state.data[4][-1] - F = 1 - np.exp(10 - 10 / Phi) - one_minus_Phi = 1 - Phi - U_bottom = self.presum + self.rhorat * Phi ** 3 * F /one_minus_Phi - return {"U at bottom": U_bottom} - @staticmethod @njit def calculate_sigma(Peclet, W_data, Peclet_min, Peclet_max): - ''' Calculate sigma following formula 8.73 from Boudreau: - "Diagenetic Models and their implementation" - - Parameters - ---------- - Peclet: ndarray(dtype=float, ndim=1) - Array along depth of the Peclet numbers. - W-data: ndarray(dtype=float, ndim=1) - Array along depth of the velocity of the pore water - (counted positive downwards) - Peclet_min: float - Lower limit for calculating 1/tanh(Peclet) - Peclet_max: float - Upper limit for calculating 1/tanh(Peclet) - - Returns - ------- - sigma: ndarray(dtype=float, ndim=1) - Array along depth with sigma values - ''' + # Assuming the arrays are 1D. sigma = np.empty(Peclet.size) for i in range(sigma.size): if np.abs(Peclet[i]) < Peclet_min: @@ -137,10 +157,25 @@ def calculate_sigma(Peclet, W_data, Peclet_min, Peclet_max): else: sigma[i] = np.cosh(Peclet[i])/np.sinh(Peclet[i]) - \ 1/Peclet[i] - return sigma - - def evolution_rate(self, state, t=0): - CA, CC, cCa, cCO3, Phi = state + return sigma + + def fun(self, t, y, progress_proxy, progress_dt, t0): + """ + For tqdm to monitor progress. + From + https://stackoverflow.com/questions/59047892/ + how-to-monitor-the-process-of-scipy-odeint """ + if self.last_t == 0.: + self.last_t = t0 + n = int((t - self.last_t) / progress_dt) + progress_proxy.update(n) + self.last_t += n * progress_dt + + CA = ScalarField(self.Depths, y[self.slices_all_fields[0]]) + CC = ScalarField(self.Depths, y[self.slices_all_fields[1]]) + cCa = ScalarField(self.Depths, y[self.slices_all_fields[2]]) + cCO3 = ScalarField(self.Depths, y[self.slices_all_fields[3]]) + Phi = ScalarField(self.Depths, y[self.slices_all_fields[4]]) two_factors = cCa * cCO3 two_factors_upp_lim = two_factors.to_scalar(lambda f: np.fmin(f,1)) @@ -158,20 +193,20 @@ def evolution_rate(self, state, t=0): (1 - two_factors_upp_lim) ** self.n2) F = 1 - np.exp(10 - 10 / Phi) - one_minus_Phi = 1 - Phi - U = self.presum + self.rhorat * Phi ** 3 * F /one_minus_Phi - - # Choose either forward or backward differencing for CA and CC - # depending on the sign of U + U = self.presum + self.rhorat * Phi ** 3 * F / (1 - Phi) + + # Instead of the default central differenced gradient from py-pde + # construct forward and backward differenced gradients and apply + # either one of them, based on the sign of U. CA_grad_back = CA.apply_operator("grad_back", self.bc_CA) 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 = ScalarField(self.Depths, 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_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 = CC.apply_operator("grad_forw", self.bc_CC) + CC_grad = ScalarField(self.Depths, np.where(U.data>0, CC_grad_back.data, \ CC_grad_forw.data)) W = self.presum - self.rhorat * Phi ** 2 * F @@ -184,182 +219,218 @@ def evolution_rate(self, state, t=0): # Implementing equation 6 from l'Heureux. denominator = 1 - 2 * np.log(Phi) - - # Fiadeiro-Veronis scheme for equations 42 and 43 - # from l'Heureux. - common_Peclet = W * self.delta_x / 2. - Peclet_cCa = common_Peclet * denominator/ self.dCa - sigma_cCa_data = LMAHeureuxPorosityDiff.calculate_sigma(\ - Peclet_cCa.data, W.data, self.Peclet_min, self.Peclet_max) - sigma_cCa = ScalarField(state.grid, sigma_cCa_data) - - Peclet_cCO3 = common_Peclet * denominator/ self.dCO3 - sigma_cCO3_data = LMAHeureuxPorosityDiff.calculate_sigma(\ - Peclet_cCO3.data, W.data, self.Peclet_min, self.Peclet_max) - sigma_cCO3 = ScalarField(state.grid, sigma_cCO3_data) - - # dPhi = self.auxcon * F * (Phi ** 3) / one_minus_Phi + # dPhi = self.auxcon * F * (Phi ** 3) / (1 - Phi) dPhi = self.dPhi_fixed - Peclet_Phi = common_Peclet / dPhi - sigma_Phi_data = LMAHeureuxPorosityDiff.calculate_sigma(\ - Peclet_Phi.data, W.data, self.Peclet_min, self.Peclet_max) - sigma_Phi = ScalarField(state.grid, sigma_Phi_data) - + if self.FV_switch: + # Fiadeiro-Veronis scheme for equations 42 and 43 + # from l'Heureux. + common_Peclet = W * self.delta_x / 2. + Peclet_cCa = common_Peclet * denominator/ self.dCa + sigma_cCa_data = LMAHeureuxPorosityDiff.calculate_sigma(\ + Peclet_cCa.data, W.data, self.Peclet_min, self.Peclet_max) + sigma_cCa = ScalarField(self.Depths, sigma_cCa_data) + + Peclet_cCO3 = common_Peclet * denominator/ self.dCO3 + sigma_cCO3_data = LMAHeureuxPorosityDiff.calculate_sigma(\ + Peclet_cCO3.data, W.data, self.Peclet_min, self.Peclet_max) + sigma_cCO3 = ScalarField(self.Depths, sigma_cCO3_data) + + Peclet_Phi = common_Peclet / dPhi + sigma_Phi_data = LMAHeureuxPorosityDiff.calculate_sigma(\ + Peclet_Phi.data, W.data, self.Peclet_min, self.Peclet_max) + sigma_Phi = ScalarField(self.Depths, sigma_Phi_data) + else: + sigma_cCa = 0 + sigma_cCO3 = 0 + sigma_Phi = 0 + + # Instead of the default central differenced gradient from py-pde + # construct forward and backward differenced gradients and give them + # appropriate weights according to a Fiadeiro-Veronis scheme. cCa_grad_back = cCa.apply_operator("grad_back", self.bc_cCa) - cCa_grad_forw = cCa.apply_operator("grad_forw", self.bc_cCa) + cCa_grad_forw = cCa.apply_operator("grad_forw", self.bc_cCa) cCa_grad = 0.5 * ((1-sigma_cCa) * cCa_grad_forw +\ (1+sigma_cCa) * cCa_grad_back) + cCa_laplace = cCa.laplace(self.bc_cCa) cCO3_grad_back = cCO3.apply_operator("grad_back", self.bc_cCO3) cCO3_grad_forw = cCO3.apply_operator("grad_forw", self.bc_cCO3) cCO3_grad = 0.5 * ((1-sigma_cCO3) * cCO3_grad_forw +\ (1+sigma_cCO3) * cCO3_grad_back) + cCO3_laplace = cCO3.laplace(self.bc_cCO3) Phi_grad_back = Phi.apply_operator("grad_back", self.bc_Phi) Phi_grad_forw = Phi.apply_operator("grad_forw", self.bc_Phi) Phi_grad = 0.5 * ((1-sigma_Phi) * Phi_grad_forw +\ (1+sigma_Phi) * Phi_grad_back) + Phi_laplace = Phi.laplace(self.bc_Phi) Phi_denom = Phi/denominator grad_Phi_denom = Phi_grad * (denominator + 2) / denominator ** 2 common_helper = coA - self.lambda_ * coC + + dcCa_dt = (cCa_grad * grad_Phi_denom + Phi_denom * cCa_laplace) \ + * self.dCa /Phi - W * cCa_grad \ + + self.Da * (1 - Phi) * (self.delta - cCa) * common_helper / Phi - dcCa_dt = (cCa_grad * grad_Phi_denom + Phi_denom * cCa.laplace(self.bc_cCa)) \ - * self.dCa /Phi -W * cCa_grad \ - + self.Da * one_minus_Phi * (self.delta - cCa) * common_helper / Phi - - dcCO3_dt = (cCO3_grad * grad_Phi_denom + Phi_denom * cCO3.laplace(self.bc_cCO3)) \ - * self.dCO3/Phi -W * cCO3_grad \ - + self.Da * one_minus_Phi * (self.delta - cCO3) * common_helper / Phi + dcCO3_dt = (cCO3_grad * grad_Phi_denom + Phi_denom * cCO3_laplace) \ + * self.dCO3/Phi - W * cCO3_grad \ + + self.Da * (1 - Phi) * (self.delta - cCO3) * common_helper / Phi dW_dx = -self.rhorat * Phi_grad * (2 * Phi * F + 10 * (F - 1)) - # This is closer to the original form of (43) from l' Heureux than - # the Matlab implementation. dPhi_dt = - (Phi * dW_dx + W * Phi_grad) \ - + dPhi * Phi.laplace(self.bc_Phi) \ - + self.Da * one_minus_Phi * common_helper - - return FieldCollection([dCA_dt, dCC_dt, dcCa_dt, dcCO3_dt, dPhi_dt]) - - def _make_pde_rhs_numba(self, state): - """ the numba-accelerated evolution equation """ - # make attributes locally available + + dPhi * Phi_laplace \ + + self.Da * (1 - Phi) * common_helper + + return FieldCollection([dCA_dt, dCC_dt, dcCa_dt, dcCO3_dt, dPhi_dt]).data.ravel() + + def fun_numba(self, t, y, progress_proxy, progress_dt, t0): + """ For tqdm to monitor progress. + From + https://stackoverflow.com/questions/59047892/ + how-to-monitor-the-process-of-scipy-odeint """ + if self.last_t == 0.: + self.last_t = t0 + n = int((t - self.last_t) / progress_dt) + progress_proxy.update(n) + self.last_t += n * progress_dt + + """ the numba-accelerated evolution equation """ KRat = self.KRat m1 = self.m1 - m2 = self.m2 + m2 = self.m2 n1 = self.n1 n2 = self.n2 nu1 = self.nu1 nu2 = self.nu2 - not_too_deep = self.not_too_deep.data - not_too_shallow = self.not_too_shallow.data + not_too_deep = self.not_too_deep + not_too_shallow = self.not_too_shallow presum = self.presum - rhorat= self.rhorat + rhorat = self.rhorat lambda_ = self.lambda_ - Da = self.Da + Da= self.Da dCa = self.dCa dCO3 = self.dCO3 delta = self.delta - # The following three numbers are also needed for Fiadeiro-Veronis. - dPhi_fixed = self.dPhi_fixed + auxcon = self.auxcon + delta_x = self.delta_x Peclet_min = self.Peclet_min Peclet_max = self.Peclet_max - delta_x = state.grid._axes_coords[0][1] - state.grid._axes_coords[0][0] - 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) - grad_back_cCO3 = state.grid.make_operator("grad_back", bc = self.bc_cCO3) - grad_forw_cCO3 = state.grid.make_operator("grad_forw", bc = self.bc_cCO3) - laplace_cCO3 = state.grid.make_operator("laplace", bc = self.bc_cCO3) - grad_back_Phi = state.grid.make_operator("grad_back", bc = self.bc_Phi) - grad_forw_Phi = state.grid.make_operator("grad_forw", bc = self.bc_Phi) - laplace_Phi = state.grid.make_operator("laplace", bc = self.bc_Phi) - - @njit - def pde_rhs(state_data, t=0): - """ compiled helper function evaluating right hand side """ - # Instead of the default central differenced gradient from py-pde - # construct forward and backward differenced gradients and apply - # 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) - cCa_grad_forw = grad_forw_cCa(cCa) - cCa_laplace = laplace_cCa(cCa) - - cCO3 = state_data[3] - cCO3_grad_back = grad_back_cCO3(cCO3) - cCO3_grad_forw = grad_forw_cCO3(cCO3) - cCO3_laplace = laplace_cCO3(cCO3) - - Phi = state_data[4] - Phi_grad_back = grad_back_Phi(Phi) - Phi_grad_forw = grad_forw_Phi(Phi) - Phi_laplace = laplace_Phi(Phi) - - rate = np.empty_like(state_data) - - # state_data.size should be the same as len(CA) or len(CC), check this. - # So the number of depths, really. - no_depths = state_data[0].size - - denominator = np.empty(no_depths) - common_helper1 = np.empty(no_depths) - common_helper2 = np.empty(no_depths) - helper_cCa_grad = np.empty(no_depths) - helper_cCO3_grad = np.empty(no_depths) - F = np.empty(no_depths) - U = np.empty(no_depths) - W = np.empty(no_depths) - two_factors = np.empty(no_depths) - two_factors_upp_lim = np.empty(no_depths) - two_factors_low_lim = np.empty(no_depths) - three_factors = np.empty(no_depths) - three_factors_upp_lim = np.empty(no_depths) - three_factors_low_lim = np.empty(no_depths) - coA = np.empty(no_depths) - coC = np.empty(no_depths) - common_helper3 = np.empty(no_depths) - dPhi = np.empty(no_depths) - dW_dx = np.empty(no_depths) - one_minus_Phi = np.empty(no_depths) - CA_grad = np.empty(no_depths) - CC_grad = np.empty(no_depths) - cCa_grad = np.empty(no_depths) - cCO3_grad = np.empty(no_depths) - Phi_grad = np.empty(no_depths) - - for i in range(no_depths): - F[i] = 1 - np.exp(10 - 10 / Phi[i]) - - U[i] = presum + rhorat * Phi[i] ** 3 * F[i]/ (1 - Phi[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] + no_depths = self.Depths.shape[0] + dPhi_fixed = self.dPhi_fixed + FV_switch = self.FV_switch + + CA = y[self.CA_sl] + CC = y[self.CC_sl] + cCa = y[self.cCa_sl] + cCO3 = y[self.cCO3_sl] + Phi = y[self.Phi_sl] + + CA_grad_back_op = self.CA_grad_back_op + CA_grad_forw_op = self.CA_grad_forw_op + CC_grad_back_op = self.CC_grad_back_op + CC_grad_forw_op = self.CC_grad_forw_op + cCa_grad_back_op = self.cCa_grad_back_op + cCa_grad_forw_op = self.cCa_grad_forw_op + cCa_laplace_op = self.cCa_laplace_op + cCO3_grad_back_op = self.cCO3_grad_back_op + cCO3_grad_forw_op = self.cCO3_grad_forw_op + cCO3_laplace_op = self.cCO3_laplace_op + Phi_grad_back_op = self.Phi_grad_back_op + Phi_grad_forw_op = self.Phi_grad_forw_op + Phi_laplace_op = self.Phi_laplace_op + + return LMAHeureuxPorosityDiff.pde_rhs(CA, CC, cCa, cCO3, Phi, + CA_grad_back_op, CA_grad_forw_op, + CC_grad_back_op, CC_grad_forw_op, + cCa_grad_back_op, cCa_grad_forw_op, + cCa_laplace_op, cCO3_grad_back_op, + cCO3_grad_forw_op, cCO3_laplace_op, + Phi_grad_back_op, Phi_grad_forw_op, + Phi_laplace_op, KRat, m1, m2, n1, + n2, nu1, nu2, not_too_deep, + not_too_shallow, presum, rhorat, + lambda_, Da, dCa, dCO3, delta, + auxcon, delta_x, Peclet_min, + Peclet_max, no_depths, dPhi_fixed, + FV_switch) - W[i] = presum - rhorat * Phi[i] ** 2 * F[i] + @staticmethod + @njit(cache=False) + def pde_rhs(CA, CC, cCa, cCO3, Phi, CA_grad_back_op, CA_grad_forw_op, + CC_grad_back_op, CC_grad_forw_op, + cCa_grad_back_op, cCa_grad_forw_op, cCa_laplace_op, + cCO3_grad_back_op, cCO3_grad_forw_op, cCO3_laplace_op, + Phi_grad_back_op, Phi_grad_forw_op, Phi_laplace_op, + KRat, m1, m2, n1, n2, nu1, nu2, not_too_deep, not_too_shallow, + presum, rhorat, lambda_, Da, dCa, dCO3, delta, auxcon, delta_x, + Peclet_min, Peclet_max, no_depths, dPhi_fixed, FV_switch): + + CA_grad_back = CA_grad_back_op(CA) + CA_grad_forw = CA_grad_forw_op(CA) + CC_grad_back = CC_grad_back_op(CC) + CC_grad_forw = CC_grad_forw_op(CC) + cCa_grad_back = cCa_grad_back_op(cCa) + cCa_grad_forw = cCa_grad_forw_op(cCa) + cCa_laplace = cCa_laplace_op(cCa) + cCO3_grad_back = cCO3_grad_back_op(cCO3) + cCO3_grad_forw = cCO3_grad_forw_op(cCO3) + cCO3_laplace = cCO3_laplace_op(cCO3) + Phi_grad_back = Phi_grad_back_op(Phi) + Phi_grad_forw = Phi_grad_forw_op(Phi) + Phi_laplace = Phi_laplace_op(Phi) + + denominator = np.empty(no_depths) + common_helper1 = np.empty(no_depths) + common_helper2 = np.empty(no_depths) + helper_cCa_grad = np.empty(no_depths) + helper_cCO3_grad = np.empty(no_depths) + F = np.empty(no_depths) + U = np.empty(no_depths) + W = np.empty(no_depths) + two_factors = np.empty(no_depths) + two_factors_upp_lim = np.empty(no_depths) + two_factors_low_lim = np.empty(no_depths) + three_factors = np.empty(no_depths) + three_factors_upp_lim = np.empty(no_depths) + three_factors_low_lim = np.empty(no_depths) + coA = np.empty(no_depths) + coC = np.empty(no_depths) + common_helper3 = np.empty(no_depths) + rate = np.empty(5 * no_depths) + dPhi = np.empty(no_depths) + dW_dx = np.empty(no_depths) + one_minus_Phi = np.empty(no_depths) + CA_grad = np.empty(no_depths) + CC_grad = np.empty(no_depths) + cCa_grad = np.empty(no_depths) + cCO3_grad = np.empty(no_depths) + Phi_grad = np.empty(no_depths) + + for i in range(no_depths): + F[i] = 1 - np.exp(10 - 10 / Phi[i]) + + U[i] = presum + rhorat * Phi[i] ** 3 * F[i]/ (1 - Phi[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] - # Implementing equation 6 from l'Heureux. - denominator[i] = 1 - 2 * np.log(Phi[i]) + # Implementing equation 6 from l'Heureux. + denominator[i] = 1 - 2 * np.log(Phi[i]) + one_minus_Phi[i] = 1 - Phi[i] + # dPhi[i] = auxcon * F[i] * (Phi[i] ** 3) / one_minus_Phi[i] + dPhi[i] = dPhi_fixed + if FV_switch: # Fiadeiro-Veronis scheme for equations 42 and 43 # from l'Heureux. Peclet_cCa = W[i] * delta_x * denominator[i]/ (2. * dCa ) @@ -368,8 +439,7 @@ def pde_rhs(state_data, t=0): elif np.abs(Peclet_cCa) > Peclet_max: sigma_cCa = np.sign(W[i]) else: - sigma_cCa = np.cosh(Peclet_cCa)/np.sinh(Peclet_cCa) - \ - 1/Peclet_cCa + sigma_cCa = np.cosh(Peclet_cCa)/np.sinh(Peclet_cCa) - 1/Peclet_cCa Peclet_cCO3 = W[i] * delta_x * denominator[i]/ (2. * dCO3) if np.abs(Peclet_cCO3) < Peclet_min: @@ -377,78 +447,148 @@ def pde_rhs(state_data, t=0): elif np.abs(Peclet_cCO3) > Peclet_max: sigma_cCO3 = np.sign(W[i]) else: - sigma_cCO3 = np.cosh(Peclet_cCO3)/np.sinh(Peclet_cCO3) - \ - 1/Peclet_cCO3 + sigma_cCO3 = np.cosh(Peclet_cCO3)/np.sinh(Peclet_cCO3) - 1/Peclet_cCO3 - one_minus_Phi[i] = 1 - Phi[i] - dPhi[i] = dPhi_fixed Peclet_Phi = W[i] * delta_x / (2. * dPhi[i]) if np.abs(Peclet_Phi) < Peclet_min: sigma_Phi = 0 elif np.abs(Peclet_Phi) > Peclet_max: sigma_Phi = np.sign(W[i]) else: - sigma_Phi = np.cosh(Peclet_Phi)/np.sinh(Peclet_Phi) - \ - 1/Peclet_Phi - - cCa_grad[i] = 0.5 * ((1-sigma_cCa) * cCa_grad_forw[i] + \ - (1+sigma_cCa) * cCa_grad_back[i]) - cCO3_grad[i] = 0.5 * ((1-sigma_cCO3) * cCO3_grad_forw[i] + \ - (1+sigma_cCO3) * cCO3_grad_back[i]) - Phi_grad[i] = 0.5 * ((1-sigma_Phi) * Phi_grad_forw [i]+ \ - (1+sigma_Phi) * Phi_grad_back[i]) - - common_helper1[i] = Phi[i]/denominator[i] - common_helper2[i] = Phi_grad[i] * (2 + denominator[i]) \ - / denominator[i] ** 2 - helper_cCa_grad[i] = dCa * (common_helper2[i] * cCa_grad[i] \ - + common_helper1[i] * cCa_laplace[i]) - helper_cCO3_grad[i] = dCO3 * (common_helper2[i] * cCO3_grad[i] \ - + common_helper1[i] * cCO3_laplace[i]) - - two_factors[i] = cCa[i] * cCO3[i] - two_factors_upp_lim[i] = min(two_factors[i],1) - two_factors_low_lim[i] = max(two_factors[i],1) - three_factors[i] = two_factors[i] * KRat - three_factors_upp_lim[i] = min(three_factors[i],1) - three_factors_low_lim[i] = max(three_factors[i],1) - - coA[i] = CA[i] * (((1 - three_factors_upp_lim[i]) ** m2) * \ - (not_too_deep[i] * not_too_shallow[i]) - nu1 * \ - (three_factors_low_lim[i] - 1) ** m1) - - coC[i] = CC[i] * (((two_factors_low_lim[i] - 1) ** n1) - nu2 * \ - (1 - two_factors_upp_lim[i]) ** n2) - - common_helper3[i] = coA[i] - lambda_* coC[i] - - dW_dx[i] = -rhorat * Phi_grad[i] * (2 * Phi[i] * F[i] + 10 * (F[i] - 1)) + sigma_Phi = np.cosh(Peclet_Phi)/np.sinh(Peclet_Phi) - 1/Peclet_Phi + else: + sigma_cCa = 0 + sigma_cCO3 = 0 + sigma_Phi = 0 + + cCa_grad[i] = 0.5 * ((1-sigma_cCa) * cCa_grad_forw[i] + \ + (1+sigma_cCa) * cCa_grad_back[i]) + cCO3_grad[i] = 0.5 * ((1-sigma_cCO3) * cCO3_grad_forw[i] + \ + (1+sigma_cCO3) * cCO3_grad_back[i]) + Phi_grad[i] = 0.5 * ((1-sigma_Phi) * Phi_grad_forw [i]+ \ + (1+sigma_Phi) * Phi_grad_back[i]) + + common_helper1[i] = Phi[i]/denominator[i] + common_helper2[i] = Phi_grad[i] * (2 + denominator[i]) \ + / denominator[i] ** 2 + helper_cCa_grad[i] = dCa * (common_helper2[i] * cCa_grad[i] \ + + common_helper1[i] * cCa_laplace[i]) + helper_cCO3_grad[i] = dCO3 * (common_helper2[i] * cCO3_grad[i] \ + + common_helper1[i] * cCO3_laplace[i]) + + two_factors[i] = cCa[i] * cCO3[i] + two_factors_upp_lim[i] = min(two_factors[i],1) + two_factors_low_lim[i] = max(two_factors[i],1) + three_factors[i] = two_factors[i] * KRat + three_factors_upp_lim[i] = min(three_factors[i],1) + three_factors_low_lim[i] = max(three_factors[i],1) + + coA[i] = CA[i] * (((1 - three_factors_upp_lim[i]) ** m2) * \ + (not_too_deep[i] * not_too_shallow[i]) - nu1 * \ + (three_factors_low_lim[i] - 1) ** m1) + + coC[i] = CC[i] * (((two_factors_low_lim[i] - 1) ** n1) - nu2 * \ + (1 - two_factors_upp_lim[i]) ** n2) + + common_helper3[i] = coA[i] - lambda_* coC[i] + + dW_dx[i] = -rhorat * Phi_grad[i] * (2 * Phi[i] * F[i] + 10 * (F[i] - 1)) - # This is dCA_dt - rate[0][i] = - U[i] * CA_grad[i] - Da * ((1 - CA[i]) \ - * coA[i] + lambda_ * CA[i] * coC[i]) - - # This is dCC_dt - rate[1][i] = - U[i] * CC_grad[i] + Da * (lambda_ * \ - (1 - CC[i]) * coC[i] + CC[i] * coA[i]) - - # This is dcCa_dt - rate[2][i] = helper_cCa_grad[i]/Phi[i] - W[i] * \ - cCa_grad[i] + Da * one_minus_Phi[i] * \ - (delta - cCa[i]) * common_helper3[i] \ - /Phi[i] - - # This is dcCO3_dt - rate[3][i] = helper_cCO3_grad[i]/Phi[i] - W[i] * \ - cCO3_grad[i] + Da * one_minus_Phi[i] * \ - (delta - cCO3[i]) * common_helper3[i] \ - /Phi[i] - - # This is dPhi_dt - rate[4][i] = - (dW_dx[i] * Phi[i] + W[i] * Phi_grad[i]) \ - + dPhi[i] * Phi_laplace[i] + Da * one_minus_Phi[i] \ - * common_helper3[i] - - return rate - - return pde_rhs \ No newline at end of file + # This is dCA_dt + rate[i] = - U[i] * CA_grad[i] - Da * ((1 - CA[i]) \ + * coA[i] + lambda_ * CA[i] * coC[i]) + + # This is dCC_dt + rate[no_depths + i] = - U[i] * CC_grad[i] + Da * (lambda_ * \ + (1 - CC[i]) * coC[i] + CC[i] * coA[i]) + + # This is dcCa_dt + rate[2 * no_depths + i] = helper_cCa_grad[i]/Phi[i] - W[i] * \ + cCa_grad[i] + Da * one_minus_Phi[i] * \ + (delta - cCa[i]) * common_helper3[i] \ + /Phi[i] + + # This is dcCO3_dt + rate[3 * no_depths + i] = helper_cCO3_grad[i]/Phi[i] - W[i] * \ + cCO3_grad[i] + Da * one_minus_Phi[i] * \ + (delta - cCO3[i]) * common_helper3[i] \ + /Phi[i] + + # This is dPhi_dt + rate[4 * no_depths + i] = - (dW_dx[i] * Phi[i] + W[i] * Phi_grad[i]) \ + + dPhi[i] * Phi_laplace[i] + Da * one_minus_Phi[i] \ + * common_helper3[i] + + return rate + + def zeros(self, t, y, progress_proxy, progress_dt, t0): + """ solve_ivp demands that I add these two extra aguments, i.e. + pbar and state, as in jac, where I need them for + tqdm progress display. + However, for this rhs calculation, they are redundant. """ + + return np.amin(y) + + def zeros_CA(self, t, y, progress_proxy, progress_dt, t0): + """ solve_ivp demands that I add these two extra aguments, i.e. + pbar and state, as in jac, where I need them for + tqdm progress display. + However, for this rhs calculation, they are redundant. """ + + return np.amin(y[self.CA_sl]) + + def zeros_CC(self, t, y, progress_proxy, progress_dt, t0): + """ solve_ivp demands that I add these two extra aguments, i.e. + pbar and state, as in jac, where I need them for + tqdm progress display. + However, for this rhs calculation, they are redundant. """ + + return np.amin(y[self.CC_sl]) + + def ones_CA_plus_CC(self, t, y, progress_proxy, progress_dt, t0): + """ solve_ivp demands that I add these two extra aguments, i.e. + pbar and state, as in jac, where I need them for + tqdm progress display. + However, for this rhs calculation, they are redundant. """ + + CA = y[self.CA_sl] + CC = y[self.CC_sl] + return np.amax(CA + CC) - 1 + + def ones_Phi(self, t, y, progress_proxy, progress_dt, t0): + """ solve_ivp demands that I add these two extra aguments, i.e. + pbar and state, as in jac, where I need them for + tqdm progress display. + However, for this rhs calculation, they are redundant. """ + + Phi = y[self.Phi_sl] + return np.amax(Phi) - 1 + + def zeros_U(self, t, y, progress_proxy, progress_dt, t0): + """ solve_ivp demands that I add these two extra aguments, i.e. + pbar and state, as in jac, where I need them for + tqdm progress display. + However, for this rhs calculation, they are redundant. """ + + Phi = y[self.Phi_sl] + F = 1 - np.exp(10 - 10 / Phi) + U = self.presum + self.rhorat * Phi ** 3 * F / (1 - Phi) + # Assume that U is positive at every depth at the beginning + # of the integration, so U will become zero first where it + # is smallest. + return np.amin(U) + + def zeros_W(self, t, y, progress_proxy, progress_dt, t0): + """ solve_ivp demands that I add these two extra aguments, i.e. + pbar and state, as in jac, where I need them for + tqdm progress display. + However, for this rhs calculation, they are redundant. """ + + Phi = y[self.Phi_sl] + F = 1 - np.exp(10 - 10 / Phi) + W = self.presum - self.rhorat * Phi ** 2 * F + # Assume that W is negative at every depth at the beginning + # of the integration, so U will become zero first where it is + # largest or least negative. + return np.amax(W) + diff --git a/marlpde/parameters.py b/marlpde/parameters.py index 5f0e5c4..021e62e 100644 --- a/marlpde/parameters.py +++ b/marlpde/parameters.py @@ -1,8 +1,4 @@ -import configparser -from pathlib import Path -from dataclasses import (dataclass, asdict, make_dataclass, fields) -from subprocess import (run) -import h5py as h5 +from dataclasses import (dataclass, make_dataclass, fields) from pint import UnitRegistry import numpy as np from scipy.sparse import lil_matrix, dia_matrix, csr_matrix @@ -118,7 +114,8 @@ def Map_Scenario(): ("n2", float, None), ("DCa", float, None), ("PhiNR", float, None), - ("N", int, None)]) + ("N", int, None), + ("FV_switch", int, None)]) def post_init(self): # The Python parameters that need additional conversion @@ -138,8 +135,12 @@ def post_init(self): self.n2 = self.n1 self.DCa = self.D0Ca self.PhiNR = self.PhiIni + # The number of cells that comprise the grid. self.N = 200 - + # It could be that F-V caused instabilities instead of resolving them, + # e.g. in the case of oscillations. + # FV_switch = 1 will use the Fiadeiro-Veronis scheme for spatial derivatives. + self.FV_switch = 1 derived_dataclass = make_dataclass("Mapped parameters", derived_fields, namespace={"__post_init__": post_init}) @@ -158,14 +159,12 @@ def jacobian_sparsity(): matrix computed by `solve_ivp` (when neither the Jacobian nor the Jacobian sparsity matrix was provided) has non-zero elements not only on the diagonals but also on adjacent positions along the - diagonals. This was confirmed by ChatGPT: - "It is indeed common for the Jacobian matrix of coupled partial - differential equations (PDEs), such as advection-diffusion equations, - to have non-zero elements beyond the main diagonal. This phenomenon, - where adjacent elements in the Jacobian are non-zero, is often due - to the spatial discretization of the differential equations." - Likewise computing a functional Jacobian is only feasible for the - diagonals, for positions near the diagonals it is very hard. + diagonals. + See equation 2.58 from "Finite Difference Methods Finite for Differential + Equations by Randall J. Leveque for an example of off-diagonal Jacobian + matrix elements. This example is about the Jacobian for solving the pde + describing the motion of a pendulum with a certain mass at the end of a + rigid (but massless) bar. Here we will choose diagonals of width 3. ''' @@ -188,7 +187,9 @@ def jacobian_sparsity(): try: assert len(offsets) == diagonals.shape[0] except AssertionError as e: + print(f'Assertion failed at the {str(e)} line. \n') print('Setup of diagonals incorrect.') + # Construct the sparse matrix. raw_matrix = lil_matrix(dia_matrix((diagonals, offsets), shape=(n, n))) # Set the Jacobian elements to zero that correspond with the derivatives @@ -203,16 +204,13 @@ class Solver(): Initialises all the parameters for the solver. So parameters like time interval, time step and tolerance. ''' - dt: float = 1.e-6 - # t_range is the integration time in units of T*. - t_range: float = 1 - 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 + first_step: float = 1e-6 + atol: float = 1e-3 + rtol: float = 1e-3 + # t_span is a tuple (begin time, end time) of integration in units of T*. + t_span: tuple = (0, 1) # solve_ivp from scipy offers six methods. They can be set here. - method: str = "LSODA" + method: str = "Radau" # 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" @@ -220,23 +218,24 @@ class Solver(): lband: int = 1 uband: int = 1 backend: str = "numba" - ret_info: bool = True - + dense_output: bool = False + jac_sparsity: csr_matrix = None + def __post_init__(self): ''' Filter out solver settings that are mutually incompatible. ''' try: - if self.solver != "scipy": - del self.__dataclass_fields__["method"] + if self.method == "LSODA": + del self.__dataclass_fields__["jac_sparsity"] + else: del self.__dataclass_fields__["lband"] del self.__dataclass_fields__["uband"] - else: - del self.__dataclass_fields__["scheme"] - del self.__dataclass_fields__["adaptive"] - if self.method != "LSODA": - del self.__dataclass_fields__["lband"] - del self.__dataclass_fields__["uband"] + + if self.method == "Radau" or self.method == "BDF": + self.jac_sparsity = jacobian_sparsity() + else: + del self.__dataclass_fields__["jac_sparsity"] except KeyError: pass @@ -244,11 +243,19 @@ def __post_init__(self): @dataclass class Tracker: ''' - Initialises all the tracking parameters, such as tracker interval. - Also indicates the quantities to be tracked, as boolean values. + Settings for tracking progress and for storing intermediate results when + solve_ivp runs. ''' - progress_tracker_interval: float = Solver().t_range / 1_000 - live_plotting: bool = False - plotting_interval: str = '0:05' - data_tracker_interval: float = 0.01 - track_U_at_bottom: bool = False + # Number of progress updates. This only affects the progress bar. + no_progress_updates: int = 100_000 + + # Number of times to evaluate, for storage. + # 2 means only initial and end values. + no_t_eval: int = 2 + + + # Array with all times that solutions from solve_ivp should be recorded. + t_eval: np.ndarray = None + + def __post_init__(self): + self.t_eval = np.linspace(*Solver().t_span, num = self.no_t_eval) diff --git a/tests/Regression_test/test_regression.py b/tests/Regression_test/test_regression.py index 216bad7..5fd174f 100644 --- a/tests/Regression_test/test_regression.py +++ b/tests/Regression_test/test_regression.py @@ -43,13 +43,13 @@ def test_integration_Scenario_A(): Scenario_parameters = asdict(Map_Scenario()) | \ {"Phi0": 0.6, "PhiIni": 0.5, "PhiNR": 0.6} # integrate_equations returns four variables, we only need the first one. - solution, _, _, _, _ = integrate_equations(asdict(Solver()), - asdict(Tracker()), - Scenario_parameters) + last_field_sol, _, _, _, _ = integrate_equations(asdict(Solver()), + asdict(Tracker()), + Scenario_parameters) # Test the final distribution of all five fields over depths - assert_allclose(solution.data, Scenario_A_data[-1, :, :], + assert_allclose(last_field_sol, Scenario_A_data[-1, :, :], rtol=rtol, atol=atol) def test_high_porosity_integration(): @@ -75,16 +75,16 @@ def test_high_porosity_integration(): # Smaller initial time step needed for these high porosities. # If not, we get a ZeroDivisionError from "wrapped_stepper". - Solver_parms = replace(Solver(), dt = 5e-7) + Solver_parms = replace(Solver(), first_step = 5e-7) # Concatenate the dict containing the Scenario parameters with the # dict containing the solver parameters (such as required tolerance). # integrate_equations returns four variables, we only need the first one. - solution, _, _, _, _ = integrate_equations(asdict(Solver_parms), - asdict(Tracker()), - Scenario_parameters) + last_field_sol, _, _, _, _ = integrate_equations(asdict(Solver_parms), + asdict(Tracker()), + Scenario_parameters) # Test the final distribution of all five fields over depths - assert_allclose(solution.data, high_porosity_data[-1, :, :], + assert_allclose(last_field_sol, high_porosity_data[-1, :, :], rtol=rtol, atol=atol) def test_cross_check_with_Matlab_output(): @@ -114,13 +114,13 @@ def test_cross_check_with_Matlab_output(): Scenario_parameters = asdict(Map_Scenario()) |\ {"Phi0": 0.5, "PhiIni": 0.5, "PhiNR": 0.5, \ "k3": 0.01, "k4": 0.01} - Xstar = Scenario_parameters["Xstar"] + Xstar = Scenario_parameters["Xstar"] max_depth = Scenario_parameters["max_depth"] - solution, _, _, _, _ = integrate_equations(asdict(Solver()), - asdict(Tracker()), - Scenario_parameters) + last_field_sol, _, _, _, _ = integrate_equations(asdict(Solver()), + asdict(Tracker()), + Scenario_parameters) Number_of_depths = Scenario_parameters["N"] @@ -138,12 +138,14 @@ def test_cross_check_with_Matlab_output(): Matlab_output_interpolated[field, :] = np.interp(Python_plotting_depths, Matlab_depths, \ Matlab_output[field, :, 0]) - # Compare the Python and Matlab output over all depths, except for close to the surface, - # since we are dealing with a boundary layer near the surface - see chapter 7, page 56 of - # Willem Hundsdorfer: "Numerical Solution of Advection-Diffusion-Reaction Equations". - # Consequently, both the Matlab and Python solutions for the concentrations and the - # porosity jump a bit up and down near the surface, due to the high Peclet numbers. - assert_allclose(solution.data[:, 2:], Matlab_output_interpolated[:, 2:], atol = atol) + # Compare the Python and Matlab output over all depths, except for close to + # the surface, since we are dealing with a boundary layer near the surface + # - see chapter 7, page 56 of Willem Hundsdorfer: "Numerical Solution of + # Advection-Diffusion-Reaction Equations". Consequently, both the Matlab and + # Python solutions for the concentrations and the porosity jump a bit up and + # down near the surface, due to the high Peclet numbers. + assert_allclose(last_field_sol[:, 2:], + Matlab_output_interpolated[:, 2:], atol = atol)