diff --git a/examples/surface_reaction_example.py b/examples/surface_reaction_example.py new file mode 100644 index 000000000..2caef8f03 --- /dev/null +++ b/examples/surface_reaction_example.py @@ -0,0 +1,169 @@ +import festim as F +import numpy as np + +import dolfinx.fem as fem +import ufl + + +class FluxFromSurfaceReaction(F.SurfaceFlux): + def __init__(self, reaction: F.SurfaceReactionBC): + super().__init__( + F.Species(), # just a dummy species here + reaction.subdomain, + ) + self.reaction = reaction.flux_bcs[0] + + def compute(self, ds): + self.value = fem.assemble_scalar( + fem.form(self.reaction.value_fenics * ds(self.surface.id)) + ) + self.data.append(self.value) + + +my_model = F.HydrogenTransportProblem() +my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 1000)) +my_mat = F.Material(name="mat", D_0=1, E_D=0) +vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat) +left = F.SurfaceSubdomain1D(id=1, x=0) +right = F.SurfaceSubdomain1D(id=2, x=1) + +my_model.subdomains = [vol, left, right] + +H = F.Species("H") +D = F.Species("D") +my_model.species = [H, D] + +my_model.temperature = 500 + +surface_reaction_hd = F.SurfaceReactionBC( + reactant=[H, D], + gas_pressure=0, + k_r0=0.01, + E_kr=0, + k_d0=0, + E_kd=0, + subdomain=right, +) + +surface_reaction_hh = F.SurfaceReactionBC( + reactant=[H, H], + gas_pressure=lambda t: ufl.conditional(ufl.gt(t, 1), 2, 0), + k_r0=0.02, + E_kr=0, + k_d0=0.03, + E_kd=0, + subdomain=right, +) + +surface_reaction_dd = F.SurfaceReactionBC( + reactant=[D, D], + gas_pressure=0, + k_r0=0.01, + E_kr=0, + k_d0=0, + E_kd=0, + subdomain=right, +) + +my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=2, species=H), + F.DirichletBC(subdomain=left, value=2, species=D), + surface_reaction_hd, + surface_reaction_hh, + surface_reaction_dd, +] + +H_flux_right = F.SurfaceFlux(H, right) +H_flux_left = F.SurfaceFlux(H, left) +D_flux_right = F.SurfaceFlux(D, right) +D_flux_left = F.SurfaceFlux(D, left) +HD_flux = FluxFromSurfaceReaction(surface_reaction_hd) +HH_flux = FluxFromSurfaceReaction(surface_reaction_hh) +DD_flux = FluxFromSurfaceReaction(surface_reaction_dd) +my_model.exports = [ + F.XDMFExport("test.xdmf", H), + H_flux_left, + H_flux_right, + D_flux_left, + D_flux_right, + HD_flux, + HH_flux, + DD_flux, +] + +my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=5, transient=True) + +my_model.settings.stepsize = 0.1 + +my_model.initialise() +my_model.run() + +import matplotlib.pyplot as plt + +plt.stackplot( + H_flux_left.t, + np.abs(H_flux_left.data), + np.abs(D_flux_left.data), + labels=["H_in", "D_in"], +) +plt.stackplot( + H_flux_right.t, + -np.abs(H_flux_right.data), + -np.abs(D_flux_right.data), + labels=["H_out", "D_out"], +) +plt.legend() +plt.xlabel("Time (s)") +plt.ylabel("Flux (atom/m^2/s)") +plt.figure() +plt.stackplot( + HD_flux.t, + np.abs(HH_flux.data), + np.abs(HD_flux.data), + np.abs(DD_flux.data), + labels=["HH", "HD", "DD"], +) +plt.legend(reverse=True) +plt.xlabel("Time (s)") +plt.ylabel("Flux (molecule/m^2/s)") + + +plt.figure() +plt.plot(H_flux_right.t, -np.array(H_flux_right.data), label="from gradient (H)") +plt.plot( + H_flux_right.t, + 2 * np.array(HH_flux.data) + np.array(HD_flux.data), + linestyle="--", + label="from reaction rates (H)", +) + +plt.plot(D_flux_right.t, -np.array(D_flux_right.data), label="from gradient (D)") +plt.plot( + D_flux_right.t, + 2 * np.array(DD_flux.data) + np.array(HD_flux.data), + linestyle="--", + label="from reaction rates (D)", +) +plt.xlabel("Time (s)") +plt.ylabel("Flux (atom/m^2/s)") +plt.legend() +plt.show() + +# check that H_flux_right == 2*HH_flux + HD_flux +H_flux_from_gradient = -np.array(H_flux_right.data) +H_flux_from_reac = 2 * np.array(HH_flux.data) + np.array(HD_flux.data) +assert np.allclose( + H_flux_from_gradient, + H_flux_from_reac, + rtol=0.5e-2, + atol=0.005, +) +# check that D_flux_right == 2*DD_flux + HD_flux +D_flux_from_gradient = -np.array(D_flux_right.data) +D_flux_from_reac = 2 * np.array(DD_flux.data) + np.array(HD_flux.data) +assert np.allclose( + D_flux_from_gradient, + D_flux_from_reac, + rtol=0.5e-2, + atol=0.005, +) diff --git a/src/festim/__init__.py b/src/festim/__init__.py index 0e833926f..db522da46 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -18,6 +18,9 @@ ) from .boundary_conditions.flux_bc import FluxBCBase, HeatFluxBC, ParticleFluxBC from .boundary_conditions.henrys_bc import HenrysBC +from .boundary_conditions.flux_bc import FluxBCBase, ParticleFluxBC, HeatFluxBC +from .boundary_conditions.surface_reaction import SurfaceReactionBC + from .boundary_conditions.sieverts_bc import SievertsBC from .exports.average_surface import AverageSurface from .exports.average_volume import AverageVolume diff --git a/src/festim/boundary_conditions/__init__.py b/src/festim/boundary_conditions/__init__.py index de6d9f19f..f229eed48 100644 --- a/src/festim/boundary_conditions/__init__.py +++ b/src/festim/boundary_conditions/__init__.py @@ -1,4 +1,11 @@ from .dirichlet_bc import FixedConcentrationBC, FixedTemperatureBC from .flux_bc import HeatFluxBC, ParticleFluxBC +from .surface_reaction import SurfaceReactionBC -__all__ = ["FixedConcentrationBC", "FixedTemperatureBC", "ParticleFluxBC", "HeatFluxBC"] +__all__ = [ + "FixedConcentrationBC", + "FixedTemperatureBC", + "ParticleFluxBC", + "SurfaceReactionBC", + "HeatFluxBC", +] diff --git a/src/festim/boundary_conditions/surface_reaction.py b/src/festim/boundary_conditions/surface_reaction.py new file mode 100644 index 000000000..0bde33e6c --- /dev/null +++ b/src/festim/boundary_conditions/surface_reaction.py @@ -0,0 +1,127 @@ +from festim.boundary_conditions import ParticleFluxBC +from festim import k_B +from dolfinx import fem +import ufl + + +class SurfaceReactionBCpartial(ParticleFluxBC): + """Boundary condition representing a surface reaction + A + B <-> C + where A, B are the reactants and C is the product + the forward reaction rate is K_r = k_r0 * exp(-E_kr / (k_B * T)) + and the backward reaction rate is K_d = k_d0 * exp(-E_kd / (k_B * T)) + The reaction rate is: + K = K_r * C_A * C_B - K_d * P_C + with C_A, C_B the concentration of species A and B, + P_C the partial pressure of species C at the surface. + + This class is used to create the flux of a single species entering the surface + Example: The flux of species A entering the surface is K. + + + Args: + reactant (list): list of F.Species objects representing the reactants + gas_pressure (float, callable): the partial pressure of the product species + k_r0 (float): the pre-exponential factor of the forward reaction rate + E_kr (float): the activation energy of the forward reaction rate (eV) + k_d0 (float): the pre-exponential factor of the backward reaction rate + E_kd (float): the activation energy of the backward reaction rate (eV) + subdomain (F.SurfaceSubdomain): the surface subdomain where the reaction occurs + species (F.Species): the species to which the flux is applied + """ + + def __init__( + self, + reactant, + gas_pressure, + k_r0, + E_kr, + k_d0, + E_kd, + subdomain, + species, + ): + self.reactant = reactant + self.gas_pressure = gas_pressure + self.k_r0 = k_r0 + self.E_kr = E_kr + self.k_d0 = k_d0 + self.E_kd = E_kd + super().__init__(subdomain=subdomain, value=None, species=species) + + def create_value_fenics(self, mesh, temperature, t: fem.Constant): + kr = self.k_r0 * ufl.exp(-self.E_kr / (k_B * temperature)) + kd = self.k_d0 * ufl.exp(-self.E_kd / (k_B * temperature)) + if callable(self.gas_pressure): + gas_pressure = self.gas_pressure(t=t) + else: + gas_pressure = self.gas_pressure + + product_of_reactants = self.reactant[0].concentration + for reactant in self.reactant[1:]: + product_of_reactants *= reactant.concentration + + self.value_fenics = kd * gas_pressure - kr * product_of_reactants + + +class SurfaceReactionBC: + """Boundary condition representing a surface reaction + A + B <-> C + where A, B are the reactants and C is the product + the forward reaction rate is K_r = k_r0 * exp(-E_kr / (k_B * T)) + and the backward reaction rate is K_d = k_d0 * exp(-E_kd / (k_B * T)) + The reaction rate is: + K = K_r * C_A * C_B - K_d * P_C + with C_A, C_B the concentration of species A and B, + P_C the partial pressure of species C at the surface. + + The flux of species A entering the surface is K. + In the special case where A=B, then the flux of particle entering the surface is 2*K + + + Args: + reactant (list): list of F.Species objects representing the reactants + gas_pressure (float, callable): the partial pressure of the product species + k_r0 (float): the pre-exponential factor of the forward reaction rate + E_kr (float): the activation energy of the forward reaction rate (eV) + k_d0 (float): the pre-exponential factor of the backward reaction rate + E_kd (float): the activation energy of the backward reaction rate (eV) + subdomain (F.SurfaceSubdomain): the surface subdomain where the reaction occurs + """ + + def __init__( + self, + reactant, + gas_pressure, + k_r0, + E_kr, + k_d0, + E_kd, + subdomain, + ): + self.reactant = reactant + self.gas_pressure = gas_pressure + self.k_r0 = k_r0 + self.E_kr = E_kr + self.k_d0 = k_d0 + self.E_kd = E_kd + self.subdomain = subdomain + + # create the flux boundary condition for each reactant + self.flux_bcs = [ + SurfaceReactionBCpartial( + reactant=self.reactant, + gas_pressure=self.gas_pressure, + k_r0=self.k_r0, + E_kr=self.E_kr, + k_d0=self.k_d0, + E_kd=self.E_kd, + subdomain=self.subdomain, + species=species, + ) + for species in self.reactant + ] + + @property + def time_dependent(self): + return False # no need to update if only using ufl.conditional objects diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 6ff1cd98f..5a7f4d7cd 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -508,10 +508,26 @@ def assign_functions_to_species(self): spe.test_function = sub_test_functions[idx] def define_boundary_conditions(self): + # @jhdark this all_bcs could be a property + # I just don't want to modify self.boundary_conditions + + # create all_bcs which includes all flux bcs from SurfaceReactionBC + all_bcs = self.boundary_conditions.copy() for bc in self.boundary_conditions: + if isinstance(bc, boundary_conditions.SurfaceReactionBC): + all_bcs += bc.flux_bcs + all_bcs.remove(bc) + + for bc in all_bcs: if isinstance(bc.species, str): # if name of species is given then replace with species object bc.species = _species.find_species_from_name(bc.species, self.species) + if isinstance(bc, boundary_conditions.ParticleFluxBC): + bc.create_value_fenics( + mesh=self.mesh.mesh, + temperature=self.temperature_fenics, + t=self.t, + ) super().define_boundary_conditions() @@ -681,6 +697,13 @@ def create_formulation(self): * bc.species.test_function * self.ds(bc.subdomain.id) ) + if isinstance(bc, boundary_conditions.SurfaceReactionBC): + for flux_bc in bc.flux_bcs: + self.formulation -= ( + flux_bc.value_fenics + * flux_bc.species.test_function + * self.ds(flux_bc.subdomain.id) + ) # check if each species is defined in all volumes if not self.settings.transient: @@ -714,8 +737,15 @@ def update_time_dependent_values(self): self.temperature_fenics.interpolate(self.temperature_expr) for bc in self.boundary_conditions: - if bc.temperature_dependent: - bc.update(t=t) + if isinstance( + bc, + ( + boundary_conditions.FixedConcentrationBC, + boundary_conditions.ParticleFluxBC, + ), + ): + if bc.temperature_dependent: + bc.update(t=t) for source in self.sources: if source.temperature_dependent: diff --git a/test/system_tests/test_surface_reactions.py b/test/system_tests/test_surface_reactions.py new file mode 100644 index 000000000..58104b4be --- /dev/null +++ b/test/system_tests/test_surface_reactions.py @@ -0,0 +1,338 @@ +import festim as F +import numpy as np + +import ufl +import dolfinx.fem as fem + + +# TODO add this to the festim package +class FluxFromSurfaceReaction(F.SurfaceFlux): + def __init__(self, reaction: F.SurfaceReactionBC): + super().__init__( + F.Species(), # just a dummy species here + reaction.subdomain, + ) + self.reaction = reaction.flux_bcs[0] + + def compute(self, ds): + self.value = fem.assemble_scalar( + fem.form(self.reaction.value_fenics * ds(self.surface.id)) + ) + self.data.append(self.value) + + +def test_2_isotopes_no_pressure(): + """ + Runs a simple 1D hydrogen transport problem with 2 isotopes + and 3 surface reactions on the right boundary. + + Then checks that the fluxes of the isotopes are consistent with the surface reactions + by computing the flux of H and D (from the gradient) and comparing it to the fluxes + computed from the surface reactions. + + Example: H + D <-> HD, H + H <-> HH, D + D <-> DD + -D grad(c_H) n = 2*kr*c_H*c_H + kr*c_H*c_D + + Also checks the mass balance between the left and right boundary + """ + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 1000)) + my_mat = F.Material(name="mat", D_0=1, E_D=0) + vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat) + left = F.SurfaceSubdomain1D(id=1, x=0) + right = F.SurfaceSubdomain1D(id=2, x=1) + + my_model.subdomains = [vol, left, right] + + H = F.Species("H") + D = F.Species("D") + my_model.species = [H, D] + + my_model.temperature = 500 + + surface_reaction_hd = F.SurfaceReactionBC( + reactant=[H, D], + gas_pressure=0, + k_r0=0.01, + E_kr=0, + k_d0=0, + E_kd=0, + subdomain=right, + ) + + surface_reaction_hh = F.SurfaceReactionBC( + reactant=[H, H], + gas_pressure=0, + k_r0=0.02, + E_kr=0, + k_d0=0, + E_kd=0, + subdomain=right, + ) + + surface_reaction_dd = F.SurfaceReactionBC( + reactant=[D, D], + gas_pressure=0, + k_r0=0.03, + E_kr=0, + k_d0=0, + E_kd=0, + subdomain=right, + ) + + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=2, species=H), + F.DirichletBC(subdomain=left, value=3, species=D), + surface_reaction_hd, + surface_reaction_hh, + surface_reaction_dd, + ] + + H_flux_right = F.SurfaceFlux(H, right) + H_flux_left = F.SurfaceFlux(H, left) + D_flux_right = F.SurfaceFlux(D, right) + D_flux_left = F.SurfaceFlux(D, left) + HD_flux = FluxFromSurfaceReaction(surface_reaction_hd) + HH_flux = FluxFromSurfaceReaction(surface_reaction_hh) + DD_flux = FluxFromSurfaceReaction(surface_reaction_dd) + my_model.exports = [ + H_flux_left, + H_flux_right, + D_flux_left, + D_flux_right, + HD_flux, + HH_flux, + DD_flux, + ] + + my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=5, transient=True) + + my_model.settings.stepsize = 0.1 + + my_model.initialise() + my_model.run() + + # TEST + + assert np.isclose( + np.abs(H_flux_right.data[-1]), + np.abs(H_flux_left.data[-1]), + rtol=0.5e-2, + atol=0.005, + ) + assert np.isclose( + np.abs(D_flux_right.data[-1]), + np.abs(D_flux_left.data[-1]), + rtol=0.5e-2, + atol=0.005, + ) + + # check that H_flux_right == 2*HH_flux + HD_flux + H_flux_from_gradient = -np.array(H_flux_right.data) + H_flux_from_reac = 2 * np.array(HH_flux.data) + np.array(HD_flux.data) + assert np.allclose( + H_flux_from_gradient, + H_flux_from_reac, + rtol=0.5e-2, + atol=0.005, + ) + + # check that D_flux_right == 2*DD_flux + HD_flux + D_flux_from_gradient = -np.array(D_flux_right.data) + D_flux_from_reac = 2 * np.array(DD_flux.data) + np.array(HD_flux.data) + assert np.allclose( + D_flux_from_gradient, + D_flux_from_reac, + rtol=0.5e-2, + atol=0.005, + ) + + +def test_2_isotopes_with_pressure(): + """ + Runs a simple 1D hydrogen transport problem with 2 isotopes + and 3 surface reactions on the right boundary. + + Then checks that the fluxes of the isotopes are consistent with the surface reactions + by computing the flux of H and D (from the gradient) and comparing it to the fluxes + computed from the surface reactions. + + Example: H + D <-> HD, H + H <-> HH, D + D <-> DD + -D grad(c_H) n = 2*kr*c_H*c_H + kr*c_H*c_D - kd*P_H2 + + Also checks the mass balance between the left and right boundary + """ + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 1000)) + my_mat = F.Material(name="mat", D_0=1, E_D=0) + vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat) + left = F.SurfaceSubdomain1D(id=1, x=0) + right = F.SurfaceSubdomain1D(id=2, x=1) + + my_model.subdomains = [vol, left, right] + + H = F.Species("H") + D = F.Species("D") + my_model.species = [H, D] + + my_model.temperature = 500 + + surface_reaction_hd = F.SurfaceReactionBC( + reactant=[H, D], + gas_pressure=0, + k_r0=0.01, + E_kr=0, + k_d0=0, + E_kd=0, + subdomain=right, + ) + + surface_reaction_hh = F.SurfaceReactionBC( + reactant=[H, H], + gas_pressure=2, + k_r0=0.02, + E_kr=0, + k_d0=0.1, + E_kd=0, + subdomain=right, + ) + + surface_reaction_dd = F.SurfaceReactionBC( + reactant=[D, D], + gas_pressure=0, + k_r0=0.03, + E_kr=0, + k_d0=0, + E_kd=0, + subdomain=right, + ) + + my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=2, species=H), + F.DirichletBC(subdomain=left, value=3, species=D), + surface_reaction_hd, + surface_reaction_hh, + surface_reaction_dd, + ] + + H_flux_right = F.SurfaceFlux(H, right) + H_flux_left = F.SurfaceFlux(H, left) + D_flux_right = F.SurfaceFlux(D, right) + D_flux_left = F.SurfaceFlux(D, left) + HD_flux = FluxFromSurfaceReaction(surface_reaction_hd) + HH_flux = FluxFromSurfaceReaction(surface_reaction_hh) + DD_flux = FluxFromSurfaceReaction(surface_reaction_dd) + my_model.exports = [ + H_flux_left, + H_flux_right, + D_flux_left, + D_flux_right, + HD_flux, + HH_flux, + DD_flux, + ] + + my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=5, transient=True) + + my_model.settings.stepsize = 0.1 + + my_model.initialise() + my_model.run() + + # TEST + + assert np.isclose( + np.abs(H_flux_right.data[-1]), + np.abs(H_flux_left.data[-1]), + rtol=0.5e-2, + atol=0.005, + ) + assert np.isclose( + np.abs(D_flux_right.data[-1]), + np.abs(D_flux_left.data[-1]), + rtol=0.5e-2, + atol=0.005, + ) + + # check that H_flux_right == 2*HH_flux + HD_flux + H_flux_from_gradient = -np.array(H_flux_right.data) + H_flux_from_reac = 2 * np.array(HH_flux.data) + np.array(HD_flux.data) + assert np.allclose( + H_flux_from_gradient, + H_flux_from_reac, + rtol=0.5e-2, + atol=0.005, + ) + + # check that D_flux_right == 2*DD_flux + HD_flux + D_flux_from_gradient = -np.array(D_flux_right.data) + D_flux_from_reac = 2 * np.array(DD_flux.data) + np.array(HD_flux.data) + assert np.allclose( + D_flux_from_gradient, + D_flux_from_reac, + rtol=0.5e-2, + atol=0.005, + ) + + +def test_pressure_varies_in_time(): + """ + Runs a problem with a surface reaction and a time-dependent pressure + on the right boundary. + + Then checks that the flux is consistent with the surface reaction + """ + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 1000)) + my_mat = F.Material(name="mat", D_0=1, E_D=0) + vol = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat) + left = F.SurfaceSubdomain1D(id=1, x=0) + right = F.SurfaceSubdomain1D(id=2, x=1) + + my_model.subdomains = [vol, left, right] + + H = F.Species("H") + my_model.species = [H] + + my_model.temperature = 500 + + t_pressure = 2 + pressure = 2 + k_d = 2 + + surface_reaction_hh = F.SurfaceReactionBC( + reactant=[H, H], + gas_pressure=lambda t: ufl.conditional(ufl.gt(t, t_pressure), pressure, 0), + k_r0=0, + E_kr=0, + k_d0=k_d, + E_kd=0, + subdomain=right, + ) + + my_model.boundary_conditions = [surface_reaction_hh] + + H_flux_right = F.SurfaceFlux(H, right) + my_model.exports = [H_flux_right] + + my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=5, transient=True) + + my_model.settings.stepsize = 0.1 + + my_model.initialise() + my_model.run() + + flux_as_array = np.array(H_flux_right.data) + time_as_array = np.array(H_flux_right.t) + + expected_flux_before_pressure = 0 + computed_flux_before_pressure = flux_as_array[time_as_array <= t_pressure] + assert np.allclose(computed_flux_before_pressure, expected_flux_before_pressure) + + expected_flux_after_pressure = -2 * k_d * pressure + computed_flux_after_pressure = flux_as_array[time_as_array > t_pressure] + print(computed_flux_after_pressure) + print(expected_flux_after_pressure) + assert np.allclose( + computed_flux_after_pressure, expected_flux_after_pressure, rtol=1e-2 + ) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 285c6dcd2..ee64ba876 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -1105,7 +1105,7 @@ def test_update_time_dependent_values_flux(bc_value, expected_values): my_model.define_meshtags_and_measures() my_model.assign_functions_to_species() my_model.define_temperature() - my_model.create_flux_values_fenics() + my_model.define_boundary_conditions() for i in range(3): # RUN @@ -1159,7 +1159,7 @@ def test_update_fluxes_with_time_dependent_temperature( my_model.define_function_spaces() my_model.assign_functions_to_species() my_model.define_meshtags_and_measures() - my_model.create_flux_values_fenics() + my_model.define_boundary_conditions() for i in range(3): # RUN