diff --git a/src/festim/__init__.py b/src/festim/__init__.py index b2095f861..6e0849464 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -59,6 +59,7 @@ ) from .hydrogen_transport_problem import ( HydrogenTransportProblem, + HydrogenTransportProblemDG, HydrogenTransportProblemDiscontinuous, HydrogenTransportProblemDiscontinuousChangeVar, ) diff --git a/src/festim/boundary_conditions/__init__.py b/src/festim/boundary_conditions/__init__.py index 1df3d9c7a..72a320d61 100644 --- a/src/festim/boundary_conditions/__init__.py +++ b/src/festim/boundary_conditions/__init__.py @@ -1,16 +1,24 @@ -from .dirichlet_bc import DirichletBCBase, FixedConcentrationBC, FixedTemperatureBC +from .dirichlet_bc import ( + DirichletBCBase, + FixedConcentrationBC, + FixedConcentrationInflowBC, + FixedTemperatureBC, +) from .flux_bc import FluxBCBase, HeatFluxBC, ParticleFluxBC from .henrys_bc import HenrysBC +from .outflow_bc import OutflowBC from .sieverts_bc import SievertsBC from .surface_reaction import SurfaceReactionBC, SurfaceReactionBCpartial __all__ = [ "DirichletBCBase", "FixedConcentrationBC", + "FixedConcentrationInflowBC", "FixedTemperatureBC", "FluxBCBase", "HeatFluxBC", "HenrysBC", + "OutflowBC", "ParticleFluxBC", "SievertsBC", "SurfaceReactionBC", diff --git a/src/festim/boundary_conditions/dirichlet_bc.py b/src/festim/boundary_conditions/dirichlet_bc.py index fde8d4527..3b00470c7 100644 --- a/src/festim/boundary_conditions/dirichlet_bc.py +++ b/src/festim/boundary_conditions/dirichlet_bc.py @@ -12,6 +12,7 @@ from festim import helpers from festim import subdomain as _subdomain +from festim.advection import VelocityField from festim.species import Species @@ -421,3 +422,53 @@ def create_value(self, function_space: fem.FunctionSpace, t: fem.Constant): helpers.get_interpolation_points(function_space.element), ) self.value_fenics.interpolate(self.bc_expr) + + +class FixedConcentrationInflowBC(FixedConcentrationBC): + """Inflow boundary condition with fixed concentration. + + This is a specific case of `FixedConcentrationBC` where the value of the + boundary condition is fixed and does not depend on time or space. It is + used for inflow boundary conditions where the concentration of the species + entering the domain is known and constant. + + Args: + subdomain (festim.Subdomain): the surface subdomain where the boundary + condition is applied + value: The value of the boundary condition. It can be a function of + space and/or time + species: The name of the species + + """ + + def __init__( + self, + subdomain: _subdomain.SurfaceSubdomain, + value: np.ndarray | fem.Constant | int | float | Callable, + species: Species, + velocity_field: fem.Function, + ): + super().__init__(subdomain, value, species, enforce_weakly=False) + + self.velocity_field = velocity_field + + @property + def velocity_field(self): + return self._velocity_field + + @velocity_field.setter + def velocity_field(self, value): + err_message = f"velocity must be a fem.Function, or callable not {type(value)}" + if value is None: + self._velocity = VelocityField(value) + elif isinstance( + value, + fem.Function, + ): + self._velocity = VelocityField(value) + elif isinstance(value, fem.Constant | fem.Expression | ufl.core.expr.Expr): + raise TypeError(err_message) + elif callable(value): + self._velocity = VelocityField(value) + else: + raise TypeError(err_message) diff --git a/src/festim/boundary_conditions/outflow_bc.py b/src/festim/boundary_conditions/outflow_bc.py new file mode 100644 index 000000000..4c24f0e7d --- /dev/null +++ b/src/festim/boundary_conditions/outflow_bc.py @@ -0,0 +1,110 @@ +import ufl +from dolfinx import fem + +from festim import subdomain as _subdomain +from festim.advection import VelocityField +from festim.species import Species + + +class OutflowBC: + """ + Outflow boundary condition class + + dot(velocity, n) * u * v * ds + + This term allows the species to flow across the boundary due to advection. It arises + from integrating the conservative form of the advection term by parts in the weak + formulation. + + Args: + velocity: The velocity field at the outflow boundary condition + species: The species for which the outflow boundary condition is applied, can be + a list of species or a single species + subdomain: The surface subdomain where the boundary condition is applied + + Attributes: + velocity: The velocity field at the outflow boundary condition + species: The species for which the outflow boundary condition is applied, can be + a list of species or a single species + subdomain: The surface subdomain where the boundary condition is applied + + """ + + velocity: VelocityField + species: list[Species] + subdomain: _subdomain.SurfaceSubdomain + + def __init__( + self, + velocity: fem.Function, + species: Species | list[Species], + subdomain: _subdomain.SurfaceSubdomain, + ): + self.subdomain = subdomain + self.velocity = velocity + self.species = species + + @property + def velocity(self): + return self._velocity + + @velocity.setter + def velocity(self, value): + err_message = f"velocity must be a fem.Function, or callable not {type(value)}" + if value is None: + self._velocity = VelocityField(value) + elif isinstance( + value, + fem.Function, + ): + self._velocity = VelocityField(value) + elif isinstance(value, fem.Constant | fem.Expression | ufl.core.expr.Expr): + raise TypeError(err_message) + elif callable(value): + self._velocity = VelocityField(value) + else: + raise TypeError(err_message) + + @property + def species(self) -> list[Species]: + return self._species + + @species.setter + def species(self, value): + if not isinstance(value, list): + value = [value] + # check that all species are of type festim.Species + for spe in value: + if not isinstance(spe, Species): + raise TypeError( + f"elements of species must be of type festim.Species not " + f"{type(spe)}" + ) + self._species = value + + @property + def subdomain(self): + return self._subdomain + + @subdomain.setter + def subdomain(self, value): + if value is None: + self._subdomain = value + elif isinstance(value, _subdomain.SurfaceSubdomain): + self._subdomain = value + else: + raise TypeError( + f"Subdomain must be a festim.SurfaceSubdomain object, not {type(value)}" + ) + + @property + def time_dependent(self): + if self.velocity is None: + raise TypeError("Value must be given to determine if its time dependent") + if isinstance(self.velocity, fem.Constant): + return False + if callable(self.velocity): + arguments = self.velocity.__code__.co_varnames + return "t" in arguments + else: + return False diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 1e6682ccb..6c20bbb08 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -46,6 +46,7 @@ __all__ = [ "HydrogenTransportProblem", + "HydrogenTransportProblemDG", "HydrogenTransportProblemDiscontinuous", ] @@ -2273,3 +2274,287 @@ def update_time_dependent_values(self): for bc in self.boundary_conditions: if isinstance(bc, boundary_conditions.FixedConcentrationBC): bc.update(self.t) + + +class HydrogenTransportProblemDG(HydrogenTransportProblem): + def __init__( + self, + mesh=None, + subdomains=None, + species=None, + reactions=None, + temperature=None, + sources=None, + initial_conditions=None, + boundary_conditions=None, + settings=None, + exports=None, + traps=None, + petsc_options: dict | None = None, + ): + super().__init__( + mesh, + subdomains, + species, + reactions, + temperature, + sources, + initial_conditions, + boundary_conditions, + settings, + exports, + traps, + petsc_options=petsc_options, + ) + + def initialise(self): + self.create_species_from_traps() + self.define_function_spaces(element_degree=self.settings.element_degree) + self.define_meshtags_and_measures() + self.assign_functions_to_species() + + self.t = fem.Constant(self.mesh.mesh, 0.0) + if self.settings.transient: + # TODO should raise error if no stepsize is provided + # TODO Should this be an attribute of festim.Stepsize? + self._dt = as_fenics_constant( + self.settings.stepsize.initial_value, self.mesh.mesh + ) + + self.create_implicit_species_value_fenics() + + self.define_temperature() + self.define_boundary_conditions() + self.convert_source_input_values_to_fenics_objects() + self.convert_advection_term_to_fenics_objects() + self.create_flux_values_fenics() + self.create_initial_conditions() + self.create_formulation() + self.create_solver() + self.initialise_exports() + + def define_function_spaces(self, element_degree: int = 1): + """Creates the function space of the modelw with a mixed element. Creates the + main solution and previous solution function u and u_n. Create global DG + function spaces of degree 0 and 1 for the global diffusion coefficient. + + Args: + element_degree: Degree order for finite element. Defaults to 1. + """ + + element_DG = basix.ufl.element( + "DG", + self.mesh.mesh.basix_cell(), + element_degree, + basix.LagrangeVariant.equispaced, + ) + + elements = [] + for spe in self.species: + if isinstance(spe, _species.Species): + elements.append(element_DG) + + element = basix.ufl.mixed_element(elements) + + self.function_space = fem.functionspace(self.mesh.mesh, element) + + # create global DG function spaces of degree 0 and 1 + element_DG0 = basix.ufl.element( + "DG", + self.mesh.mesh.basix_cell(), + 0, + basix.LagrangeVariant.equispaced, + ) + element_DG1 = basix.ufl.element( + "DG", + self.mesh.mesh.basix_cell(), + 1, + basix.LagrangeVariant.equispaced, + ) + self.V_DG_0 = fem.functionspace(self.mesh.mesh, element_DG0) + self.V_DG_1 = fem.functionspace(self.mesh.mesh, element_DG1) + + self.u = fem.Function(self.function_space) + self.u_n = fem.Function(self.function_space) + + def create_formulation(self): + """Creates the formulation of the model.""" + + self.formulation = 0 + + # add diffusion and time derivative for each species + for spe in self.species: + u = spe.solution + u_n = spe.prev_solution + v = spe.test_function + + for vol in self.volume_subdomains: + if spe.mobile: + D = vol.material.get_diffusion_coefficient( + self.mesh.mesh, self.temperature_fenics, spe + ) + match self.mesh.coordinate_system: + case CoordinateSystem.CARTESIAN: + self.formulation += ufl.dot( + D * ufl.grad(u), ufl.grad(v) + ) * self.dx(vol.id) + case CoordinateSystem.CYLINDRICAL: + r = ufl.SpatialCoordinate(self.mesh.mesh)[0] + self.formulation += ( + r + * ufl.dot(D * ufl.grad(u), ufl.grad(v / r)) + * self.dx(vol.id) + ) + case CoordinateSystem.SPHERICAL: + r = ufl.SpatialCoordinate(self.mesh.mesh)[0] + self.formulation += ( + r**2 + * ufl.dot(D * ufl.grad(u), ufl.grad(v / r**2)) + * self.dx(vol.id) + ) + case _: + raise NotImplementedError( + f"Unknown coordinate system {self.mesh.coordinate_system!s}" # noqa: E501 + ) + + # Add SIPG diffusion + self.formulation += ( + -D + * ufl.inner(ufl.avg(ufl.grad(u)), ufl.jump(v, self.mesh.n)) + * self.dS + ) + self.formulation += ( + -D + * ufl.inner(ufl.jump(u, self.mesh.n), ufl.avg(ufl.grad(v))) + * self.dS + ) + PENALY = 100 + self.formulation += ( + D + * (PENALY / ufl.avg(self.mesh.h)) + * ufl.inner(ufl.jump(u, self.mesh.n), ufl.jump(v, self.mesh.n)) + * self.dS + ) + + if self.settings.transient: + self.formulation += ((u - u_n) / self.dt) * v * self.dx(vol.id) + + for reaction in self.reactions: + for reactant in reaction.reactant: + if isinstance(reactant, _species.Species): + self.formulation += ( + reaction.reaction_term(self.temperature_fenics) + * reactant.test_function + * self.dx(reaction.volume.id) + ) + + # product + if isinstance(reaction.product, list): + products = reaction.product + else: + products = [reaction.product] + for product in products: + self.formulation += ( + -reaction.reaction_term(self.temperature_fenics) + * product.test_function + * self.dx(reaction.volume.id) + ) + # add sources + for source in self.sources: + self.formulation -= ( + source.value.fenics_object + * source.species.test_function + * self.dx(source.volume.id) + ) + + # add fluxes + for bc in self.boundary_conditions: + if isinstance(bc, boundary_conditions.ParticleFluxBC): + self.formulation -= ( + bc.value_fenics + * 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) + ) + if isinstance(bc, boundary_conditions.FixedConcentrationBC): + if not bc.enforce_weakly: + raise ValueError( + "FixedConcentrationBC must be enforced weakly for DG formulation" + ) + u = bc.species.solution + v = bc.species.test_function + self.formulation += bc.weak_formulation(u, v, self.ds) + + if isinstance(bc, boundary_conditions.FixedConcentrationInflowBC): + if not bc.enforce_weakly: + raise ValueError( + "FixedConcentrationInflowBC must be enforced weakly " + "for DG formulation" + ) + u = bc.species.solution + v = bc.species.test_function + vel = bc.velocity.fenics_object + + lmbda = ufl.conditional(ufl.gt(ufl.dot(vel, self.mesh.n), 0), 1, 0) + + self.formulation += bc.weak_formulation(u, v, self.ds) + + self.formulation += -ufl.inner( + (1 - lmbda) * ufl.dot(vel, self.mesh.n) * u, v + ) * self.ds(bc.subdomain.id) + + if isinstance(bc, boundary_conditions.OutflowBC): + for species in bc.species: + conc = species.solution + v = species.test_function + vel = bc.velocity.fenics_object + + lmbda = ufl.conditional(ufl.gt(ufl.dot(vel, self.mesh.n), 0), 1, 0) + + self.formulation += ( + ufl.inner(lmbda * ufl.dot(vel, self.mesh.n) * u, v) * self.ds + ) + + for adv_term in self.advection_terms: + for species in adv_term.species: + conc = species.solution + v = species.test_function + vel = adv_term.velocity.fenics_object + + lmbda = ufl.conditional(ufl.gt(ufl.dot(vel, self.mesh.n), 0), 1, 0) + + # Advection, bulk + self.formulation += -ufl.inner(vel * conc, ufl.grad(v)) * self.dx( + adv_term.subdomain.id + ) + + # Advection, interior upwind flux + self.formulation += ( + ufl.inner(2 * ufl.avg(lmbda * vel * conc), ufl.jump(v, self.mesh.n)) + * self.dS + ) + + # check if each species is defined in all volumes + if not self.settings.transient: + for spe in self.species: + # if species mobile, already defined in diffusion term + if not spe.mobile: + not_defined_in_volume = self.volume_subdomains.copy() + for vol in self.volume_subdomains: + # check reactions + for reaction in self.reactions: + if vol == reaction.volume: + if vol in not_defined_in_volume: + not_defined_in_volume.remove(vol) + + # add c = 0 to formulation where needed + for vol in not_defined_in_volume: + self.formulation += ( + spe.solution * spe.test_function * self.dx(vol.id) + ) diff --git a/src/festim/mesh/mesh.py b/src/festim/mesh/mesh.py index 76fa9003f..0a3d983f1 100644 --- a/src/festim/mesh/mesh.py +++ b/src/festim/mesh/mesh.py @@ -62,6 +62,7 @@ class Mesh: vdim: int fdim: int n: ufl.FacetNormal + h: ufl.CellDiameter def __init__( self, @@ -127,6 +128,12 @@ def n(self): raise RuntimeError("Mesh is not defined") return ufl.FacetNormal(self._mesh) + @property + def h(self): + if self._mesh is None: + raise RuntimeError("Mesh is not defined") + return ufl.CellDiameter(self._mesh) + def define_meshtags(self, surface_subdomains, volume_subdomains, interfaces=None): """Defines the facet and volume meshtags of the mesh. diff --git a/src/festim/problem.py b/src/festim/problem.py index 7ee91e070..892e67b49 100644 --- a/src/festim/problem.py +++ b/src/festim/problem.py @@ -118,6 +118,7 @@ def define_meshtags_and_measures(self): self.dx = ufl.Measure( "dx", domain=self.mesh.mesh, subdomain_data=self.volume_meshtags ) + self.dS = ufl.Measure("dS", domain=self.mesh.mesh) def define_boundary_conditions(self): """Defines the dirichlet boundary conditions of the model."""