From 29faa2fee7539977f9e1c41a86441c1a722b7e94 Mon Sep 17 00:00:00 2001 From: Thomas Roy Date: Thu, 7 Mar 2024 14:54:23 -0800 Subject: [PATCH 1/5] homog_params input to define bulk reactions --- docs/source/user_guide/homog_params.rst | 86 ++++++++++ docs/source/user_guide/index.rst | 1 + docs/source/user_guide/physical_params.rst | 2 +- echemfem/solver.py | 50 ++++++ examples/carbonate_homog_params.py | 183 +++++++++++++++++++++ 5 files changed, 321 insertions(+), 1 deletion(-) create mode 100644 docs/source/user_guide/homog_params.rst create mode 100644 examples/carbonate_homog_params.py diff --git a/docs/source/user_guide/homog_params.rst b/docs/source/user_guide/homog_params.rst new file mode 100644 index 0000000..0f9265a --- /dev/null +++ b/docs/source/user_guide/homog_params.rst @@ -0,0 +1,86 @@ +Homogeneous Reaction Parameters +=============================== + +The parameters of homogeneous/bulk reactions can be provided through +:attr:`echemfem.EchemSolver.homog_params`, a list containing one dictionary for +each reaction. Below is a list of different keys that should appear in each +dictionary + +* :Key: ``"stoichiometry"`` + :Type: :py:class:`dict` + :Description: This entry defines the stoichiometry of the reaction. Each key in this dictionary should be the name of a species (as defined in the ``"name"`` key values of ``conc_params``). The corresponding value is the stoichemetry coefficient for the species in this reaction. It should be an integer: negative for a reactant, positive for a product. +* :Key: ``"forward rate constant"`` + :Type: :py:class:`float` or Firedrake expression + :Description: Rate constant for the forward reaction. +* :Key: ``"backward rate constant"`` or ``"equilibrium constant"`` + :Type: :py:class:`float` or Firedrake expression + :Description: Rate constant for the backward reaction or equilibrium constant of the reaction. Only of one these is required. +* :Key: ``"reference concentration"`` (optional) + :Type: :py:class:`float` or Firedrake expression + :Description: Value of 1M in the appropriate units. By setting this, the rate constants are assumed to have the same units (typically 1/s), and the equilibrium constant to be nondimensional. + +Assuming an arbitrary homogeneous reaction system where the forward and +backward rate constants for reaction :math:`i` are :math:`k_i^f` and +:math:`k_i^b`, respectively, and the stoichiometry coefficient of species +:math:`j` in reaction :math:`i` is :math:`s_{j,i}`, then, according to the law +of mass action, the reaction for species :math:`j` is given by + +.. math:: + + R_j = \sum_i s_{j, i} \left( k_i^f \prod_{n, s_{n, i} < 0} c_n^{-s_{n, i}} - k_i^b \prod_{n, s_{n, i} > 0} c_n^{s_{n, i}} \right), + +where :math:`c_n` is the concentration of species :math:`n`. If the equilibrium +constant :math:`K_i` is provided instead of the backward rate constant, then it +is automatically recovered via :math:`k_i^b = \dfrac{k_i^f}{K_i}`. + +Note that here the rate constants can have different units, depending on the +number of reactants or products. It is also common to write the reactions in +terms of dimensionless "activities" instead of concentrations. Then, the rate +constants have the same units (typically 1/s) and the equilibrium constants are +dimensionless. With this definition of rate constants, we instead write the +reaction for species :math:`j` as + + +.. math:: + + R_j = \sum_i s_{j, i} c_\mathrm{ref} \left( k_i^f \prod_{n, s_{n, i} < 0} a_n^{-s_{n, i}} - k_i^b \prod_{n, s_{n, i} > 0} a_n^{s_{n, i}} \right), + +where :math:`a_n = c_n / c_\mathrm{ref}` is the activity of species :math:`n` +and :math:`c_\mathrm{ref} = 1\text{M}` is a reference concentration, which +needs to be provided in the correct units to ``homog_params``. + +Here is an example of a reaction system that can be implement via +``homog_params``. Simplified CO2 - bicarbonate homogeneous reactions: + +.. math:: + + \mathrm{CO}_2 + \mathrm{OH}^- \xrightarrow[k_1^b]{k_1^f} \mathrm{HCO}_3^- + + \mathrm{HCO}_3^- + \mathrm{OH}^- \xrightarrow[k_2^b]{k_2^f} \mathrm{CO}_3^{2-} + +Here, we assume the dilute solution case where :math:`\mathrm{H}_2\mathrm{O}` +concentration is not tracked. In ``conc_params``, we will have entries for the +other species in this reactions, with ``"name"`` key values: ``"CO2"``, +``"OH"``, ``"HCO3"``, and ``"CO3"``. We get the following ``homog_params`` +list: + +.. code:: + + [{"stoichiometry": {"CO": -1, + "OH": -1, + "HCO3": 1}, + "forward rate constant": k1f, + "backward rate constant": k1b + }, + {"stoichiometry": {"HCO3": -1, + "OH": -1, + "CO3": 1}, + "forward rate constant": k2f, + "backward rate constant": k2b + }] + +As an alternative to the ``homog_params`` interface, the reactions can be +written directly and passed as a function in ``physical_params["bulk +reaction"]`` (:doc:`physical_params`). See ``examples/carbonate.py`` and +``examples/carbonate_homog_params.py`` for two equivalent implementations of +the bicarbonate bulk reactions. diff --git a/docs/source/user_guide/index.rst b/docs/source/user_guide/index.rst index 545dc3e..68f3a8a 100644 --- a/docs/source/user_guide/index.rst +++ b/docs/source/user_guide/index.rst @@ -25,6 +25,7 @@ To create a customized solver, the user should set the following inputs: physical_params boundary_conditions echem_params + homog_params Here is a barebone example of how a user might define their own solver. Several concrete examples can be found in `echemfem/examples diff --git a/docs/source/user_guide/physical_params.rst b/docs/source/user_guide/physical_params.rst index 8a723ff..d7684b8 100644 --- a/docs/source/user_guide/physical_params.rst +++ b/docs/source/user_guide/physical_params.rst @@ -103,7 +103,7 @@ Below is a list of other keys that can appear in the dictionary * :Key: ``"bulk reaction"`` :Type: a function - :Description: Homogeneous/bulk reactions to be added to the right-hand side of mass conservation equations. + :Description: Homogeneous/bulk reactions to be added to the right-hand side of mass conservation equations. These can instead be set using ``homog_params`` (:doc:`homog_params`). Args: u: solution state. The value of the different concentrations can be recovered through ``u([self.i_c["species name"])`` within a ``echemfem.EchemSolver`` object. diff --git a/echemfem/solver.py b/echemfem/solver.py index b6d12f1..ac35ddd 100644 --- a/echemfem/solver.py +++ b/echemfem/solver.py @@ -60,6 +60,7 @@ def __init__( physical_params, mesh, echem_params=[], + homog_params=[], gas_params=[], stats_file="", overwrite_stats_file=False, @@ -72,6 +73,7 @@ def __init__( self.physical_params = physical_params self.conc_params = conc_params self.echem_params = echem_params + self.homog_params = homog_params self.gas_params = gas_params self.num_c = len(conc_params) self.num_g = len(gas_params) @@ -380,6 +382,9 @@ def _internal_dS(subdomain_id=None, domain=None): elif self.flow["advection"]: self.set_velocity() + if self.homog_params: + self.setup_bulk_reactions() + self.setup_forms(us, v) self.u = u @@ -1205,6 +1210,51 @@ def init_solver_parameters( solver_parameters.update(pc) self.solver_parameters = solver_parameters + def setup_bulk_reactions(self): + """ Setting up the homogeneous bulk reactions using homog_params + + Setup using the law of mass action + """ + def bulk_reaction(u): + reactions = [0] * self.num_liquid + reaction_f = [0] * len(self.homog_params) + reaction_b = [0] * len(self.homog_params) + for idx, reaction in enumerate(self.homog_params): + reaction_f[idx] = reaction["forward rate constant"] + if reaction.get("backward rate constant"): + reaction_b[idx] = reaction["backward rate constant"] + else: + reaction_b[idx] = reaction["forward rate constant"] / \ + reaction["equilibrium constant"] + if reaction.get("reference concentration"): + c_ref = reaction["reference concentration"] + else: + c_ref = 1.0 + for name in reaction["stoichiometry"]: + s = reaction["stoichiometry"][name] + a = u[self.i_c[name]] / c_ref + if s < 0: + reaction_f[idx] *= a ** (-s) + if s > 0: + reaction_b[idx] *= a ** s + + for i in self.idx_c: + name = self.conc_params[i]["name"] + for idx, reaction in enumerate(self.homog_params): + if name in reaction["stoichiometry"]: + s = reaction["stoichiometry"][name] + if reaction.get("reference concentration"): + c_ref = reaction["reference concentration"] + else: + c_ref = 1.0 + reactions[i] += s * c_ref * (reaction_f[idx] - reaction_b[idx]) + + return reactions + if self.physical_params.get("bulk reaction"): + print = PETSc.Sys.Print + print("WARNING: Overwriting bulk reaction from physical_params with the one in homog_params") + self.physical_params["bulk reaction"] = bulk_reaction + def effective_diffusion(self, D, phase="liquid"): """ Bruggeman correlation for effective diffusion in a porous medium """ diff --git a/examples/carbonate_homog_params.py b/examples/carbonate_homog_params.py new file mode 100644 index 0000000..aff74f0 --- /dev/null +++ b/examples/carbonate_homog_params.py @@ -0,0 +1,183 @@ +from firedrake import * +from echemfem import EchemSolver, IntervalBoundaryLayerMesh + + +class CarbonateSolver(EchemSolver): + """ + Steady-state version of example from + Gupta, N., Gattrell, M. and MacDougall, B., 2006. Calculation for the + cathode surface concentrations in the electrochemical reduction of CO2 in + KHCO3 solutions. Journal of applied electrochemistry, 36(2), pp.161-172. + but using the bicarbonate bulk reactions from + Schulz, K.G., Riebesell, U., Rost, B., Thoms, S. and Zeebe, R.E., 2006. + Determination of the rate constants for the carbon dioxide to bicarbonate + inter-conversion in pH-buffered seawater systems. Marine chemistry, + 100(1-2), pp.53-65. + """ + + def __init__(self): + + delta = 0.0001 + mesh = IntervalBoundaryLayerMesh(100, delta, 100, 1e-7) + + # Reaction Rates in SI units + k1f = 2.23 # m3/mol.s + k1r = 9.71e-5 # 1/s + k2f = 3.71e-2 # 1/s + k2r = 2.67e1 # m3/mol.s + k3f = 3.06e5 # 1/s + k3r = 6.0e6 # m3/mol.s + k4f = 5.0e7 # m3/mol.s + k4r = 59.44 # 1/s + k5f = 2.31e7 # m3/mols. (Note - typo in Schulz) + k5r = 1.4 # mol/m3.s + # Bulk Conditions + C_HCO3_bulk = 500. # mol/m3 + C_CO32_bulk = 6.5 # mol/m3 + C_OH_bulk = k3f / k3r * C_CO32_bulk / C_HCO3_bulk # mol/m3 + C_H_bulk = (k5r / k5f) / C_OH_bulk # mol/m3 + C_CO2_bulk = (k1r / k1f) * C_HCO3_bulk / C_OH_bulk # mol/m3 + + homog_params = [] + + homog_params.append({"stoichiometry": {"CO2": -1, + "OH": -1, + "HCO3": 1}, + "forward rate constant": k1f, + "backward rate constant": k1r + }) + + homog_params.append({"stoichiometry": {"CO2": -1, + "H": 1, + "HCO3": 1}, + "forward rate constant": k2f, + "backward rate constant": k2r + }) + + homog_params.append({"stoichiometry": {"CO3": -1, + "OH": 1, + "HCO3": 1}, + "forward rate constant": k3f, + "backward rate constant": k3r + }) + + homog_params.append({"stoichiometry": {"CO3": -1, + "H": -1, + "HCO3": 1}, + "forward rate constant": k4f, + "backward rate constant": k4r + }) + + homog_params.append({"stoichiometry": {"OH": -1, + "H": -1}, + "forward rate constant": k5f, + "backward rate constant": k5r + }) + + conc_params = [] + + conc_params.append({"name": "CO2", + "diffusion coefficient": 19.1e-10, # m^2/s + "bulk": C_CO2_bulk, # mol/m3 + }) + + conc_params.append({"name": "OH", + "diffusion coefficient": 54.0e-10, # m^2/s + "bulk": C_OH_bulk, # mol/m3 + }) + + conc_params.append({"name": "H", + "diffusion coefficient": 96.0e-10, # m^2/s + "bulk": C_H_bulk, # mol/m3 + }) + + conc_params.append({"name": "CO3", + "diffusion coefficient": 7.0e-10, # m^2/s + "bulk": C_CO32_bulk, # mol/m3 + }) + + conc_params.append({"name": "HCO3", + "diffusion coefficient": 9.4e-10, # m^2/s + "bulk": C_HCO3_bulk, # mol/m3 + }) + + physical_params = {"flow": ["diffusion"], + "F": 96485., # C/mol + } + + super().__init__(conc_params, physical_params, mesh, family="CG", homog_params=homog_params) + + def neumann(self, C, conc_params, u): + name = conc_params["name"] + if name in ["HCO3", "CO3", "H"]: + return Constant(0) + j = 50. # A/m^2 + F = self.physical_params["F"] + zeffCH4 = 8. + zeffC2H4 = 12. + zeffCO = 2. + zeffHCOO = 2. + zeffH2 = 2. + cefCH4 = 0.25 + cefC2H4 = 0.20 + cefCO = 0.05 + cefHCOO = 0.10 + cefH2 = 0.40 + if name == "CO2": + print(-(j / F) * (cefHCOO / zeffHCOO + cefCO / zeffCO + + cefCH4 / zeffCH4 + 2 * cefC2H4 / zeffC2H4)) + return -(j / F) * (cefHCOO / zeffHCOO + cefCO / zeffCO + + cefCH4 / zeffCH4 + 2 * cefC2H4 / zeffC2H4) + if name == "OH": + return (j / F) * (cefHCOO / zeffHCOO + 2 * cefCO / zeffCO + + 8 * cefCH4 / zeffCH4 + 12 * cefC2H4 / zeffC2H4 + + 2 * cefH2 / zeffH2) # note error in Gupta et al. + + def set_boundary_markers(self): + self.boundary_markers = {"bulk dirichlet": (1,), # U_solid = U_app, C = C_0 + "neumann": (2,), + } + + +solver = CarbonateSolver() +solver.setup_solver() +solver.solve() +C_CO2, C_OH, C_H, C_CO3, C_HCO3 = solver.u.subfunctions +# OH boundary layer +x = solver.mesh.coordinates +C_OH_bl = Function(solver.V).assign(C_OH).dat.data[100:] +x_bl = x.dat.data[100:] + +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +from matplotlib.ticker import FormatStrFormatter +filename = "carbonate.png" +fig = plt.figure(constrained_layout=True, figsize=(16, 8)) +spec = gridspec.GridSpec(ncols=3, nrows=2, figure=fig) +ax1 = fig.add_subplot(spec[0, 0]) +ax2 = fig.add_subplot(spec[0, 1]) +ax3 = fig.add_subplot(spec[1, 0]) +ax4 = fig.add_subplot(spec[1, 1]) +ax5 = fig.add_subplot(spec[0, 2]) +ax6 = fig.add_subplot(spec[1, 2]) + +plot(C_CO2, axes=ax1) +ax1.set(xlabel='distance (m)', + ylabel='CO$_2$ concentration (M)') +plot(C_HCO3, axes=ax2) +ax2.set(xlabel='distance (m)', + ylabel='HCO$_3$ concentration (M)') +plot(C_CO3, axes=ax3) +ax3.set(xlabel='distance (m)', + ylabel='CO$_3$ concentration (M)') +plot(C_OH, axes=ax4) +ax4.set(xlabel='distance (m)', + ylabel='OH concentration (M)') +plot(C_H, axes=ax5) +ax5.set(xlabel='distance (m)', + ylabel='H concentration (M)') +plt.plot(x_bl, C_OH_bl, axes=ax6, color='k', linewidth=2) +ax6.set(xlabel='distance (m)', + ylabel='OH concentration (M)') + +plt.savefig(filename) From d0d377bb278ebca360217f4e2ad3ebe03552e5bf Mon Sep 17 00:00:00 2001 From: Thomas Roy Date: Fri, 8 Mar 2024 09:42:36 -0800 Subject: [PATCH 2/5] fixed SMPNP --- echemfem/solver.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/echemfem/solver.py b/echemfem/solver.py index ac35ddd..dec1507 100644 --- a/echemfem/solver.py +++ b/echemfem/solver.py @@ -1574,15 +1574,16 @@ def mass_conservation_form( NA = self.physical_params["Avogadro constant"] phi = 0.0 a0 = 2.3e-10 # Size of water molecule + aj = conc_params.get("solvated diameter") + betaj = (aj/a0)**3 for i in range(n_c): ai = self.conc_params[i].get("solvated diameter") if ai is not None and ai != 0.0: - betai = (ai/a0)**3 - phi += betai * NA * ai ** 3 * u[i] + phi += NA * ai ** 3 * u[i] # size-modified (SMPNP) correction (cf. Eq. 4 in doi:10.26434/chemrxiv-2022-h2mrp) - mu_ex = -ln(1 - phi) + mu_ex = - betaj * ln(1 - phi) a += D * C * inner(grad(ln(C) + mu_ex), grad(test_fn)) * self.dx() From 5168a31c9645b985e42c6afc1def1b5b5e4e27e4 Mon Sep 17 00:00:00 2001 From: Thomas Roy Date: Thu, 21 Mar 2024 13:46:00 -0700 Subject: [PATCH 3/5] Add test for transient solver --- tests/test_heat_equation.py | 83 +++++++++++++++++++++++++++++++++++++ 1 file changed, 83 insertions(+) create mode 100644 tests/test_heat_equation.py diff --git a/tests/test_heat_equation.py b/tests/test_heat_equation.py new file mode 100644 index 0000000..4eb3427 --- /dev/null +++ b/tests/test_heat_equation.py @@ -0,0 +1,83 @@ +import pytest +from firedrake import * +from echemfem import TransientEchemSolver +import numpy as np + + +class HeatEquationSolver(TransientEchemSolver): + def __init__(self, N, family="DG"): + mesh = UnitSquareMesh(N, N, quadrilateral=True) + conc_params = [] + + self.time = Constant(0) # for now, need to define this here + t = self.time + x, y = SpatialCoordinate(mesh) + C1ex = (sin(x) + cos(y)) * exp(t) + C2ex = (cos(x) + sin(y)) * exp(-t) + self.C1ex = C1ex + self.C2ex = C2ex + + def f(C): + f1 = C1ex - div(grad(C1ex)) + f2 = -C2ex - div(grad(C2ex)) + return [f1, f2] + + conc_params.append({"name": "C1", + "diffusion coefficient": 1.0, + "bulk": C1ex, + }) + + conc_params.append({"name": "C2", + "diffusion coefficient": 1.0, + "bulk": C2ex, + }) + physical_params = {"flow": ["diffusion"], + "bulk reaction": f, + } + + super().__init__(conc_params, physical_params, mesh, family=family) + + def set_boundary_markers(self): + self.boundary_markers = {"bulk dirichlet": (1, 2, 3, 4,), + } + + +def test_convergence(): + err_old = 1e6 + for i in range(4): + solver = HeatEquationSolver(2**(i + 1)) + solver.setup_solver() + times = np.linspace(0, 1, 1+2**(2*(i+1))) + solver.solve(times) + c1, c2 = solver.u.subfunctions + err = errornorm(solver.C1ex, c1) + errornorm(solver.C2ex, c2) + assert err < 0.27 * err_old + err_old = err + + +def test_convergence_CG(): + err_old = 1e6 + for i in range(5): + solver = HeatEquationSolver(2**(i + 1), family="CG") + solver.setup_solver() + times = np.linspace(0, 1, 1+2**(2*i)) + solver.solve(times) + c1, c2 = solver.u.subfunctions + err = errornorm(solver.C1ex, c1) + errornorm(solver.C2ex, c2) + assert err < 0.26 * err_old + err_old = err + + +def test_convergence_BE(): + err_old = 1e6 + solver = HeatEquationSolver(16, family="CG") + for i in range(5): + solver.time.assign(0) + solver.setup_solver() + solver.u_old.assign(solver.u) + times = np.linspace(0, 1, 1+2**(i+1)) + solver.solve(times) + c1, c2 = solver.u.subfunctions + err = errornorm(solver.C1ex, c1) + errornorm(solver.C2ex, c2) + assert err < 0.52 * err_old + err_old = err From 89f1eebbf0997b1ba0c1bcead0e9c0ec6918ff41 Mon Sep 17 00:00:00 2001 From: Thomas Roy Date: Wed, 3 Apr 2024 09:33:13 -0700 Subject: [PATCH 4/5] bulk reactions for eliminated conc --- echemfem/solver.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/echemfem/solver.py b/echemfem/solver.py index dec1507..2f02965 100644 --- a/echemfem/solver.py +++ b/echemfem/solver.py @@ -435,7 +435,7 @@ def steady_forms(self, us, v): if gas_params: Form += self.dissolutions(us, v, conc_params, gas_params) - # TODO: bulk_reactions and dissolutions for eliminated concentrations + # TODO: dissolutions for eliminated concentrations # charge conservation if self.flow["electroneutrality"]: a, bc = self.charge_conservation_form(us, @@ -1216,7 +1216,7 @@ def setup_bulk_reactions(self): Setup using the law of mass action """ def bulk_reaction(u): - reactions = [0] * self.num_liquid + reactions = [0] * self.num_c reaction_f = [0] * len(self.homog_params) reaction_b = [0] * len(self.homog_params) for idx, reaction in enumerate(self.homog_params): From 3f20368fecc6c656bf3682dae7751a597c4724c0 Mon Sep 17 00:00:00 2001 From: Thomas Roy Date: Wed, 3 Apr 2024 15:08:40 -0700 Subject: [PATCH 5/5] Remove untested Darcy unit tests --- ...est_advection_diffusion_migration_darcy.py | 121 ---------- ...tion_diffusion_migration_darcy_twophase.py | 212 ------------------ 2 files changed, 333 deletions(-) delete mode 100644 tests/test_advection_diffusion_migration_darcy.py delete mode 100644 tests/test_advection_diffusion_migration_darcy_twophase.py diff --git a/tests/test_advection_diffusion_migration_darcy.py b/tests/test_advection_diffusion_migration_darcy.py deleted file mode 100644 index a4c5628..0000000 --- a/tests/test_advection_diffusion_migration_darcy.py +++ /dev/null @@ -1,121 +0,0 @@ -import pytest -from firedrake import * -from echemfem import EchemSolver - - -class AdvectionDiffusionMigrationSolver(EchemSolver): - def __init__(self, N, vavg): - mesh = UnitSquareMesh(N, N, quadrilateral=True) - conc_params = [] - - x, y = SpatialCoordinate(mesh) - C1ex = cos(x) + sin(y) + 3 - C2ex = C1ex - U1ex = sin(x) + cos(y) + 3 - U2ex = sin(x) + cos(y) + 3 - pex = (0.5 * cos(x) + sin(y) + 2) * vavg - vel = -grad(pex) - self.C1ex = C1ex - self.C2ex = C2ex - self.U1ex = U1ex - self.U2ex = U2ex - self.pex = pex - D1 = 0.5 - D2 = 1.0 - z1 = 2.0 - z2 = -2.0 - K = 1.0 - - def f(C): - f1 = div((vel - self.effective_diffusion(D1) * z1 * grad(U1ex)) * C1ex) - \ - div(self.effective_diffusion(D1) * grad(C1ex)) - f2 = div((vel - self.effective_diffusion(D2) * z2 * grad(U1ex)) * C2ex) - \ - div(self.effective_diffusion(D2) * grad(C2ex)) - f3 = div(- self.effective_diffusion(K, phase="solid") * grad(U2ex)) - return [f1, f2, f3] - - def fp(): - return -div(grad(pex) + self.effective_diffusion(D1) * z1 * grad(U1ex) - * C1ex + self.effective_diffusion(D2) * z2 * grad(U1ex) * C2ex) - - conc_params.append({"name": "C1", - "diffusion coefficient": 0.5, - "z": 2., - "bulk": C1ex, - "molar mass": 1.0, - }) - - conc_params.append({"name": "C2", - "diffusion coefficient": 1.0, - "z": -2., - "bulk": C2ex, - "molar mass": 1.0, - }) - physical_params = {"flow": ["diffusion", "advection", "migration", "electroneutrality", "porous", "darcy"], - "F": 1.0, # C/mol - "R": 1.0, # J/K/mol - "T": 1.0, # K - "U_app": U1ex, # V - "bulk reaction": f, - "v_avg": vavg, - "porosity": 0.5, - "solid conductivity": 1., - "specific surface area": 1., - "standard potential": 0., - "p_gas": pex, - "permeability": 1.0, - "liquid viscosity": 1.0, - "liquid density": 1.0, - "pressure source": fp, # only for testing - } - - super().__init__(conc_params, physical_params, mesh) - - def set_boundary_markers(self): - self.boundary_markers = {"inlet": (1, 2, 3, 4,), - "bulk dirichlet": (1, 2, 3, 4,), - "applied": (1, 2, 3, 4,), - "outlet": (1, 2, 3, 4,), - "liquid applied": (1, 2, 3, 4), - "gas": (1, 2, 3, 4,), - } - - -def test_convergence_low_peclet(): - err_old = 1e6 - errp_old = 1e6 - errC_old = 1e6 - errU1_old = 1e6 - errU2_old = 1e6 - for i in range(6): - solver = AdvectionDiffusionMigrationSolver(2**(i + 1), 1.0) - solver.setup_solver() - solver.solve() - c1, U1, U2, p = solver.u.subfunctions - err = errornorm(solver.C1ex, - c1) + errornorm(solver.U1ex, - U1) + errornorm(solver.U2ex, - U2) + errornorm(solver.pex, - p) - assert err < 0.29 * err_old - err_old = err - - -def test_convergence_high_peclet(): - err_old = 1e6 - errp_old = 1e6 - errC_old = 1e6 - errU1_old = 1e6 - errU2_old = 1e6 - for i in range(3, 6): - solver = AdvectionDiffusionMigrationSolver(2**(i + 1), 1e3) - solver.setup_solver() - solver.solve() - c1, U1, U2, p = solver.u.subfunctions - err = errornorm(solver.C1ex, - c1) + errornorm(solver.U1ex, - U1) + errornorm(solver.U2ex, - U2) + errornorm(solver.pex, - p) - assert err < 0.35 * err_old - err_old = err diff --git a/tests/test_advection_diffusion_migration_darcy_twophase.py b/tests/test_advection_diffusion_migration_darcy_twophase.py deleted file mode 100644 index 47d2bca..0000000 --- a/tests/test_advection_diffusion_migration_darcy_twophase.py +++ /dev/null @@ -1,212 +0,0 @@ -import pytest -from firedrake import * -from echemfem import EchemSolver - - -class AdvectionDiffusionMigrationSolver(EchemSolver): - def __init__(self, N, vavg): - mesh = UnitSquareMesh(N, N, quadrilateral=True) - conc_params = [] - - x, y = SpatialCoordinate(mesh) - C1ex = cos(x) + sin(y) + 3 - C2ex = C1ex - X1ex = 0.2 * (sin(x) + cos(y)) + 0.5 - X2ex = 1 - X1ex - U1ex = sin(x) + cos(y) + 3 - U2ex = sin(x) + cos(y) + 3 - pex = (0.5 * cos(x) + sin(y) + 2) * vavg - pgex = pex - # pgex = cos(x) + 0.5 * sin(y) + 3 - Sex = self.saturation(pex, pgex) - krl = self.relative_permeability(Sex) - krg = self.relative_permeability(1 - Sex) - vel = -krl * grad(pex) - velg = -krg * grad(pgex) - self.C1ex = C1ex - self.C2ex = C2ex - self.X1ex = X1ex - self.X2ex = X2ex - self.U1ex = U1ex - self.U2ex = U2ex - self.pex = pex - self.pgex = pgex - D1 = 0.5 - D2 = 1.0 - z1 = 2.0 - z2 = -2.0 - K = 1.0 - Mn = 1 / (X1ex + X2ex / 2.) - rhog = Mn * pgex - sumyD = X2ex / 2.0 - D1M = (1. - X1ex) / sumyD / Mn - sumyD = X1ex / 1.0 - D2M = (1. - X2ex) / sumyD / Mn - D1K = 2. * 0.1 / 3. * sqrt(8. / pi / 1.) - D2K = 2. * 0.1 / 3. * sqrt(8. / pi / 2.) - DX1 = 1 / (1 / D1M + 1 / D1K) - DX2 = 1 / (1 / D2M + 1 / D2K) - - def f(C): - f1 = div((vel - self.effective_diffusion(D1) * z1 * grad(U1ex)) * C1ex) - \ - div(self.effective_diffusion(D1) * grad(C1ex)) - f2 = div((vel - self.effective_diffusion(D2) * z2 * grad(U1ex)) * C2ex) - \ - div(self.effective_diffusion(D2) * grad(C2ex)) - f3 = div(- self.effective_diffusion(K, phase="solid") * grad(U2ex)) - return [f1, f2, f3] - - def fg(X): - f1 = div(velg * X1ex * rhog) - div(rhog * self.effective_diffusion(DX1, - phase="gas") * (grad(X1ex) + X1ex * grad(Mn) / Mn)) - f2 = div(velg * X2ex * rhog) - div(rhog * self.effective_diffusion(DX2, - phase="gas") * (grad(X2ex) + X2ex * grad(Mn) / Mn)) - return [f1, f2] - - def fp(): - return -div(krl * grad(pex) + D1 * z1 * grad(U1ex) - * C1ex + D2 * z2 * grad(U1ex) * C2ex) - - # def fpg(): - # return -div(krg * rhog * grad(pgex)) - - conc_params.append({"name": "C1", - "diffusion coefficient": 0.5, - "z": 2., - "bulk": C1ex, - "molar mass": 1.0, - }) - - conc_params.append({"name": "C2", - "diffusion coefficient": 1.0, - "z": -2., - "bulk": C2ex, - "molar mass": 1.0, - }) - physical_params = {"flow": ["diffusion", "advection", "migration", "electroneutrality", "porous", "darcy"], - "F": 1.0, # C/mol - "R": 1.0, # J/K/mol - "T": 1.0, # K - "U_app": U1ex, # V - "bulk reaction": f, - "porosity": 0.5, - "solid conductivity": 1., - "specific surface area": 1., - "standard potential": 0., - "p_gas": pex, - "permeability": 1.0, - "liquid viscosity": 1.0, - "liquid density": 1.0, - "pressure source": fp, # only for testing - # "gas pressure source": fpg, #only for testing - "gas source": fg, # only for testing - "gas viscosity": 1.0, - "average pore radius": 0.1, - } - gas_params = [] - gas_params.append({"name": "X1", - "diffusion coefficient": {"X2": 1.0}, - "gas": X1ex, - "molar mass": 1.0, - }) - - gas_params.append({"name": "X2", - "diffusion coefficient": {"X1": 1.0}, - "gas": X2ex, - "n": 2, - "molar mass": 2.0, - }) - - super().__init__(conc_params, physical_params, mesh, gas_params=gas_params) - - def set_boundary_markers(self): - self.boundary_markers = {"inlet": (1, 2, 3, 4,), - "applied": (1, 2, 3, 4,), - "outlet": (1, 2, 3, 4,), - "bulk dirichlet": (1, 2, 3, 4), - "liquid applied": (1, 2, 3, 4), - "gas": (1, 2, 3, 4), - "gas inlet": (1, 2, 3, 4), - "gas outlet": (1, 2, 3, 4), - } - - def saturation(self, pl, pg): - # eye-ball approximation of the saturation curve from Weng Weber - S_wr = 0.25 # Residual saturation of water - S_gr = 0.01 # Residual saturation of gas - c = 5e-2 - p0 = 0 - pc = pl - pg # capillary pressure - return ((1. - S_gr) - S_wr) * tanh(c * (pc - p0)) / 2. + 0.64 - # return 0.64 - - def set_gas_density(self, X, pg, gas_params): - R = self.physical_params["R"] - T = self.physical_params["T"] - - # average molar mass of the mixture - Mn = 0. - Xj = 1.0 - for i in range(self.num_g): - if gas_params[i].get("eliminated") is None: - Mn += X[i] / gas_params[i]["molar mass"] - Xj -= X[i] - else: - j = i - Mn += Xj / gas_params[j]["molar mass"] - Mn = 1. / Mn - self.Mn = Mn - self.Xj = Xj - - # ideal gas law - self.rhog = Mn * pg / R / T - if self.flow["diffusion"] and ( - self.boundary_markers.get("gas") is not None): - Mn_gas = 0.0 - for i in range(self.num_g): - Mn_gas += gas_params[i]["gas"] / gas_params[i]["molar mass"] - Mn_gas = 1. / Mn_gas - self.Mn_gas = Mn_gas - - -def test_convergence_low_peclet(): - err_old = 1e6 - errp_old = 1e6 - errC_old = 1e6 - errU1_old = 1e6 - errU2_old = 1e6 - for i in range(6): - solver = AdvectionDiffusionMigrationSolver(2**(i + 1), 1.0) - solver.setup_solver() - solver.solve() - c1, X1, U1, U2, p, pg = solver.u.subfunctions - err = errornorm(solver.C1ex, - c1) + errornorm(solver.U1ex, - U1) + errornorm(solver.U2ex, - U2) + errornorm(solver.pex, - p) + errornorm(solver.X1ex, - X1) + errornorm(solver.pgex, - pg) - assert err < 0.29 * err_old - err_old = err - - -def test_convergence_high_peclet(): - err_old = 1e6 - errp_old = 1e6 - errC_old = 1e6 - errU1_old = 1e6 - errU2_old = 1e6 - for i in range(4, 6): - solver = AdvectionDiffusionMigrationSolver(2**(i + 1), 1e2) - solver.setup_solver() - solver.solve() - c1, X1, U1, U2, p, pg = solver.u.subfunctions - err = errornorm(solver.C1ex, - c1) + errornorm(solver.U1ex, - U1) + errornorm(solver.U2ex, - U2) + errornorm(solver.pex, - p) + errornorm(solver.X1ex, - X1) + errornorm(solver.pgex, - pg) - assert err < 0.29 * err_old - err_old = err