From 1e686fe574554cbcc7bf2d554a2dc98609e5dfd7 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 10 Apr 2024 11:48:50 -0400 Subject: [PATCH 01/19] first skeleton of class --- .../boundary_conditions/surface_reaction.py | 37 +++++++++++++++++++ 1 file changed, 37 insertions(+) create mode 100644 festim/boundary_conditions/surface_reaction.py diff --git a/festim/boundary_conditions/surface_reaction.py b/festim/boundary_conditions/surface_reaction.py new file mode 100644 index 000000000..0472b9aac --- /dev/null +++ b/festim/boundary_conditions/surface_reaction.py @@ -0,0 +1,37 @@ +import festim as F +import fenics as f + + +class SurfaceReactionBC(F.ParticleFluxBC): + 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 + + value = self.reaction_rate() + + super().__init__(subdomain, value, species) + + def reaction_rate(self, T): + kr = self.k_r0 * f.exp(-self.E_kr / (F.k_B * T)) + kd = self.k_d0 * f.exp(-self.E_kd / (F.k_B * T)) + if callable( + self.gas_pressure + ): # need to check if gas_pressure is a function of time only + return lambda t: kd * self.gas_pressure(t) - kr * self.reactant + else: + return kd * self.gas_pressure - kr * self.reactant From 2e84a71c4932cc850f77ddb893b50499ee098640 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 10 Apr 2024 14:12:53 -0400 Subject: [PATCH 02/19] ufl --- festim/boundary_conditions/surface_reaction.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/festim/boundary_conditions/surface_reaction.py b/festim/boundary_conditions/surface_reaction.py index 0472b9aac..109da65ef 100644 --- a/festim/boundary_conditions/surface_reaction.py +++ b/festim/boundary_conditions/surface_reaction.py @@ -1,5 +1,5 @@ import festim as F -import fenics as f +import ufl class SurfaceReactionBC(F.ParticleFluxBC): @@ -27,8 +27,8 @@ def __init__( super().__init__(subdomain, value, species) def reaction_rate(self, T): - kr = self.k_r0 * f.exp(-self.E_kr / (F.k_B * T)) - kd = self.k_d0 * f.exp(-self.E_kd / (F.k_B * T)) + kr = self.k_r0 * ufl.exp(-self.E_kr / (F.k_B * T)) + kd = self.k_d0 * ufl.exp(-self.E_kd / (F.k_B * T)) if callable( self.gas_pressure ): # need to check if gas_pressure is a function of time only From 361d7761162c232b952db6305a038e2025422a92 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 12 Apr 2024 10:58:22 -0400 Subject: [PATCH 03/19] initial draft --- festim/__init__.py | 1 + .../boundary_conditions/surface_reaction.py | 34 +++++++++++++------ festim/hydrogen_transport_problem.py | 7 ++++ 3 files changed, 31 insertions(+), 11 deletions(-) diff --git a/festim/__init__.py b/festim/__init__.py index a294d9752..8bb118248 100644 --- a/festim/__init__.py +++ b/festim/__init__.py @@ -27,6 +27,7 @@ from .boundary_conditions.sieverts_bc import SievertsBC from .boundary_conditions.henrys_bc import HenrysBC from .boundary_conditions.flux_bc import FluxBCBase, ParticleFluxBC, HeatFluxBC +from .boundary_conditions.surface_reaction import SurfaceReactionBC from .material import Material diff --git a/festim/boundary_conditions/surface_reaction.py b/festim/boundary_conditions/surface_reaction.py index 109da65ef..ad7fdfdf1 100644 --- a/festim/boundary_conditions/surface_reaction.py +++ b/festim/boundary_conditions/surface_reaction.py @@ -2,6 +2,7 @@ import ufl +# TODO instead of inheriting from ParticleFluxBC, it should be composed of several ParticleFluxBCs class SurfaceReactionBC(F.ParticleFluxBC): def __init__( self, @@ -21,17 +22,28 @@ def __init__( self.E_kr = E_kr self.k_d0 = k_d0 self.E_kd = E_kd + name_to_species = {f"c{i}": species for i, species in enumerate(self.reactant)} - value = self.reaction_rate() + value = self.create_value() - super().__init__(subdomain, value, species) + super().__init__( + subdomain, value, species, species_dependent_value=name_to_species + ) - def reaction_rate(self, T): - kr = self.k_r0 * ufl.exp(-self.E_kr / (F.k_B * T)) - kd = self.k_d0 * ufl.exp(-self.E_kd / (F.k_B * T)) - if callable( - self.gas_pressure - ): # need to check if gas_pressure is a function of time only - return lambda t: kd * self.gas_pressure(t) - kr * self.reactant - else: - return kd * self.gas_pressure - kr * self.reactant + def create_value(self): + kwargs = {f"c{i}": None for i, _ in enumerate(self.reactant)} + kwargs["T"] = None + if callable(self.gas_pressure): + kwargs["t"] = None + + def value(**kwargs): + T = kwargs["T"] + kr = self.k_r0 * ufl.exp(-self.E_kr / (F.k_B * T)) + kd = self.k_d0 * ufl.exp(-self.E_kd / (F.k_B * T)) + if callable(self.gas_pressure): + gas_pressure = self.gas_pressure(t=kwargs["t"]) + else: + gas_pressure = self.gas_pressure + return kd * gas_pressure - kr * self.reactant + + return value diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index f14826feb..bbf28d07b 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -701,6 +701,13 @@ def create_formulation(self): * bc.species.test_function * self.ds(bc.subdomain.id) ) + if isinstance(bc, F.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: From acff010beb6da27fcd07855cd4bebb64735d812b Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 12 Apr 2024 15:43:30 -0400 Subject: [PATCH 04/19] example + working implementation --- .../boundary_conditions/surface_reaction.py | 81 +++++++++++++------ festim/hydrogen_transport_problem.py | 27 ++++--- surface_reaction_example.py | 52 ++++++++++++ 3 files changed, 123 insertions(+), 37 deletions(-) create mode 100644 surface_reaction_example.py diff --git a/festim/boundary_conditions/surface_reaction.py b/festim/boundary_conditions/surface_reaction.py index ad7fdfdf1..4a90d150a 100644 --- a/festim/boundary_conditions/surface_reaction.py +++ b/festim/boundary_conditions/surface_reaction.py @@ -1,9 +1,9 @@ import festim as F +from dolfinx import fem import ufl -# TODO instead of inheriting from ParticleFluxBC, it should be composed of several ParticleFluxBCs -class SurfaceReactionBC(F.ParticleFluxBC): +class SurfaceReactionFlux(F.ParticleFluxBC): def __init__( self, reactant, @@ -15,6 +15,42 @@ def __init__( 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 reaction_rate(self, temperature, t: fem.Constant): + kr = self.k_r0 * ufl.exp(-self.E_kr / (F.k_B * temperature)) + kd = self.k_d0 * ufl.exp(-self.E_kd / (F.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 + + return kd * gas_pressure - kr * product_of_reactants + + def create_value_fenics(self, mesh, temperature, t: fem.Constant): + self.value_fenics = self.reaction_rate(temperature, t) + + +class SurfaceReactionBC: + def __init__( + self, + reactant, + gas_pressure, + k_r0, + E_kr, + k_d0, + E_kd, + subdomain, + ): self.reactant = reactant self.gas_pressure = gas_pressure @@ -22,28 +58,25 @@ def __init__( self.E_kr = E_kr self.k_d0 = k_d0 self.E_kd = E_kd - name_to_species = {f"c{i}": species for i, species in enumerate(self.reactant)} + self.subdomain = subdomain - value = self.create_value() + self.flux_bcs = self.create_flux_bcs() - super().__init__( - subdomain, value, species, species_dependent_value=name_to_species - ) + @property + def time_dependent(self): + return False # no need to update if only using ufl.conditional objects - def create_value(self): - kwargs = {f"c{i}": None for i, _ in enumerate(self.reactant)} - kwargs["T"] = None - if callable(self.gas_pressure): - kwargs["t"] = None - - def value(**kwargs): - T = kwargs["T"] - kr = self.k_r0 * ufl.exp(-self.E_kr / (F.k_B * T)) - kd = self.k_d0 * ufl.exp(-self.E_kd / (F.k_B * T)) - if callable(self.gas_pressure): - gas_pressure = self.gas_pressure(t=kwargs["t"]) - else: - gas_pressure = self.gas_pressure - return kd * gas_pressure - kr * self.reactant - - return value + def create_flux_bcs(self): + return [ + SurfaceReactionFlux( + 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 + ] diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index bbf28d07b..94b93938a 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -248,7 +248,6 @@ def initialise(self): self.define_temperature() self.define_boundary_conditions() self.create_source_values_fenics() - self.create_flux_values_fenics() self.create_initial_conditions() self.create_formulation() self.create_solver() @@ -514,12 +513,26 @@ def define_meshtags_and_measures(self): def define_boundary_conditions(self): """Defines the dirichlet boundary conditions of the model""" for bc in self.boundary_conditions: + if isinstance(bc, F.SurfaceReactionBC): + for flux_bc in bc.flux_bcs: + flux_bc.create_value_fenics( + mesh=self.mesh.mesh, + temperature=self.temperature_fenics, + t=self.t, + ) + continue if isinstance(bc.species, str): # if name of species is given then replace with species object bc.species = F.find_species_from_name(bc.species, self.species) if isinstance(bc, F.DirichletBC): form = self.create_dirichletbc_form(bc) self.bc_forms.append(form) + if isinstance(bc, F.ParticleFluxBC): + bc.create_value_fenics( + mesh=self.mesh.mesh, + temperature=self.temperature_fenics, + t=self.t, + ) def create_dirichletbc_form(self, bc): """Creates a dirichlet boundary condition form @@ -590,18 +603,6 @@ def create_source_values_fenics(self): t=self.t, ) - def create_flux_values_fenics(self): - """For each particle flux create the value_fenics""" - for bc in self.boundary_conditions: - # create value_fenics for all F.ParticleFluxBC objects - if isinstance(bc, F.ParticleFluxBC): - - bc.create_value_fenics( - mesh=self.mesh.mesh, - temperature=self.temperature_fenics, - t=self.t, - ) - def create_initial_conditions(self): """For each initial condition, create the value_fenics and assign it to the previous solution of the condition's species""" diff --git a/surface_reaction_example.py b/surface_reaction_example.py new file mode 100644 index 000000000..14924b13e --- /dev/null +++ b/surface_reaction_example.py @@ -0,0 +1,52 @@ +import festim as F +import numpy as np +import ufl + + +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 = 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_manual = F.ParticleFluxBC( + subdomain=right, + value=lambda c1, T: -2 * 0.01 * ufl.exp(-0 / F.k_B / T) * c1**2, + species=H, + species_dependent_value={"c1": H}, +) + +my_model.boundary_conditions = [ + F.DirichletBC(subdomain=left, value=5, species=H), + F.DirichletBC(subdomain=left, value=5, species=D), + surface_reaction, + # surface_reaction_manual, +] + +my_model.exports = [F.XDMFExport("test.xdmf", H)] + +my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=10, transient=True) + +my_model.settings.stepsize = 0.1 + +my_model.initialise() +my_model.run() From 9544d7417e9028e96628a784015bbb77c397ae7a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 12 Apr 2024 16:00:53 -0400 Subject: [PATCH 05/19] refactoring --- .../boundary_conditions/surface_reaction.py | 24 +++++++------------ festim/hydrogen_transport_problem.py | 12 ++++------ 2 files changed, 14 insertions(+), 22 deletions(-) diff --git a/festim/boundary_conditions/surface_reaction.py b/festim/boundary_conditions/surface_reaction.py index 4a90d150a..8506cf8dd 100644 --- a/festim/boundary_conditions/surface_reaction.py +++ b/festim/boundary_conditions/surface_reaction.py @@ -3,7 +3,7 @@ import ufl -class SurfaceReactionFlux(F.ParticleFluxBC): +class SurfaceReactionBCpartial(F.ParticleFluxBC): def __init__( self, reactant, @@ -23,7 +23,7 @@ def __init__( self.E_kd = E_kd super().__init__(subdomain=subdomain, value=None, species=species) - def reaction_rate(self, temperature, t: fem.Constant): + def create_value_fenics(self, mesh, temperature, t: fem.Constant): kr = self.k_r0 * ufl.exp(-self.E_kr / (F.k_B * temperature)) kd = self.k_d0 * ufl.exp(-self.E_kd / (F.k_B * temperature)) if callable(self.gas_pressure): @@ -34,10 +34,7 @@ def reaction_rate(self, temperature, t: fem.Constant): for reactant in self.reactant[1:]: product_of_reactants *= reactant.concentration - return kd * gas_pressure - kr * product_of_reactants - - def create_value_fenics(self, mesh, temperature, t: fem.Constant): - self.value_fenics = self.reaction_rate(temperature, t) + self.value_fenics = kd * gas_pressure - kr * product_of_reactants class SurfaceReactionBC: @@ -60,15 +57,8 @@ def __init__( self.E_kd = E_kd self.subdomain = subdomain - self.flux_bcs = self.create_flux_bcs() - - @property - def time_dependent(self): - return False # no need to update if only using ufl.conditional objects - - def create_flux_bcs(self): - return [ - SurfaceReactionFlux( + self.flux_bcs = [ + SurfaceReactionBCpartial( reactant=self.reactant, gas_pressure=self.gas_pressure, k_r0=self.k_r0, @@ -80,3 +70,7 @@ def create_flux_bcs(self): ) 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/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 94b93938a..0252b4506 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -512,15 +512,13 @@ def define_meshtags_and_measures(self): def define_boundary_conditions(self): """Defines the dirichlet boundary conditions of the model""" + all_bcs = self.boundary_conditions.copy() for bc in self.boundary_conditions: if isinstance(bc, F.SurfaceReactionBC): - for flux_bc in bc.flux_bcs: - flux_bc.create_value_fenics( - mesh=self.mesh.mesh, - temperature=self.temperature_fenics, - t=self.t, - ) - continue + 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 = F.find_species_from_name(bc.species, self.species) From ef70734c67e33b5f205485013af2e1d15931dc79 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 12 Apr 2024 16:31:55 -0400 Subject: [PATCH 06/19] documentation --- .../boundary_conditions/surface_reaction.py | 51 +++++++++++++++++++ festim/hydrogen_transport_problem.py | 6 ++- 2 files changed, 56 insertions(+), 1 deletion(-) diff --git a/festim/boundary_conditions/surface_reaction.py b/festim/boundary_conditions/surface_reaction.py index 8506cf8dd..4b5faef81 100644 --- a/festim/boundary_conditions/surface_reaction.py +++ b/festim/boundary_conditions/surface_reaction.py @@ -4,6 +4,31 @@ class SurfaceReactionBCpartial(F.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, @@ -30,6 +55,7 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant): 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 @@ -38,6 +64,30 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant): 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, @@ -57,6 +107,7 @@ def __init__( self.E_kd = E_kd self.subdomain = subdomain + # create the flux boundary condition for each reactant self.flux_bcs = [ SurfaceReactionBCpartial( reactant=self.reactant, diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 0252b4506..e860d1573 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -511,7 +511,11 @@ def define_meshtags_and_measures(self): ) def define_boundary_conditions(self): - """Defines the dirichlet boundary conditions of the model""" + """Create forms for DirichletBC and value_fenics for ParticleFluxBC""" + # @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, F.SurfaceReactionBC): From b1f2332b3f74d035f1b7bfccacd26aca0d4bf48d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 18 Apr 2024 08:44:46 -0400 Subject: [PATCH 07/19] fixed tests --- test/test_h_transport_problem.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 23f3ebf69..25e9f9e2b 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -964,7 +964,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 @@ -1018,7 +1018,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 @@ -1031,8 +1031,8 @@ def test_update_fluxes_with_time_dependent_temperature( assert np.isclose(computed_value, expected_values[i]) -def test_create_source_values_fenics_multispecies(): - """Test that the create_flux_values_fenics method correctly sets the value_fenics +def test_define_boundary_conditions_multispecies(): + """Test that the define_boundary_conditions method correctly sets the value_fenics attribute in a multispecies case""" # BUILD my_vol = F.VolumeSubdomain1D(id=1, borders=[0, 4], material=dummy_mat) @@ -1056,7 +1056,7 @@ def test_create_source_values_fenics_multispecies(): my_model.define_temperature() # RUN - my_model.create_flux_values_fenics() + my_model.define_boundary_conditions() # TEST assert np.isclose(my_model.boundary_conditions[0].value_fenics.value, 5) From 4a55c640c7e7f1345d2e7d335555cbc68fe22be9 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 May 2024 13:54:26 +0200 Subject: [PATCH 08/19] added HH and DD reactions to surface --- surface_reaction_example.py | 30 +++++++++++++++++++++--------- 1 file changed, 21 insertions(+), 9 deletions(-) diff --git a/surface_reaction_example.py b/surface_reaction_example.py index 14924b13e..1a26aa729 100644 --- a/surface_reaction_example.py +++ b/surface_reaction_example.py @@ -1,7 +1,5 @@ import festim as F import numpy as np -import ufl - my_model = F.HydrogenTransportProblem() my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, 1000)) @@ -18,7 +16,7 @@ my_model.temperature = 500 -surface_reaction = F.SurfaceReactionBC( +surface_reaction_hd = F.SurfaceReactionBC( reactant=[H, D], gas_pressure=0, k_r0=0.01, @@ -28,18 +26,32 @@ subdomain=right, ) -surface_reaction_manual = F.ParticleFluxBC( +surface_reaction_hh = F.SurfaceReactionBC( + reactant=[H, H], + gas_pressure=0, + k_r0=0.01, + 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.01, + E_kr=0, + k_d0=0, + E_kd=0, subdomain=right, - value=lambda c1, T: -2 * 0.01 * ufl.exp(-0 / F.k_B / T) * c1**2, - species=H, - species_dependent_value={"c1": H}, ) my_model.boundary_conditions = [ F.DirichletBC(subdomain=left, value=5, species=H), F.DirichletBC(subdomain=left, value=5, species=D), - surface_reaction, - # surface_reaction_manual, + surface_reaction_hd, + surface_reaction_hh, + surface_reaction_dd, ] my_model.exports = [F.XDMFExport("test.xdmf", H)] From 575323f8a295f49aa4ee4dc80cb4ba76b08ff095 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 May 2024 14:06:16 +0200 Subject: [PATCH 09/19] added derived quantities --- surface_reaction_example.py | 31 +++++++++++++++++++++++++++++-- 1 file changed, 29 insertions(+), 2 deletions(-) diff --git a/surface_reaction_example.py b/surface_reaction_example.py index 1a26aa729..a8b031a9c 100644 --- a/surface_reaction_example.py +++ b/surface_reaction_example.py @@ -48,13 +48,23 @@ my_model.boundary_conditions = [ F.DirichletBC(subdomain=left, value=5, species=H), - F.DirichletBC(subdomain=left, value=5, species=D), + F.DirichletBC(subdomain=left, value=2, species=D), surface_reaction_hd, surface_reaction_hh, surface_reaction_dd, ] -my_model.exports = [F.XDMFExport("test.xdmf", H)] +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) +my_model.exports = [ + F.XDMFExport("test.xdmf", H), + H_flux_left, + H_flux_right, + D_flux_left, + D_flux_right, +] my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=10, transient=True) @@ -62,3 +72,20 @@ 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.show() From 769558709dfa91ad6f33c038c454c5c699cab3be Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 May 2024 15:58:49 +0200 Subject: [PATCH 10/19] added system test --- test/system_tests/test_surface_reactions.py | 120 ++++++++++++++++++++ 1 file changed, 120 insertions(+) create mode 100644 test/system_tests/test_surface_reactions.py diff --git a/test/system_tests/test_surface_reactions.py b/test/system_tests/test_surface_reactions.py new file mode 100644 index 000000000..39ea78b46 --- /dev/null +++ b/test/system_tests/test_surface_reactions.py @@ -0,0 +1,120 @@ +import festim as F +import numpy as np + +import dolfinx.fem as fem + + +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(): + 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 + + # check that H_flux_right == 2*HH_flux + HD_flux + H_flux_from_gradient = np.abs(H_flux_right.data) + H_flux_from_reac = 2 * np.abs(HH_flux.data) + np.abs(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.abs(D_flux_right.data) + D_flux_from_reac = 2 * np.abs(DD_flux.data) + np.abs(HD_flux.data) + assert np.allclose( + D_flux_from_gradient, + D_flux_from_reac, + rtol=0.5e-2, + atol=0.005, + ) From 06dbc340061a4b663921fd5e78ffa9dc922c65b7 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 May 2024 16:02:08 +0200 Subject: [PATCH 11/19] added plots --- surface_reaction_example.py | 69 +++++++++++++++++++++++++++++++++++-- 1 file changed, 67 insertions(+), 2 deletions(-) diff --git a/surface_reaction_example.py b/surface_reaction_example.py index a8b031a9c..b74821719 100644 --- a/surface_reaction_example.py +++ b/surface_reaction_example.py @@ -1,6 +1,24 @@ import festim as F import numpy as np +import dolfinx.fem as fem + + +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) @@ -47,7 +65,7 @@ ) my_model.boundary_conditions = [ - F.DirichletBC(subdomain=left, value=5, species=H), + F.DirichletBC(subdomain=left, value=2, species=H), F.DirichletBC(subdomain=left, value=2, species=D), surface_reaction_hd, surface_reaction_hh, @@ -58,15 +76,21 @@ 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=10, transient=True) +my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=5, transient=True) my_model.settings.stepsize = 0.1 @@ -88,4 +112,45 @@ labels=["H_out", "D_out"], ) plt.legend() + +plt.figure() +plt.stackplot( + HD_flux.t, + np.abs(HD_flux.data), + np.abs(HH_flux.data), + np.abs(DD_flux.data), + labels=["HD", "HH", "DD"], +) +plt.legend() + + +plt.figure() +plt.plot(H_flux_right.t, np.abs(H_flux_right.data)) +plt.plot( + H_flux_right.t, 2 * np.abs(HH_flux.data) + np.abs(HD_flux.data), linestyle="--" +) + +plt.plot(D_flux_right.t, np.abs(D_flux_right.data)) +plt.plot( + D_flux_right.t, 2 * np.abs(DD_flux.data) + np.abs(HD_flux.data), linestyle="--" +) + +# check that H_flux_right == 2*HH_flux + HD_flux +H_flux_from_gradient = np.abs(H_flux_right.data) +H_flux_from_reac = 2 * np.abs(HH_flux.data) + np.abs(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.abs(D_flux_right.data) +D_flux_from_reac = 2 * np.abs(DD_flux.data) + np.abs(HD_flux.data) +assert np.allclose( + D_flux_from_gradient, + D_flux_from_reac, + rtol=0.5e-2, + atol=0.005, +) plt.show() From 2612b3567d93e83ea880f0ddbcad352f2bb018c1 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 May 2024 16:17:15 +0200 Subject: [PATCH 12/19] add two system tests --- test/system_tests/test_surface_reactions.py | 164 +++++++++++++++++++- 1 file changed, 159 insertions(+), 5 deletions(-) diff --git a/test/system_tests/test_surface_reactions.py b/test/system_tests/test_surface_reactions.py index 39ea78b46..7dd5ca285 100644 --- a/test/system_tests/test_surface_reactions.py +++ b/test/system_tests/test_surface_reactions.py @@ -4,6 +4,7 @@ import dolfinx.fem as fem +# TODO add this to the festim package class FluxFromSurfaceReaction(F.SurfaceFlux): def __init__(self, reaction: F.SurfaceReactionBC): super().__init__( @@ -19,7 +20,20 @@ def compute(self, ds): self.data.append(self.value) -def test(): +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) @@ -99,9 +113,149 @@ def test(): # 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.abs(H_flux_right.data) - H_flux_from_reac = 2 * np.abs(HH_flux.data) + np.abs(HD_flux.data) + 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, @@ -110,8 +264,8 @@ def test(): ) # check that D_flux_right == 2*DD_flux + HD_flux - D_flux_from_gradient = np.abs(D_flux_right.data) - D_flux_from_reac = 2 * np.abs(DD_flux.data) + np.abs(HD_flux.data) + 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, From 170578681e81a5352a581450229544e3ee404efa Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 May 2024 16:25:22 +0200 Subject: [PATCH 13/19] labels to axes and legend --- surface_reaction_example.py | 40 ++++++++++++++++++++++++------------- 1 file changed, 26 insertions(+), 14 deletions(-) diff --git a/surface_reaction_example.py b/surface_reaction_example.py index b74821719..077de3ac3 100644 --- a/surface_reaction_example.py +++ b/surface_reaction_example.py @@ -47,7 +47,7 @@ def compute(self, ds): surface_reaction_hh = F.SurfaceReactionBC( reactant=[H, H], gas_pressure=0, - k_r0=0.01, + k_r0=0.02, E_kr=0, k_d0=0, E_kd=0, @@ -112,32 +112,45 @@ def compute(self, ds): 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(HD_flux.data), np.abs(HH_flux.data), + np.abs(HD_flux.data), np.abs(DD_flux.data), - labels=["HD", "HH", "DD"], + labels=["HH", "HD", "DD"], ) -plt.legend() +plt.legend(reverse=True) +plt.xlabel("Time (s)") +plt.ylabel("Flux (molecule/m^2/s)") plt.figure() -plt.plot(H_flux_right.t, np.abs(H_flux_right.data)) +plt.plot(H_flux_right.t, -np.array(H_flux_right.data), label="from gradient (H)") plt.plot( - H_flux_right.t, 2 * np.abs(HH_flux.data) + np.abs(HD_flux.data), linestyle="--" + 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.abs(D_flux_right.data)) +plt.plot(D_flux_right.t, -np.array(D_flux_right.data), label="from gradient (D)") plt.plot( - D_flux_right.t, 2 * np.abs(DD_flux.data) + np.abs(HD_flux.data), linestyle="--" + 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.abs(H_flux_right.data) -H_flux_from_reac = 2 * np.abs(HH_flux.data) + np.abs(HD_flux.data) +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, @@ -145,12 +158,11 @@ def compute(self, ds): atol=0.005, ) # check that D_flux_right == 2*DD_flux + HD_flux -D_flux_from_gradient = np.abs(D_flux_right.data) -D_flux_from_reac = 2 * np.abs(DD_flux.data) + np.abs(HD_flux.data) +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, ) -plt.show() From 5d71100734130d2de8d6af77f0a9235998ff37e9 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 May 2024 16:25:45 +0200 Subject: [PATCH 14/19] moved to examples folder --- .../surface_reaction_example.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename surface_reaction_example.py => examples/surface_reaction_example.py (100%) diff --git a/surface_reaction_example.py b/examples/surface_reaction_example.py similarity index 100% rename from surface_reaction_example.py rename to examples/surface_reaction_example.py From 8d3c4da9e1f1a0dc7b08b7523e9ade9c2da6fdad Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 May 2024 16:49:04 +0200 Subject: [PATCH 15/19] time dependent pressure --- examples/surface_reaction_example.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/surface_reaction_example.py b/examples/surface_reaction_example.py index 077de3ac3..2caef8f03 100644 --- a/examples/surface_reaction_example.py +++ b/examples/surface_reaction_example.py @@ -2,6 +2,7 @@ import numpy as np import dolfinx.fem as fem +import ufl class FluxFromSurfaceReaction(F.SurfaceFlux): @@ -46,10 +47,10 @@ def compute(self, ds): surface_reaction_hh = F.SurfaceReactionBC( reactant=[H, H], - gas_pressure=0, + gas_pressure=lambda t: ufl.conditional(ufl.gt(t, 1), 2, 0), k_r0=0.02, E_kr=0, - k_d0=0, + k_d0=0.03, E_kd=0, subdomain=right, ) From 00815e785f0eb82429ee2df7a39bf4c1ca5f7067 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 May 2024 10:45:20 +0200 Subject: [PATCH 16/19] added additional test --- test/system_tests/test_surface_reactions.py | 62 +++++++++++++++++++++ 1 file changed, 62 insertions(+) diff --git a/test/system_tests/test_surface_reactions.py b/test/system_tests/test_surface_reactions.py index 7dd5ca285..03dbc64b4 100644 --- a/test/system_tests/test_surface_reactions.py +++ b/test/system_tests/test_surface_reactions.py @@ -1,6 +1,7 @@ import festim as F import numpy as np +import ufl import dolfinx.fem as fem @@ -272,3 +273,64 @@ def test_2_isotopes_with_pressure(): 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] + assert np.allclose( + computed_flux_after_pressure, expected_flux_after_pressure, rtol=1e-2 + ) From ad3b10b1a0c736ea3cd1e77835ef70de82208a10 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 15 Oct 2024 11:20:51 +0800 Subject: [PATCH 17/19] Check if BC is a BC.... --- festim/hydrogen_transport_problem.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 274d2a5d0..03dc1f36e 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -682,8 +682,9 @@ 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, (F.FixedConcentrationBC, F.ParticleFluxBC)): + if bc.temperature_dependent: + bc.update(t=t) for source in self.sources: if source.temperature_dependent: From 7b5899cb35014e285b0754ab7421111201a21cac Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 1 Nov 2024 14:21:30 -0400 Subject: [PATCH 18/19] fix tests --- src/festim/boundary_conditions/__init__.py | 9 ++++- src/festim/hydrogen_transport_problem.py | 43 +++++++++++++++++---- test/system_tests/test_surface_reactions.py | 2 + test/test_h_transport_problem.py | 2 +- 4 files changed, 46 insertions(+), 10 deletions(-) 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/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 4d426b73c..5a7f4d7cd 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -508,13 +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.FixedConcentrationBC): - 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.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() @@ -684,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: @@ -717,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 index 03dbc64b4..58104b4be 100644 --- a/test/system_tests/test_surface_reactions.py +++ b/test/system_tests/test_surface_reactions.py @@ -331,6 +331,8 @@ def test_pressure_varies_in_time(): 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 9a0152590..ee64ba876 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -1197,7 +1197,7 @@ def test_create_flux_values_fenics_multispecies(): my_model.define_temperature() # RUN - my_model.define_boundary_conditions() + my_model.create_flux_values_fenics() # TEST assert np.isclose(my_model.boundary_conditions[0].value_fenics.value, 5) From 9203ac7fe1a5039545df1ff287d5497252b13e48 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 1 Nov 2024 14:35:58 -0400 Subject: [PATCH 19/19] ruff --- src/festim/boundary_conditions/surface_reaction.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/festim/boundary_conditions/surface_reaction.py b/src/festim/boundary_conditions/surface_reaction.py index d4a6cf502..0bde33e6c 100644 --- a/src/festim/boundary_conditions/surface_reaction.py +++ b/src/festim/boundary_conditions/surface_reaction.py @@ -99,7 +99,6 @@ def __init__( E_kd, subdomain, ): - self.reactant = reactant self.gas_pressure = gas_pressure self.k_r0 = k_r0