diff --git a/.github/workflows/ci_docker.yml b/.github/workflows/ci_docker.yml index 7eec50a65..4dbf7bfd7 100644 --- a/.github/workflows/ci_docker.yml +++ b/.github/workflows/ci_docker.yml @@ -20,7 +20,8 @@ jobs: - name: Install local package and dependencies run: | - pip install .[test] + python -m pip install scipy # needed in scifem, can be removed when scifem adds it to their dependencies + python -m pip install .[test] - name: Overload adios4dolfinx if: ${{ matrix.container_version == 'nightly' }} diff --git a/pyproject.toml b/pyproject.toml index fdf0ab022..50046b102 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,9 @@ [build-system] -requires = ["setuptools >= 61", "wheel", "setuptools-scm[toml] >= 7.0.5"] +requires = [ + "setuptools>=77.0.3", # Ensuring license field is correctly parsed + "wheel", + "setuptools-scm[toml] >= 7.0.5", +] build-backend = "setuptools.build_meta" [project] diff --git a/src/festim/boundary_conditions/dirichlet_bc.py b/src/festim/boundary_conditions/dirichlet_bc.py index 4ae5604cd..fde8d4527 100644 --- a/src/festim/boundary_conditions/dirichlet_bc.py +++ b/src/festim/boundary_conditions/dirichlet_bc.py @@ -3,8 +3,10 @@ import numpy as np import numpy.typing as npt import ufl +import ufl.argument import ufl.core import ufl.core.expr +import ufl.indexed from dolfinx import fem from dolfinx import mesh as _mesh @@ -19,6 +21,10 @@ class DirichletBCBase: Args: subdomain: The surface subdomain where the boundary condition is applied value: The value of the boundary condition + enforce_weakly: Whether to enforce the boundary condition weakly using Nitsche's + method. Defaults to False. + penalty: The penalty parameter to use if ``enforce_weakly`` is True. + Defaults to None Attributes: subdomain: The surface subdomain where the boundary condition is applied @@ -26,6 +32,9 @@ class DirichletBCBase: value_fenics: The value of the boundary condition in fenics format bc_expr: The expression of the boundary condition that is used to update the `value_fenics` + enforce_weakly: Whether to enforce the boundary condition weakly using Nitsche's + method. + penalty: The penalty parameter to use if ``enforce_weakly`` is True. """ subdomain: _subdomain.SurfaceSubdomain @@ -45,9 +54,13 @@ def __init__( self, subdomain: _subdomain.SurfaceSubdomain, value: np.ndarray | fem.Constant | int | float | Callable, + enforce_weakly: bool = False, + penalty: float | None = None, ): self.subdomain = subdomain self.value = value + self.enforce_weakly = enforce_weakly + self.penalty = penalty self.value_fenics = None self.bc_expr = None @@ -138,6 +151,47 @@ def update(self, t: float): elif self.bc_expr is not None: self.value_fenics.interpolate(self.bc_expr) + def weak_formulation( + self, + u: fem.Function | ufl.indexed.Indexed, + v: ufl.argument.Argument | ufl.indexed.Indexed, + ds: ufl.Measure, + ) -> ufl.core.expr.Expr: + """ + Returns the Nitsche weak formulation for the BC + This follows the dolfinx tutorial + https://jsdokken.com/dolfinx-tutorial/chapter1/nitsche.html + + Args: + u: the solution function associated to the species for which the BC + is applied + v: the test function + ds: the surface measure + + Returns: + the weak formulation + """ + mesh = ds.ufl_domain() + n = ufl.FacetNormal(mesh) + h = ufl.Circumradius(mesh) # FIXME this doesn't work for rectangles + alpha = self.penalty + assert alpha is not None, ( + "Penalty parameter must be given for weakly enforced Dirichlet BCs" + ) + assert self.value_fenics is not None, ( + "value_fenics must be defined for weakly enforced Dirichlet BCs" + ) + + form = -ufl.inner(n, ufl.grad(u)) * v * ds(self.subdomain.id) + + form += ( + +ufl.inner(n, ufl.grad(v)) * (u - self.value_fenics) * ds(self.subdomain.id) + ) + form += ( + -alpha / h * ufl.inner((u - self.value_fenics), v) * ds(self.subdomain.id) + ) + return form + class FixedConcentrationBC(DirichletBCBase): """ @@ -147,6 +201,10 @@ class FixedConcentrationBC(DirichletBCBase): value: The value of the boundary condition. It can be a function of space and/or time species: The name of the species + enforce_weakly: Whether to enforce the boundary condition weakly using Nitsche's + method. Defaults to False. + penalty: The penalty parameter to use if ``enforce_weakly`` is True. + Defaults to None Attributes: temperature_dependent (bool): True if the value of the bc is @@ -180,9 +238,13 @@ def __init__( subdomain: _subdomain.SurfaceSubdomain, value: np.ndarray | fem.Constant | int | float | Callable, species: Species, + enforce_weakly: bool = False, + penalty: float | None = None, ): self.species = species - super().__init__(subdomain, value) + super().__init__( + subdomain, value, enforce_weakly=enforce_weakly, penalty=penalty + ) @property def temperature_dependent(self): diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 7b181366e..1e6682ccb 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -878,7 +878,7 @@ def create_formulation(self): * self.dx(source.volume.id) ) - # add fluxes + # add boundary conditions (fluxes and weak dirichlet) for bc in self.boundary_conditions: if isinstance(bc, boundary_conditions.ParticleFluxBC): self.formulation -= ( @@ -894,6 +894,11 @@ def create_formulation(self): * self.ds(flux_bc.subdomain.id) ) + if isinstance(bc, boundary_conditions.FixedConcentrationBC): + if bc.enforce_weakly: + u = bc.species.solution + v = bc.species.test_function + self.formulation += bc.weak_formulation(u, v, self.ds) for adv_term in self.advection_terms: # create vector functionspace based on the elements in the mesh @@ -1549,6 +1554,11 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain): if subdomain == self.surface_to_volume[bc.subdomain]: v = bc.species.subdomain_to_test_function[subdomain] form -= bc.value_fenics * v * self.ds(bc.subdomain.id) + if isinstance(bc, boundary_conditions.FixedConcentrationBC): + if bc.enforce_weakly: + u = bc.species.subdomain_to_solution[subdomain] + v = bc.species.subdomain_to_test_function[subdomain] + form += bc.weak_formulation(u, v, self.ds) # add volumetric sources for source in self.sources: diff --git a/src/festim/problem.py b/src/festim/problem.py index 675ef7474..7ee91e070 100644 --- a/src/festim/problem.py +++ b/src/festim/problem.py @@ -124,7 +124,8 @@ def define_boundary_conditions(self): for bc in self.boundary_conditions: if isinstance(bc, F.DirichletBCBase): form = self.create_dirichletbc_form(bc) - self.bc_forms.append(form) + if not bc.enforce_weakly: + self.bc_forms.append(form) def get_petsc_options(self) -> dict[str, Any]: """Gets the PETSc options to pass to the NewtonProblem solver. Default options diff --git a/test/system_tests/test_misc.py b/test/system_tests/test_misc.py index 8d678f43f..40a046c12 100644 --- a/test/system_tests/test_misc.py +++ b/test/system_tests/test_misc.py @@ -3,10 +3,12 @@ import dolfinx import numpy as np import pytest +import ufl import festim as F from .test_multi_mat_penalty import generate_mesh +from .tools import error_L2 def test_petsc_options(): @@ -389,3 +391,90 @@ def test_sub_temperature_as_function_mixed_domain_not_updated(): my_model.run() assert not np.allclose(avg_surf.data, 4.0) + + +@pytest.mark.parametrize( + "model_class", [F.HydrogenTransportProblem, F.HydrogenTransportProblemDiscontinuous] +) +def test_MMS_weak_dirichlet( + model_class: type[F.HydrogenTransportProblem] + | type[F.HydrogenTransportProblemDiscontinuous], +): + """ + MMS test with one mobile species at steady state, with Dirichlet BCs enforced + weakly. + """ + + def u_exact(mod): + return lambda x: 1 + mod.sin(2 * mod.pi * x[0]) + mod.cos(2 * mod.pi * x[1]) + + H_analytical_ufl = u_exact(ufl) + H_analytical_np = u_exact(np) + + mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 50, 50) + + V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1)) + T = dolfinx.fem.Function(V) + + D_0 = 1 + E_D = 0.1 + + def T_expr(x): + return 500 + 100 * x[0] + + T.interpolate(T_expr) + D = D_0 * ufl.exp(-E_D / (F.k_B * T)) + + my_model = model_class() + my_model.mesh = F.Mesh(mesh=mesh) + + my_mat = F.Material(name="mat", D_0=D_0, E_D=E_D) + vol = F.VolumeSubdomain(id=1, material=my_mat) + left = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[0], 0)) + right = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 1)) + my_model.subdomains = [vol, left, right] + + H = F.Species("H") + H2 = F.Species("D") + my_model.species = [H, H2] + + if isinstance(my_model, F.HydrogenTransportProblemDiscontinuous): + H.subdomains = [vol] + H2.subdomains = [vol] + + my_model.temperature = T_expr + + my_model.boundary_conditions = [ + F.FixedConcentrationBC( + subdomain=left, + value=H_analytical_ufl, + species=H, + enforce_weakly=True, + penalty=300, + ), + F.FixedConcentrationBC( + subdomain=right, + value=H_analytical_ufl, + species=H, + enforce_weakly=True, + penalty=300, + ), + ] + + x = ufl.SpatialCoordinate(mesh) + f = -ufl.div(D * ufl.grad(H_analytical_ufl(x))) + my_model.sources = [F.ParticleSource(value=f, volume=vol, species=H)] + + my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False) + + my_model.initialise() + my_model.run() + + if isinstance(my_model, F.HydrogenTransportProblemDiscontinuous): + H_computed = H.subdomain_to_post_processing_solution[vol] + else: + H_computed = H.post_processing_solution + + L2_error = error_L2(H_computed, H_analytical_np) + + assert L2_error < 2e-3