From 1d441a6eb1572a32c06291f3129d119366bfce03 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 18 Oct 2023 11:28:35 -0400 Subject: [PATCH] Address review comments --- festim/boundary_conditions/dirichlet_bc.py | 33 ++++++++++++-- festim/helpers.py | 53 ---------------------- festim/hydrogen_transport_problem.py | 16 +++---- test/test_dirichlet_bc.py | 21 ++------- 4 files changed, 42 insertions(+), 81 deletions(-) diff --git a/festim/boundary_conditions/dirichlet_bc.py b/festim/boundary_conditions/dirichlet_bc.py index 6c1dff021..efb4d4dba 100644 --- a/festim/boundary_conditions/dirichlet_bc.py +++ b/festim/boundary_conditions/dirichlet_bc.py @@ -22,6 +22,8 @@ class DirichletBC: species (str): the name of the species value_fenics (fem.Function or fem.Constant): the value of the boundary condition in fenics format + bc_expr (fem.Expression): the expression of the boundary condition that is used to + update the value_fenics Usage: >>> from festim import DirichletBC @@ -60,21 +62,43 @@ def define_surface_subdomain_dofs(self, facet_meshtags, mesh, function_space): condition Args: - mesh (festim.Mesh): the domain mesh + facet_meshtags (ddolfinx.mesh.MeshTags): the mesh tags of the + surface facets + mesh (dolfinx.mesh.Mesh): the mesh + function_space (dolfinx.fem.FunctionSpace): the function space """ bc_facets = facet_meshtags.find(self.subdomain.id) bc_dofs = fem.locate_dofs_topological(function_space, mesh.fdim, bc_facets) return bc_dofs - def create_value(self, mesh, function_space, temperature, t): + def create_value( + self, mesh, function_space: fem.FunctionSpace, temperature, t: fem.Constant + ): + """Creates the value of the boundary condition as a fenics object and sets it to + self.value_fenics. + If the value is a constant, it is converted to a fenics.Constant. + If the value is a function of t, it is converted to a fenics.Constant. + Otherwise, it is converted to a fenics.Function and the + expression of the function is stored in self.bc_expr. + + Args: + mesh (dolfinx.mesh.Mesh) : the mesh + function_space (dolfinx.fem.FunctionSpace): the function space + temperature (float): the temperature + t (dolfinx.fem.Constant): the time + """ x = ufl.SpatialCoordinate(mesh) + if isinstance(self.value, (int, float)): self.value_fenics = F.as_fenics_constant(mesh=mesh, value=self.value) + elif callable(self.value): arguments = self.value.__code__.co_varnames + if "t" in arguments and "x" not in arguments and "T" not in arguments: + # only t is an argument self.value_fenics = F.as_fenics_constant( - mesh=mesh, value=self.value(t=float(t)) + mesh=mesh, value=self.value(t=t.value) ) else: self.value_fenics = fem.Function(function_space) @@ -85,6 +109,9 @@ def create_value(self, mesh, function_space, temperature, t): kwargs["x"] = x if "T" in arguments: kwargs["T"] = temperature + + # store the expression of the boundary condition + # to update the value_fenics later self.bc_expr = fem.Expression( self.value(**kwargs), function_space.element.interpolation_points(), diff --git a/festim/helpers.py b/festim/helpers.py index 96e0f479f..9b17fb92b 100644 --- a/festim/helpers.py +++ b/festim/helpers.py @@ -1,7 +1,4 @@ from dolfinx import fem -import ufl -import numpy as np -from petsc4py.PETSc import ScalarType def as_fenics_constant(value, mesh): @@ -25,53 +22,3 @@ def as_fenics_constant(value, mesh): raise TypeError( f"Value must be a float, an int or a dolfinx.Constant, not {type(value)}" ) - - -# class SpaceTimeDependentExpression: -# def __init__(self, function, t=0): -# self.t = t -# self.function = function -# self.values = None - -# def __call__(self, x): -# if self.values is None: -# self.values = np.zeros(x.shape[1], dtype=ScalarType) -# self.values = np.full(x.shape[1], self.function(x=x, t=self.t)) -# return self.values - - -# def convert_to_appropriate_obj(object, function_space, mesh): -# """Converts a value to a dolfinx.Constant or a dolfinx.Function -# depending on the type of the value - -# Args: -# object (callable or float): the value to convert -# function_space (dolfinx.fem.FunctionSpace): the function space of the domain -# mesh (dolfinx.mesh.mesh): the mesh of the domain - -# Returns: -# dolfinx.Constant or dolfinx.Function: the converted value -# festim.SpaceTimeDependentExpression or None: the expression if the value is -# space and time dependent, None otherwise -# """ -# x = ufl.SpatialCoordinate(mesh) -# if isinstance(object, (int, float)): -# # case 1 pressure isn't space dependent or only time dependent: -# return as_fenics_constant(mesh=mesh, value=object), None -# # case 2 pressure is space dependent -# elif callable(object): -# arguments = object.__code__.co_varnames -# if "t" in arguments and "x" in arguments: -# expr = fem.Expression( -# object(x=x, t=0), function_space.element.interpolation_points() -# ) -# fenics_obj = fem.Function(function_space) -# fenics_obj.interpolate(expr) -# return fenics_obj, expr -# elif "x" in arguments: -# fenics_obj = fem.Function(function_space) -# fenics_obj.interpolate(object) -# return fenics_obj, None - -# elif "t" in arguments: -# return as_fenics_constant(mesh=mesh, value=object(t=0)), None diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py index 6f64535ce..ae2fd9b64 100644 --- a/festim/hydrogen_transport_problem.py +++ b/festim/hydrogen_transport_problem.py @@ -41,8 +41,8 @@ class HydrogenTransportProblem: ds (dolfinx.fem.ds): the surface measure of the model function_space (dolfinx.fem.FunctionSpace): the function space of the model - facet_meshtags (dolfinx.cpp.mesh.MeshTags): the facet tags of the model - volume_meshtags (dolfinx.cpp.mesh.MeshTags): the volume tags of the + facet_meshtags (dolfinx.mesh.MeshTags): the facet tags of the model + volume_meshtags (dolfinx.mesh.MeshTags): the volume tags of the model formulation (ufl.form.Form): the formulation of the model solver (dolfinx.nls.newton.NewtonSolver): the solver of the model @@ -283,23 +283,23 @@ def run(self, final_time: float): progress = tqdm.autonotebook.tqdm( desc="Solving H transport problem", total=final_time ) - while float(self.t) < final_time: - progress.update(float(self.dt)) - self.t.value += float(self.dt) + while self.t.value < final_time: + progress.update(self.dt.value) + self.t.value += self.dt.value # update boundary conditions for bc in self.boundary_conditions: - bc.update(float(self.t)) + bc.update(self.t.value) self.solver.solve(self.u) - mobile_xdmf.write_function(self.u, float(self.t)) + mobile_xdmf.write_function(self.u, self.t.value) surface_flux = form(D * dot(grad(cm), n) * self.ds(2)) flux = assemble_scalar(surface_flux) flux_values.append(flux) - times.append(float(self.t)) + times.append(self.t.value) # update previous solution self.u_n.x.array[:] = self.u.x.array[:] diff --git a/test/test_dirichlet_bc.py b/test/test_dirichlet_bc.py index 0a003ab95..cb299ac39 100644 --- a/test/test_dirichlet_bc.py +++ b/test/test_dirichlet_bc.py @@ -14,6 +14,7 @@ def test_init(): + """Test that the attributes are set correctly""" # create a DirichletBC object subdomain = F.SurfaceSubdomain1D(1, x=0) value = 1.0 @@ -29,6 +30,9 @@ def test_init(): def test_value_fenics(): + """Test that the value_fenics attribute can be set to a valid value + and that an invalid type throws an error + """ # create a DirichletBC object subdomain = F.SurfaceSubdomain1D(1, x=0) value = 1.0 @@ -47,23 +51,6 @@ def test_value_fenics(): bc.value_fenics = "invalid" -# def test_define_surface_subdomain_dofs(): -# # create a DirichletBC object -# subdomain = F.SurfaceSubdomain1D(1, x=0) -# value = 1.0 -# species = "test" -# bc = F.DirichletBC(subdomain, value, species) - -# # create facet meshtags, mesh, and function space objects -# facet_meshtags = ..... -# function_space = fem.FunctionSpace(mesh, ("Lagrange", 1)) - -# # call the method being tested -# bc_dofs = bc.define_surface_subdomain_dofs(facet_meshtags, mesh, function_space) - -# assert bc_dofs ... - - def test_callable_for_value(): """Test that the value attribute can be a callable function of x and t"""