diff --git a/src/festim/boundary_conditions/flux_bc.py b/src/festim/boundary_conditions/flux_bc.py index c840e1314..5ae48c777 100644 --- a/src/festim/boundary_conditions/flux_bc.py +++ b/src/festim/boundary_conditions/flux_bc.py @@ -187,6 +187,7 @@ def __init__(self, subdomain, value, species, species_dependent_value={}): super().__init__(subdomain=subdomain, value=value) self.species = species self.species_dependent_value = species_dependent_value + self._volume_subdomain = None def create_value_fenics(self, mesh, temperature, t: fem.Constant): """Creates the value of the boundary condition as a fenics object and sets it to @@ -227,7 +228,12 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant): kwargs["T"] = temperature for name, species in self.species_dependent_value.items(): - kwargs[name] = species.concentration + if species.concentration: + kwargs[name] = species.concentration + else: # probably in discontinuous case + kwargs[name] = species.subdomain_to_solution[ + self._volume_subdomain + ] self.value_fenics = self.value(**kwargs) diff --git a/src/festim/boundary_conditions/surface_reaction.py b/src/festim/boundary_conditions/surface_reaction.py index bb54de609..ca8efe92e 100644 --- a/src/festim/boundary_conditions/surface_reaction.py +++ b/src/festim/boundary_conditions/surface_reaction.py @@ -58,9 +58,22 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant): else: gas_pressure = self.gas_pressure - product_of_reactants = self.reactant[0].concentration + # initialise with the concentration of the first reactant + if self.reactant[0].concentration: + product_of_reactants = self.reactant[0].concentration + else: # probably used in HydrogenTransportProblemDiscontinuous + product_of_reactants = self.reactant[0].subdomain_to_solution[ + self._volume_subdomain + ] + + # then compute the product of the concentrations of the other reactants for reactant in self.reactant[1:]: - product_of_reactants *= reactant.concentration + if reactant.concentration: + product_of_reactants *= reactant.concentration + else: # probably used in HydrogenTransportProblemDiscontinuous + product_of_reactants *= reactant.subdomain_to_solution[ + self._volume_subdomain + ] self.value_fenics = kd * gas_pressure - kr * product_of_reactants diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 5497c8aa4..e783f1a39 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -296,6 +296,17 @@ def volume_meshtags(self, value): else: raise TypeError("value must be of type dolfinx.mesh.MeshTags") + @property + def _unpacked_bcs(self): + """Returns all boundary conditions, including fluxes from surface reactions""" + all_boundary_conditions = [] + for bc in self.boundary_conditions: + if isinstance(bc, boundary_conditions.SurfaceReactionBC): + all_boundary_conditions.extend(bc.flux_bcs) + else: + all_boundary_conditions.append(bc) + return all_boundary_conditions + def initialise(self): self.create_species_from_traps() self.define_function_spaces(element_degree=self.settings.element_degree) @@ -641,17 +652,9 @@ 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 + """Defines the boundary conditions of the model""" - # 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: + for bc in self._unpacked_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) @@ -1322,6 +1325,13 @@ def convert_advection_term_to_fenics_objects(self): function_space=V, t=self.t ) + def define_boundary_conditions(self): + for bc in self._unpacked_bcs: + if isinstance(bc, boundary_conditions.ParticleFluxBC): + bc._volume_subdomain = self.surface_to_volume[bc.subdomain] + + super().define_boundary_conditions() + def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain): """ Creates the variational formulation for each subdomain and stores it in diff --git a/test/test_complex_reaction.py b/test/test_complex_reaction.py index 3efc8ddbb..bcbbb7923 100644 --- a/test/test_complex_reaction.py +++ b/test/test_complex_reaction.py @@ -174,3 +174,43 @@ def test_reaction(k, p, c_A_0): model = model_test_reaction(stepsize=0.5, k=k, p=p, c_A_0=c_A_0) error = compute_error(model) assert error < 1e-2 + + +def test_surface_reaction_in_discontinuous_case(): + """Test the SurfaceReactionBC in a discontinuous problem setup""" + my_model = F.HydrogenTransportProblemDiscontinuous() + + my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=1000)) + + material = F.Material(D_0=1e-2, E_D=0, K_S_0=2.0, E_K_S=0) + subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=material) + left_surface = F.SurfaceSubdomain1D(id=1, x=0) + right_surface = F.SurfaceSubdomain1D(id=2, x=1) + + my_model.subdomains = [subdomain, left_surface, right_surface] + my_model.surface_to_volume = { + right_surface: subdomain, + left_surface: subdomain, + } + + H = F.Species("H", mobile=True, subdomains=[subdomain]) + my_model.species = [H] + + my_model.boundary_conditions = [ + F.FixedConcentrationBC(left_surface, value=100, species=H), + F.SurfaceReactionBC( + reactant=[H, H], + gas_pressure=1, + k_r0=1, + E_kr=0.2, + k_d0=0, + E_kd=0, + subdomain=right_surface, + ), + ] + + my_model.temperature = 500 + + my_model.settings = F.Settings(atol=1e-5, rtol=1e-5, transient=False) + + my_model.initialise()