From 06b1dabc89f2f6004ba3734acd280a7830b28c59 Mon Sep 17 00:00:00 2001 From: jhdark Date: Wed, 22 Jan 2025 10:57:16 -0500 Subject: [PATCH 01/40] new Convert to fenics class --- src/festim/__init__.py | 2 +- .../boundary_conditions/dirichlet_bc.py | 218 +++++++++--------- src/festim/boundary_conditions/flux_bc.py | 138 +---------- src/festim/helpers.py | 178 ++++++++++++++ src/festim/hydrogen_transport_problem.py | 178 +++++++------- src/festim/initial_condition.py | 63 +---- src/festim/problem.py | 8 +- src/festim/source.py | 143 +----------- test/test_h_transport_problem.py | 28 ++- 9 files changed, 407 insertions(+), 549 deletions(-) diff --git a/src/festim/__init__.py b/src/festim/__init__.py index be4cdd6ed..5edea51b3 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -36,7 +36,7 @@ from .exports.vtx import VTXSpeciesExport, VTXTemperatureExport from .exports.xdmf import XDMFExport from .heat_transfer_problem import HeatTransferProblem -from .helpers import as_fenics_constant +from .helpers import as_fenics_constant, ConvertToFenicsObject from .hydrogen_transport_problem import ( HTransportProblemDiscontinuous, HydrogenTransportProblem, diff --git a/src/festim/boundary_conditions/dirichlet_bc.py b/src/festim/boundary_conditions/dirichlet_bc.py index 0a5af01d6..67d8dc670 100644 --- a/src/festim/boundary_conditions/dirichlet_bc.py +++ b/src/festim/boundary_conditions/dirichlet_bc.py @@ -49,38 +49,38 @@ def __init__( self.subdomain = subdomain self.value = value - self.value_fenics = None - self.bc_expr = None - - @property - def value_fenics(self): - return self._value_fenics - - @value_fenics.setter - def value_fenics(self, value: None | fem.Function | fem.Constant | np.ndarray): - if value is None: - self._value_fenics = value - return - if not isinstance(value, (fem.Function, fem.Constant, np.ndarray)): - # FIXME: Should we allow sending in a callable here? - raise TypeError( - "Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not" - + f"{type(value)}" - ) - self._value_fenics = value - - @property - def time_dependent(self) -> bool: - """Returns true if the value of the boundary condition is time dependent""" - if self.value is None: - return False - if isinstance(self.value, fem.Constant): - return False - if callable(self.value): - arguments = self.value.__code__.co_varnames - return "t" in arguments - else: - return False + self.value_fenics = helpers.ConvertToFenicsObject(value) + # self.bc_expr = None + + # @property + # def value_fenics(self): + # return self._value_fenics + + # @value_fenics.setter + # def value_fenics(self, value: None | fem.Function | fem.Constant | np.ndarray): + # if value is None: + # self._value_fenics = value + # return + # if not isinstance(value, (fem.Function, fem.Constant, np.ndarray)): + # # FIXME: Should we allow sending in a callable here? + # raise TypeError( + # "Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not" + # + f"{type(value)}" + # ) + # self._value_fenics = value + + # @property + # def time_dependent(self) -> bool: + # """Returns true if the value of the boundary condition is time dependent""" + # if self.value is None: + # return False + # if isinstance(self.value, fem.Constant): + # return False + # if callable(self.value): + # arguments = self.value.__code__.co_varnames + # return "t" in arguments + # else: + # return False def define_surface_subdomain_dofs( self, @@ -120,18 +120,18 @@ def define_surface_subdomain_dofs( return bc_dofs - def update(self, t: float): - """Updates the boundary condition value + # def update(self, t: float): + # """Updates the boundary condition value - Args: - t (float): the time - """ - if callable(self.value): - arguments = self.value.__code__.co_varnames - if isinstance(self.value_fenics, fem.Constant) and "t" in arguments: - self.value_fenics.value = self.value(t=t) - else: - self.value_fenics.interpolate(self.bc_expr) + # Args: + # t (float): the time + # """ + # if callable(self.value): + # arguments = self.value.__code__.co_varnames + # if isinstance(self.value_fenics, fem.Constant) and "t" in arguments: + # self.value_fenics.value = self.value(t=t) + # else: + # self.value_fenics.interpolate(self.bc_expr) class FixedConcentrationBC(DirichletBCBase): @@ -174,72 +174,72 @@ def __init__( self.species = species super().__init__(subdomain, value) - @property - def temperature_dependent(self): - if self.value is None: - return False - if isinstance(self.value, fem.Constant): - return False - if callable(self.value): - arguments = self.value.__code__.co_varnames - return "T" in arguments - else: - return False - - def create_value( - self, - function_space: fem.FunctionSpace, - temperature: float | fem.Constant, - t: float | 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 `dolfinx.fem.Constant`. - If the value is a function of t, it is converted to `dolfinx.fem.Constant`. - Otherwise, it is converted to a `dolfinx.fem.Function`.Function and the - expression of the function is stored in `bc_expr`. - - Args: - function_space (dolfinx.fem.FunctionSpace): the function space - temperature: The temperature - t (dolfinx.fem.Constant): the time - """ - mesh = function_space.mesh - x = ufl.SpatialCoordinate(mesh) - - if isinstance(self.value, (int, float)): - self.value_fenics = helpers.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 - if not isinstance(self.value(t=float(t)), (float, int)): - raise ValueError( - "self.value should return a float or an int, not " - + f"{type(self.value(t=float(t)))} " - ) - self.value_fenics = helpers.as_fenics_constant( - mesh=mesh, value=self.value(t=float(t)) - ) - else: - self.value_fenics = fem.Function(function_space) - kwargs = {} - if "t" in arguments: - kwargs["t"] = t - if "x" in arguments: - 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(), - ) - self.value_fenics.interpolate(self.bc_expr) + # @property + # def temperature_dependent(self): + # if self.value is None: + # return False + # if isinstance(self.value, fem.Constant): + # return False + # if callable(self.value): + # arguments = self.value.__code__.co_varnames + # return "T" in arguments + # else: + # return False + + # def create_value( + # self, + # function_space: fem.FunctionSpace, + # temperature: float | fem.Constant, + # t: float | 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 `dolfinx.fem.Constant`. + # If the value is a function of t, it is converted to `dolfinx.fem.Constant`. + # Otherwise, it is converted to a `dolfinx.fem.Function`.Function and the + # expression of the function is stored in `bc_expr`. + + # Args: + # function_space (dolfinx.fem.FunctionSpace): the function space + # temperature: The temperature + # t (dolfinx.fem.Constant): the time + # """ + # mesh = function_space.mesh + # x = ufl.SpatialCoordinate(mesh) + + # if isinstance(self.value, (int, float)): + # self.value_fenics = helpers.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 + # if not isinstance(self.value(t=float(t)), (float, int)): + # raise ValueError( + # "self.value should return a float or an int, not " + # + f"{type(self.value(t=float(t)))} " + # ) + # self.value_fenics = helpers.as_fenics_constant( + # mesh=mesh, value=self.value(t=float(t)) + # ) + # else: + # self.value_fenics = fem.Function(function_space) + # kwargs = {} + # if "t" in arguments: + # kwargs["t"] = t + # if "x" in arguments: + # 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(), + # ) + # self.value_fenics.interpolate(self.bc_expr) # alias for FixedConcentrationBC diff --git a/src/festim/boundary_conditions/flux_bc.py b/src/festim/boundary_conditions/flux_bc.py index 99349a0f1..616425bed 100644 --- a/src/festim/boundary_conditions/flux_bc.py +++ b/src/festim/boundary_conditions/flux_bc.py @@ -34,100 +34,7 @@ def __init__(self, subdomain, value): self.subdomain = subdomain self.value = value - self.value_fenics = None - self.bc_expr = None - - @property - def value_fenics(self): - return self._value_fenics - - @value_fenics.setter - def value_fenics(self, value): - if value is None: - self._value_fenics = value - return - if not isinstance( - value, (fem.Function, fem.Constant, np.ndarray, ufl.core.expr.Expr) - ): - raise TypeError( - f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, np.ndarray or ufl.core.expr.Expr not {type(value)}" - ) - self._value_fenics = value - - @property - def time_dependent(self): - if self.value is None: - raise TypeError("Value must be given to determine if its time dependent") - if isinstance(self.value, fem.Constant): - return False - if callable(self.value): - arguments = self.value.__code__.co_varnames - return "t" in arguments - else: - return False - - @property - def temperature_dependent(self): - if self.value is None: - return False - if isinstance(self.value, fem.Constant): - return False - if callable(self.value): - arguments = self.value.__code__.co_varnames - return "T" in arguments - else: - return False - - 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 - 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 ufl Expression - - Args: - mesh (dolfinx.mesh.Mesh) : the mesh - 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 - if not isinstance(self.value(t=float(t)), (float, int)): - raise ValueError( - f"self.value should return a float or an int, not {type(self.value(t=float(t)))} " - ) - self.value_fenics = F.as_fenics_constant( - mesh=mesh, value=self.value(t=float(t)) - ) - else: - kwargs = {} - if "t" in arguments: - kwargs["t"] = t - if "x" in arguments: - kwargs["x"] = x - if "T" in arguments: - kwargs["T"] = temperature - - self.value_fenics = self.value(**kwargs) - - def update(self, t): - """Updates the flux bc value - - Args: - t (float): the time - """ - if callable(self.value): - arguments = self.value.__code__.co_varnames - if isinstance(self.value_fenics, fem.Constant) and "t" in arguments: - self.value_fenics.value = self.value(t=t) + self.value_fenics = F.ConvertToFenicsObject(value) class ParticleFluxBC(FluxBCBase): @@ -176,49 +83,6 @@ def __init__(self, subdomain, value, species, species_dependent_value={}): self.species = species self.species_dependent_value = species_dependent_value - 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 - 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 ufl Expression - - Args: - mesh (dolfinx.mesh.Mesh) : the mesh - 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 - if not isinstance(self.value(t=float(t)), (float, int)): - raise ValueError( - f"self.value should return a float or an int, not {type(self.value(t=float(t)))} " - ) - self.value_fenics = F.as_fenics_constant( - mesh=mesh, value=self.value(t=float(t)) - ) - else: - kwargs = {} - if "t" in arguments: - kwargs["t"] = t - if "x" in arguments: - kwargs["x"] = x - if "T" in arguments: - kwargs["T"] = temperature - - for name, species in self.species_dependent_value.items(): - kwargs[name] = species.concentration - - self.value_fenics = self.value(**kwargs) - class HeatFluxBC(FluxBCBase): """ diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 70448dcb3..2c90f6fb0 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -1,5 +1,8 @@ import dolfinx from dolfinx import fem +import numpy as np +from collections.abc import Callable +import ufl def as_fenics_constant( @@ -25,3 +28,178 @@ def as_fenics_constant( raise TypeError( f"Value must be a float, an int or a dolfinx.Constant, not {type(value)}" ) + + +class ConvertToFenicsObject: + """ + u = value + + Args: + value: The value of the user input + + Attributes: + value: The value of the user input + fenics_object: The value of the user input in fenics format + fenics_interpolation_expression: The expression of the user input that is used to + update the `fenics_object` + + """ + + input_value: ( + np.ndarray + | fem.Constant + | int + | float + | Callable[[np.ndarray], np.ndarray] + | Callable[[np.ndarray, float], np.ndarray] + | Callable[[float], float] + ) + + def __init__(self, input_value): + self, + self.input_value = input_value + + self.mapped_function = None + self.fenics_interpolation_expression = None + self.fenics_object = None + + @property + def time_dependent(self) -> bool: + """Returns true if the value given is time dependent""" + if self.input_value is None: + return False + if isinstance(self.input_value, fem.Constant): + return False + if callable(self.input_value): + arguments = self.input_value.__code__.co_varnames + return "t" in arguments + else: + return False + + @property + def temperature_dependent(self) -> bool: + """Returns true if the value given is temperature dependent""" + if self.input_value is None: + return False + if isinstance(self.input_value, fem.Constant): + return False + if callable(self.input_value): + arguments = self.input_value.__code__.co_varnames + return "T" in arguments + else: + return False + + def as_fenics_constant(self, value, mesh: dolfinx.mesh.Mesh) -> fem.Constant: + """Converts a value to a dolfinx.Constant + + Args: + value: the value to convert + mesh: the mesh of the domiain + + Raises: + TypeError: if the value is not a float, an int or a dolfinx.Constant + """ + if isinstance(value, (float, int)): + self.fenics_object = fem.Constant( + mesh, dolfinx.default_scalar_type(float(value)) + ) + else: + raise TypeError( + f"Value must be a float, an int or a dolfinx.Constant, not {type(value)}" + ) + + def as_mapped_function(self, mesh=None, t=None, temperature=None): + arguments = self.input_value.__code__.co_varnames + + kwargs = {} + if "t" in arguments: + kwargs["t"] = t + if "x" in arguments: + x = ufl.SpatialCoordinate(mesh) + kwargs["x"] = x + if "T" in arguments: + kwargs["T"] = temperature + + self.fenics_object = self.input_value(**kwargs) + + def as_fenics_interpolation_expression( + self, function_space, temperature=None, t=None, mesh=None + ): + + self.as_mapped_function(mesh=mesh, t=t, temperature=temperature) + # store the expression + self.fenics_interpolation_expression = fem.Expression( + self.fenics_object, + function_space.element.interpolation_points(), + ) + + def as_fenics_function(self, function_space, t, mesh=None, temperature=None): + + if self.fenics_interpolation_expression is None: + self.as_fenics_interpolation_expression( + mesh=mesh, function_space=function_space, temperature=temperature, t=t + ) + + self.fenics_object = fem.Function(function_space) + self.fenics_object.interpolate(self.fenics_interpolation_expression) + + def convert_value(self, mesh=None, function_space=None, t=None, temperature=None): + """Converts the value to a fenics object + + Args: + function_space: the function space of the fenics object + t: the time + temperature: the temperature + """ + if isinstance(self.input_value, fem.Constant): + self.fenics_object = self.input_value + elif isinstance(self.input_value, fem.Expression): + self.fenics_interpolation_expression = self.input_value + elif isinstance(self.input_value, (fem.Function, ufl.core.expr.Expr)): + self.fenics_object = self.input_value + elif isinstance(self.input_value, (float, int)): + if mesh is None: + raise ValueError("Mesh must be provided to create a constant") + self.as_fenics_constant(mesh=mesh, value=self.input_value) + elif callable(self.input_value): + args = self.input_value.__code__.co_varnames + # if only t is an argument, create constant object + if "t" in args and "x" not in args and "T" not in args: + if not isinstance(self.input_value(t=float(t)), (float, int)): + raise ValueError( + "self.value should return a float or an int, not " + + f"{type(self.input_value(t=float(t)))} " + ) + self.as_fenics_constant(mesh=mesh, value=self.input_value(t=float(t))) + elif self.temperature_dependent: + self.as_fenics_function( + mesh=mesh, + function_space=function_space, + t=t, + temperature=temperature, + ) + else: + self.as_fenics_function(mesh=mesh, function_space=function_space, t=t) + else: + raise TypeError( + f"Value must be a float, an int or a callable, not {type(self.input_value)}" + ) + + def update(self, t): + """Updates the value + + Args: + t (float): the time + """ + if callable(self.input_value): + arguments = self.input_value.__code__.co_varnames + if isinstance(self.fenics_object, fem.Constant) and "t" in arguments: + self.fenics_object.value = float(self.input_value(t=t)) + elif isinstance(self.fenics_object, fem.Function): + if self.fenics_interpolation_expression is not None: + self.fenics_object.interpolate(self.fenics_interpolation_expression) + + # if isinstance(self.fenics_object, fem.Constant): + # self.fenics_object.value = self.input_value(t=t) + # elif isinstance(self.fenics_object, fem.Function): + # self.fenics_object.interpolate(self.fenics_interpolation_expression) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 09e34f9a4..158557748 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -9,6 +9,7 @@ import ufl from dolfinx import fem from scifem import NewtonSolver +import numpy as np import festim.boundary_conditions import festim.problem @@ -178,6 +179,9 @@ def temperature(self, value): elif isinstance(value, (float, int, fem.Constant, fem.Function)): self._temperature = value elif callable(value): + arguments = value.__code__.co_varnames + if "T" in arguments: + raise ValueError("Temperature cannot be a function of temperature, T") self._temperature = value else: raise TypeError( @@ -195,9 +199,10 @@ def temperature_fenics(self, value): return elif not isinstance( value, - (fem.Constant, fem.Function), + festim.ConvertToFenicsObject, ): raise TypeError("Value must be a fem.Constant or fem.Function") + self._temperature_fenics = value @property @@ -312,54 +317,24 @@ def define_temperature(self): if self.temperature is None: raise ValueError("the temperature attribute needs to be defined") - # if temperature is a float or int, create a fem.Constant - elif isinstance(self.temperature, (float, int)): - self.temperature_fenics = as_fenics_constant( - self.temperature, self.mesh.mesh - ) - # if temperature is a fem.Constant or function, pass it to temperature_fenics - elif isinstance(self.temperature, (fem.Constant, fem.Function)): - self.temperature_fenics = self.temperature + self.temperature_fenics = festim.ConvertToFenicsObject( + input_value=self.temperature + ) - # if temperature is callable, process accordingly - elif callable(self.temperature): - arguments = self.temperature.__code__.co_varnames - if "t" in arguments and "x" not in arguments: - if not isinstance(self.temperature(t=float(self.t)), (float, int)): - raise ValueError( - f"self.temperature should return a float or an int, not " - f"{type(self.temperature(t=float(self.t)))} " - ) - # only t is an argument - self.temperature_fenics = as_fenics_constant( - mesh=self.mesh.mesh, value=self.temperature(t=float(self.t)) - ) - else: - x = ufl.SpatialCoordinate(self.mesh.mesh) - degree = 1 - element_temperature = basix.ufl.element( - basix.ElementFamily.P, - self.mesh.mesh.basix_cell(), - degree, - basix.LagrangeVariant.equispaced, - ) - function_space_temperature = fem.functionspace( - self.mesh.mesh, element_temperature - ) - self.temperature_fenics = fem.Function(function_space_temperature) - kwargs = {} - if "t" in arguments: - kwargs["t"] = self.t - if "x" in arguments: - kwargs["x"] = x - - # store the expression of the temperature - # to update the temperature_fenics later - self.temperature_expr = fem.Expression( - self.temperature(**kwargs), - function_space_temperature.element.interpolation_points(), - ) - self.temperature_fenics.interpolate(self.temperature_expr) + degree = 1 + element_temperature = basix.ufl.element( + basix.ElementFamily.P, + self.mesh.mesh.basix_cell(), + degree, + basix.LagrangeVariant.equispaced, + ) + function_space_temperature = fem.functionspace( + self.mesh.mesh, element_temperature + ) + + self.temperature_fenics.convert_value( + mesh=self.mesh.mesh, function_space=function_space_temperature, t=self.t + ) def initialise_exports(self): """Defines the export writers of the model, if field is given as @@ -443,7 +418,9 @@ def define_D_global(self, species): D = fem.Function(self.V_DG_1) expr = D_0 * ufl.exp( - -E_D / as_fenics_constant(k_B, self.mesh.mesh) / self.temperature_fenics + -E_D + / as_fenics_constant(k_B, self.mesh.mesh) + / self.temperature_fenics.fenics_object ) D_expr = fem.Expression(expr, self.V_DG_1.element.interpolation_points()) D.interpolate(D_expr) @@ -538,11 +515,20 @@ def define_boundary_conditions(self): # 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, - ) + if isinstance(bc.value, (int, float)): + bc.value_fenics.as_fenics_constant( + mesh=self.mesh.mesh, + value=bc.value, + ) + elif isinstance(bc.value, (fem.Function, ufl.core.expr.Expr)): + bc.value_fenics.fenics_object = bc.value + + else: + bc.value_fenics.as_mapped_function( + mesh=self.mesh.mesh, + temperature=self.temperature_fenics.fenics_object, + t=self.t, + ) super().define_boundary_conditions() @@ -562,14 +548,17 @@ def create_dirichletbc_form(self, bc): else: function_space_value = bc.species.collapsed_function_space - bc.create_value( - temperature=self.temperature_fenics, + bc.value_fenics.convert_value( + mesh=self.mesh.mesh, + temperature=self.temperature_fenics.fenics_object, function_space=function_space_value, t=self.t, ) # get dofs - if self.multispecies and isinstance(bc.value_fenics, (fem.Function)): + if self.multispecies and isinstance( + bc.value_fenics.fenics_object, (fem.Function) + ): function_space_dofs = ( bc.species.sub_function_space, bc.species.collapsed_function_space, @@ -583,14 +572,16 @@ def create_dirichletbc_form(self, bc): ) # create form - if not self.multispecies and isinstance(bc.value_fenics, (fem.Function)): + if not self.multispecies and isinstance( + bc.value_fenics.fenics_object, (fem.Function) + ): # no need to pass the functionspace since value_fenics is already a Function function_space_form = None else: function_space_form = bc.species.sub_function_space form = fem.dirichletbc( - value=bc.value_fenics, + value=bc.value_fenics.fenics_object, dofs=bc_dofs, V=function_space_form, ) @@ -602,18 +593,27 @@ def create_source_values_fenics(self): for source in self.sources: # create value_fenics for all F.ParticleSource objects if isinstance(source, festim.source.ParticleSource): - source.create_value_fenics( - mesh=self.mesh.mesh, - temperature=self.temperature_fenics, - t=self.t, - ) + if isinstance(source.value, (int, float)): + source.value_fenics.as_fenics_constant( + mesh=self.mesh.mesh, + value=source.value, + ) + elif isinstance(source.value, (fem.Function, ufl.core.expr.Expr)): + source.value_fenics.fenics_object = self.value + elif callable(source.value): + source.value_fenics.as_mapped_function( + mesh=self.mesh.mesh, + temperature=self.temperature_fenics.fenics_object, + 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, boundary_conditions.ParticleFluxBC): - bc.create_value_fenics( + + bc.value_fenics.convert_value( mesh=self.mesh.mesh, temperature=self.temperature_fenics, t=self.t, @@ -640,18 +640,27 @@ def create_initial_conditions(self): else: function_space_value = condition.species.collapsed_function_space - condition.create_expr_fenics( - mesh=self.mesh.mesh, - temperature=self.temperature_fenics, - function_space=function_space_value, - ) + if isinstance(condition.value, (int, float)): + condition.value_fenics.fenics_interpolation_expression = ( + lambda x: np.full(x.shape[1], condition.value) + ) + else: + condition.value_fenics.as_fenics_interpolation_expression( + mesh=self.mesh.mesh, + temperature=self.temperature_fenics.fenics_object, + function_space=function_space_value, + ) # assign to previous solution of species if not self.multispecies: - condition.species.prev_solution.interpolate(condition.expr_fenics) + condition.species.prev_solution.interpolate( + condition.value_fenics.fenics_interpolation_expression + ) else: idx = self.species.index(condition.species) - self.u_n.sub(idx).interpolate(condition.expr_fenics) + self.u_n.sub(idx).interpolate( + condition.value_fenics.fenics_interpolation_expression + ) def create_formulation(self): """Creates the formulation of the model""" @@ -666,7 +675,7 @@ def create_formulation(self): for vol in self.volume_subdomains: D = vol.material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature_fenics, spe + self.mesh.mesh, self.temperature_fenics.fenics_object, spe ) if spe.mobile: self.formulation += ufl.dot(D * ufl.grad(u), ufl.grad(v)) * self.dx( @@ -680,7 +689,7 @@ def create_formulation(self): for reactant in reaction.reactant: if isinstance(reactant, festim.species.Species): self.formulation += ( - reaction.reaction_term(self.temperature_fenics) + reaction.reaction_term(self.temperature_fenics.fenics_object) * reactant.test_function * self.dx(reaction.volume.id) ) @@ -692,14 +701,14 @@ def create_formulation(self): products = [reaction.product] for product in products: self.formulation += ( - -reaction.reaction_term(self.temperature_fenics) + -reaction.reaction_term(self.temperature_fenics.fenics_object) * product.test_function * self.dx(reaction.volume.id) ) # add sources for source in self.sources: self.formulation -= ( - source.value_fenics + source.value_fenics.fenics_object * source.species.test_function * self.dx(source.volume.id) ) @@ -708,14 +717,14 @@ def create_formulation(self): for bc in self.boundary_conditions: if isinstance(bc, boundary_conditions.ParticleFluxBC): self.formulation -= ( - bc.value_fenics + bc.value_fenics.fenics_object * 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.value_fenics.fenics_object * flux_bc.species.test_function * self.ds(flux_bc.subdomain.id) ) @@ -751,10 +760,7 @@ def update_time_dependent_values(self): if not self.temperature_time_dependent: return - if isinstance(self.temperature_fenics, fem.Constant): - self.temperature_fenics.value = self.temperature(t=t) - elif isinstance(self.temperature_fenics, fem.Function): - self.temperature_fenics.interpolate(self.temperature_expr) + self.temperature_fenics.update(t=t) for bc in self.boundary_conditions: if isinstance( @@ -764,12 +770,12 @@ def update_time_dependent_values(self): boundary_conditions.ParticleFluxBC, ), ): - if bc.temperature_dependent: - bc.update(t=t) + if bc.value_fenics.temperature_dependent: + bc.value_fenics.update(t=t) for source in self.sources: - if source.temperature_dependent: - source.update(t=t) + if source.value_fenics.temperature_dependent: + source.value_fenics.update(t=t) def post_processing(self): """Post processes the model""" diff --git a/src/festim/initial_condition.py b/src/festim/initial_condition.py index b4947c521..ff9c621fb 100644 --- a/src/festim/initial_condition.py +++ b/src/festim/initial_condition.py @@ -1,6 +1,7 @@ import numpy as np import ufl from dolfinx import fem +import festim as F # TODO rename this to InitialConcentration and create a new base class @@ -30,69 +31,11 @@ def __init__(self, value, species): self.value = value self.species = species - self.expr_fenics = None - - def create_expr_fenics(self, mesh, temperature, function_space): - """Creates the expr_fenics of the initial condition. - If the value is a float or int, a function is created with an array with - the shape of the mesh and all set to the value. - Otherwise, it is converted to a fem.Expression. - - Args: - mesh (dolfinx.mesh.Mesh) : the mesh - temperature (float): the temperature - function_space(dolfinx.fem.FunctionSpaceBase): the function space of the species - """ - x = ufl.SpatialCoordinate(mesh) - - if isinstance(self.value, (int, float)): - self.expr_fenics = lambda x: np.full(x.shape[1], self.value) - - elif callable(self.value): - arguments = self.value.__code__.co_varnames - kwargs = {} - if "t" in arguments: - raise ValueError("Initial condition cannot be a function of time.") - if "x" in arguments: - kwargs["x"] = x - if "T" in arguments: - kwargs["T"] = temperature - - self.expr_fenics = fem.Expression( - self.value(**kwargs), - function_space.element.interpolation_points(), - ) + self.value_fenics = F.ConvertToFenicsObject(value) class InitialTemperature: def __init__(self, value) -> None: self.value = value - self.expr_fenics = None - - def create_expr_fenics(self, mesh, function_space): - """Creates the expr_fenics of the initial condition. - If the value is a float or int, a function is created with an array with - the shape of the mesh and all set to the value. - Otherwise, it is converted to a fem.Expression. - - Args: - mesh (dolfinx.mesh.Mesh) : the mesh - function_space(dolfinx.fem.FunctionSpace): the function space of the species - """ - x = ufl.SpatialCoordinate(mesh) - - if isinstance(self.value, (int, float)): - self.expr_fenics = lambda x: np.full(x.shape[1], self.value) - - elif callable(self.value): - arguments = self.value.__code__.co_varnames - kwargs = {} - if "t" in arguments: - raise ValueError("Initial condition cannot be a function of time.") - if "x" in arguments: - kwargs["x"] = x - self.expr_fenics = fem.Expression( - self.value(**kwargs), - function_space.element.interpolation_points(), - ) + self.value_fenics = F.ConvertToFenicsObject(value) diff --git a/src/festim/problem.py b/src/festim/problem.py index d1066ef79..cfff53269 100644 --- a/src/festim/problem.py +++ b/src/festim/problem.py @@ -189,9 +189,9 @@ def iterate(self): def update_time_dependent_values(self): t = float(self.t) for bc in self.boundary_conditions: - if bc.time_dependent: - bc.update(t=t) + if bc.value_fenics.time_dependent: + bc.value_fenics.update(t=t) for source in self.sources: - if source.time_dependent: - source.update(t=t) + if source.value_fenics.time_dependent: + source.value_fenics.update(t=t) diff --git a/src/festim/source.py b/src/festim/source.py index a792dda49..153273bd9 100644 --- a/src/festim/source.py +++ b/src/festim/source.py @@ -40,8 +40,7 @@ def __init__(self, value, volume): self.value = value self.volume = volume - self.value_fenics = None - self.source_expr = None + self.value_fenics = F.ConvertToFenicsObject(value) @property def volume(self): @@ -54,46 +53,6 @@ def volume(self, value): raise TypeError("volume must be of type festim.VolumeSubdomain") self._volume = value - @property - def value_fenics(self): - return self._value_fenics - - @value_fenics.setter - def value_fenics(self, value): - if value is None: - self._value_fenics = value - return - if not isinstance( - value, (fem.Function, fem.Constant, np.ndarray, ufl.core.expr.Expr) - ): - raise TypeError( - f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, np.ndarray or a ufl.core.expr.Expr, not {type(value)}" - ) - self._value_fenics = value - - @property - def time_dependent(self): - if self.value is None: - return False - if isinstance(self.value, fem.Constant): - return False - if callable(self.value): - arguments = self.value.__code__.co_varnames - return "t" in arguments - else: - return False - - def update(self, t): - """Updates the source value - - Args: - t (float): the time - """ - if callable(self.value): - arguments = self.value.__code__.co_varnames - if isinstance(self.value_fenics, fem.Constant) and "t" in arguments: - self.value_fenics.value = self.value(t=t) - class ParticleSource(SourceBase): def __init__(self, value, volume, species): @@ -112,107 +71,7 @@ def species(self, value): self._species = value - @property - def temperature_dependent(self): - if self.value is None: - return False - if isinstance(self.value, fem.Constant): - return False - if callable(self.value): - arguments = self.value.__code__.co_varnames - return "T" in arguments - else: - return False - - def create_value_fenics(self, mesh, temperature, t: fem.Constant): - """Creates the value of the source 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 ufl Expression - - Args: - mesh (dolfinx.mesh.Mesh) : the mesh - 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 isinstance(self.value, (fem.Function, ufl.core.expr.Expr)): - self.value_fenics = 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 - if not isinstance(self.value(t=float(t)), (float, int)): - raise ValueError( - f"self.value should return a float or an int, not {type(self.value(t=float(t)))} " - ) - self.value_fenics = F.as_fenics_constant( - mesh=mesh, value=self.value(t=float(t)) - ) - else: - kwargs = {} - if "t" in arguments: - kwargs["t"] = t - if "x" in arguments: - kwargs["x"] = x - if "T" in arguments: - kwargs["T"] = temperature - - self.value_fenics = self.value(**kwargs) - class HeatSource(SourceBase): def __init__(self, value, volume): super().__init__(value, volume) - - def create_value_fenics( - self, - mesh: dolfinx.mesh.Mesh, - t: fem.Constant, - ): - """Creates the value of the source 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 ufl.Expression - - Args: - mesh (dolfinx.mesh.Mesh) : the mesh - 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 isinstance(self.value, (fem.Function, ufl.core.expr.Expr)): - self.value_fenics = self.value - - elif callable(self.value): - arguments = self.value.__code__.co_varnames - - if "t" in arguments and "x" not in arguments: - # only t is an argument - if not isinstance(self.value(t=float(t)), (float, int)): - raise ValueError( - f"self.value should return a float or an int, not {type(self.value(t=float(t)))} " - ) - self.value_fenics = F.as_fenics_constant( - mesh=mesh, value=self.value(t=float(t)) - ) - else: - kwargs = {} - if "t" in arguments: - kwargs["t"] = t - if "x" in arguments: - kwargs["x"] = x - # TODO could the source be dependend on T? why not? - - self.value_fenics = self.value(**kwargs) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index b6639d8a4..1c0887643 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -107,7 +107,7 @@ def test_define_temperature(input, expected_type): my_model.define_temperature() # TEST - assert isinstance(my_model.temperature_fenics, expected_type) + assert isinstance(my_model.temperature_fenics.fenics_object, expected_type) @pytest.mark.parametrize( @@ -129,7 +129,7 @@ def test_define_temperature_error_if_ufl_conditional_t_only(input): with pytest.raises( ValueError, - match="self.temperature should return a float or an int, not ", + match="self.value should return a float or an int, not", ): my_model.define_temperature() @@ -258,6 +258,7 @@ def test_define_D_global_different_temperatures(): my_model.define_function_spaces() my_model.define_meshtags_and_measures() + my_model.t = 1 my_model.define_temperature() D_computed, D_expr = my_model.define_D_global(H) @@ -293,6 +294,7 @@ def test_define_D_global_different_materials(): my_model.define_function_spaces() my_model.define_meshtags_and_measures() + my_model.t = 0 my_model.define_temperature() D_computed, D_expr = my_model.define_D_global(H) @@ -340,6 +342,7 @@ def test_initialise_exports_multiple_exports_same_species(): my_model.define_function_spaces() my_model.define_meshtags_and_measures() + my_model.t = 0 my_model.define_temperature() my_model.initialise_exports() @@ -402,6 +405,7 @@ def test_define_D_global_multispecies(): my_model.define_function_spaces() my_model.define_meshtags_and_measures() + my_model.t = 0 my_model.define_temperature() D_A_computed, D_A_expr = my_model.define_D_global(A) @@ -677,8 +681,8 @@ def test_update_time_dependent_values_source(source_value, expected_values): my_model.update_time_dependent_values() # TEST - if isinstance(my_model.sources[0].value_fenics, fem.Constant): - computed_value = float(my_model.sources[0].value_fenics) + if isinstance(my_model.sources[0].value_fenics.fenics_object, fem.Constant): + computed_value = float(my_model.sources[0].value_fenics.fenics_object) assert np.isclose(computed_value, expected_values[i]) @@ -731,8 +735,8 @@ def test_update_sources_with_time_dependent_temperature( my_model.update_time_dependent_values() # TEST - if isinstance(my_model.sources[0].value_fenics, fem.Constant): - computed_value = float(my_model.sources[0].value_fenics) + if isinstance(my_model.sources[0].value_fenics.fenics_object, fem.Constant): + computed_value = float(my_model.sources[0].value_fenics.fenics_object) assert np.isclose(computed_value, expected_values[i]) @@ -763,8 +767,8 @@ def test_create_source_values_fenics_multispecies(): my_model.create_source_values_fenics() # TEST - assert np.isclose(my_model.sources[0].value_fenics.value, 5) - assert np.isclose(my_model.sources[1].value_fenics.value, 11) + assert np.isclose(float(my_model.sources[0].value_fenics.fenics_object), 5) + assert np.isclose(float(my_model.sources[1].value_fenics.fenics_object), 11) # TODO replace this by a proper MMS test @@ -1229,5 +1233,9 @@ def test_create_flux_values_fenics_multispecies(): my_model.create_flux_values_fenics() # TEST - assert np.isclose(my_model.boundary_conditions[0].value_fenics.value, 5) - assert np.isclose(my_model.boundary_conditions[1].value_fenics.value, 11) + assert np.isclose( + float(my_model.boundary_conditions[0].value_fenics.fenics_object), 5 + ) + assert np.isclose( + float(my_model.boundary_conditions[1].value_fenics.fenics_object), 11 + ) From ac60724213483c355b4d7769bedecf8954ede1c0 Mon Sep 17 00:00:00 2001 From: jhdark Date: Wed, 22 Jan 2025 17:10:38 -0500 Subject: [PATCH 02/40] Rename to Value, work with temperature only classes --- src/festim/__init__.py | 8 +- .../boundary_conditions/dirichlet_bc.py | 205 ++++---------- src/festim/boundary_conditions/flux_bc.py | 22 +- src/festim/heat_transfer_problem.py | 66 +++-- src/festim/helpers.py | 254 ++++++++++++------ src/festim/hydrogen_transport_problem.py | 138 +++++----- src/festim/initial_condition.py | 34 ++- src/festim/problem.py | 8 +- src/festim/source.py | 23 +- test/test_h_transport_problem.py | 56 ++-- test/test_heat_transfer_problem.py | 7 +- 11 files changed, 447 insertions(+), 374 deletions(-) diff --git a/src/festim/__init__.py b/src/festim/__init__.py index 5edea51b3..d0f4c0f58 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -36,7 +36,13 @@ from .exports.vtx import VTXSpeciesExport, VTXTemperatureExport from .exports.xdmf import XDMFExport from .heat_transfer_problem import HeatTransferProblem -from .helpers import as_fenics_constant, ConvertToFenicsObject +from .helpers import ( + as_fenics_constant, + as_mapped_function, + as_fenics_interpolation_expression, + as_fenics_interp_expr_and_function, + Value, +) from .hydrogen_transport_problem import ( HTransportProblemDiscontinuous, HydrogenTransportProblem, diff --git a/src/festim/boundary_conditions/dirichlet_bc.py b/src/festim/boundary_conditions/dirichlet_bc.py index 67d8dc670..80e834363 100644 --- a/src/festim/boundary_conditions/dirichlet_bc.py +++ b/src/festim/boundary_conditions/dirichlet_bc.py @@ -49,38 +49,22 @@ def __init__( self.subdomain = subdomain self.value = value - self.value_fenics = helpers.ConvertToFenicsObject(value) - # self.bc_expr = None - - # @property - # def value_fenics(self): - # return self._value_fenics - - # @value_fenics.setter - # def value_fenics(self, value: None | fem.Function | fem.Constant | np.ndarray): - # if value is None: - # self._value_fenics = value - # return - # if not isinstance(value, (fem.Function, fem.Constant, np.ndarray)): - # # FIXME: Should we allow sending in a callable here? - # raise TypeError( - # "Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not" - # + f"{type(value)}" - # ) - # self._value_fenics = value - - # @property - # def time_dependent(self) -> bool: - # """Returns true if the value of the boundary condition is time dependent""" - # if self.value is None: - # return False - # if isinstance(self.value, fem.Constant): - # return False - # if callable(self.value): - # arguments = self.value.__code__.co_varnames - # return "t" in arguments - # else: - # return False + @property + def value(self): + return self._value + + @value.setter + def value(self, value): + if value is None: + self._value = value + elif isinstance(value, (float, int, fem.Constant, fem.Function)): + self._value = helpers.Value(value) + elif callable(value): + self._value = helpers.Value(value) + else: + raise TypeError( + "Value must be a float, int, fem.Constant, fem.Function, or callable" + ) def define_surface_subdomain_dofs( self, @@ -120,19 +104,6 @@ def define_surface_subdomain_dofs( return bc_dofs - # def update(self, t: float): - # """Updates the boundary condition value - - # Args: - # t (float): the time - # """ - # if callable(self.value): - # arguments = self.value.__code__.co_varnames - # if isinstance(self.value_fenics, fem.Constant) and "t" in arguments: - # self.value_fenics.value = self.value(t=t) - # else: - # self.value_fenics.interpolate(self.bc_expr) - class FixedConcentrationBC(DirichletBCBase): """ @@ -174,122 +145,42 @@ def __init__( self.species = species super().__init__(subdomain, value) - # @property - # def temperature_dependent(self): - # if self.value is None: - # return False - # if isinstance(self.value, fem.Constant): - # return False - # if callable(self.value): - # arguments = self.value.__code__.co_varnames - # return "T" in arguments - # else: - # return False - - # def create_value( - # self, - # function_space: fem.FunctionSpace, - # temperature: float | fem.Constant, - # t: float | 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 `dolfinx.fem.Constant`. - # If the value is a function of t, it is converted to `dolfinx.fem.Constant`. - # Otherwise, it is converted to a `dolfinx.fem.Function`.Function and the - # expression of the function is stored in `bc_expr`. - - # Args: - # function_space (dolfinx.fem.FunctionSpace): the function space - # temperature: The temperature - # t (dolfinx.fem.Constant): the time - # """ - # mesh = function_space.mesh - # x = ufl.SpatialCoordinate(mesh) - - # if isinstance(self.value, (int, float)): - # self.value_fenics = helpers.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 - # if not isinstance(self.value(t=float(t)), (float, int)): - # raise ValueError( - # "self.value should return a float or an int, not " - # + f"{type(self.value(t=float(t)))} " - # ) - # self.value_fenics = helpers.as_fenics_constant( - # mesh=mesh, value=self.value(t=float(t)) - # ) - # else: - # self.value_fenics = fem.Function(function_space) - # kwargs = {} - # if "t" in arguments: - # kwargs["t"] = t - # if "x" in arguments: - # 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(), - # ) - # self.value_fenics.interpolate(self.bc_expr) - # alias for FixedConcentrationBC DirichletBC = FixedConcentrationBC class FixedTemperatureBC(DirichletBCBase): - def create_value(self, function_space: fem.FunctionSpace, 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 `dolfinx.fem.Constant`. - If the value is a function of t, it is converted to a `dolfinx.fem.Constant`. - Otherwise, it is converted to a` dolfinx.fem.Function` and the - expression of the function is stored in `bc_expr`. + """ + 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 - Args: - function_space: the function space - t: the time - """ - mesh = function_space.mesh - x = ufl.SpatialCoordinate(mesh) - - if isinstance(self.value, (int, float)): - self.value_fenics = helpers.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: - # only t is an argument - if not isinstance(self.value(t=float(t)), (float, int)): - raise ValueError( - "self.value should return a float or an int, not " - + f"{type(self.value(t=float(t)))} " - ) - self.value_fenics = helpers.as_fenics_constant( - mesh=mesh, value=self.value(t=float(t)) - ) - else: - self.value_fenics = fem.Function(function_space) - kwargs = {} - if "t" in arguments: - kwargs["t"] = t - if "x" in arguments: - kwargs["x"] = x - - # 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(), - ) - self.value_fenics.interpolate(self.bc_expr) + Examples: + + .. highlight:: python + .. code-block:: python + + from festim import FixedTemperatureBC + FixedTemperatureBC(subdomain=my_subdomain, value=1) + FixedTemperatureBC(subdomain=my_subdomain, + value=lambda x: 1 + x[0]) + FixedTemperatureBC(subdomain=my_subdomain, + value=lambda t: 1 + t) + FixedTemperatureBC(subdomain=my_subdomain, + value=lambda x, t: 1 + x[0] + t) + + """ + + def __init__( + self, + subdomain: _subdomain.SurfaceSubdomain, + value: np.ndarray | fem.Constant | int | float | Callable, + ): + super().__init__(subdomain, value) + + if self.value.temperature_dependent: + raise ValueError( + "Temperature dependent boundary conditions are not supported for FixedTemperatureBC" + ) diff --git a/src/festim/boundary_conditions/flux_bc.py b/src/festim/boundary_conditions/flux_bc.py index 616425bed..3a72c1fc0 100644 --- a/src/festim/boundary_conditions/flux_bc.py +++ b/src/festim/boundary_conditions/flux_bc.py @@ -1,5 +1,3 @@ -import numpy as np -import ufl from dolfinx import fem import festim as F @@ -34,7 +32,22 @@ def __init__(self, subdomain, value): self.subdomain = subdomain self.value = value - self.value_fenics = F.ConvertToFenicsObject(value) + @property + def value(self): + return self._value + + @value.setter + def value(self, value): + if value is None: + self._value = value + elif isinstance(value, (float, int, fem.Constant, fem.Function)): + self._value = F.Value(value) + elif callable(value): + self._value = F.Value(value) + else: + raise TypeError( + "Value must be a float, int, fem.Constant, fem.Function, or callable" + ) class ParticleFluxBC(FluxBCBase): @@ -118,3 +131,6 @@ class HeatFluxBC(FluxBCBase): def __init__(self, subdomain, value): super().__init__(subdomain=subdomain, value=value) + + if self.value.temperature_dependent: + raise ValueError("Heat flux cannot be temperature dependent") diff --git a/src/festim/heat_transfer_problem.py b/src/festim/heat_transfer_problem.py index 9c49fce6c..90b1207c0 100644 --- a/src/festim/heat_transfer_problem.py +++ b/src/festim/heat_transfer_problem.py @@ -3,6 +3,8 @@ from dolfinx import fem from dolfinx.io import VTXWriter +import numpy as np + from festim import boundary_conditions, exports, helpers, problem from festim import source as _source @@ -111,8 +113,9 @@ def create_dirichletbc_form(self, bc): dolfinx.fem.bcs.DirichletBC: A representation of the boundary condition for modifying linear systems. """ - bc.create_value( + bc.value.convert_input_value( function_space=self.function_space, + mesh=self.mesh.mesh, t=self.t, ) @@ -121,11 +124,11 @@ def create_dirichletbc_form(self, bc): function_space=self.function_space, ) - if isinstance(bc.value_fenics, (fem.Function)): - form = fem.dirichletbc(value=bc.value_fenics, dofs=bc_dofs) + if isinstance(bc.value.fenics_object, (fem.Function)): + form = fem.dirichletbc(value=bc.value.fenics_object, dofs=bc_dofs) else: form = fem.dirichletbc( - value=bc.value_fenics, dofs=bc_dofs, V=self.function_space + value=bc.value.fenics_object, dofs=bc_dofs, V=self.function_space ) return form @@ -134,19 +137,30 @@ def create_source_values_fenics(self): """For each source create the value_fenics""" for source in self.sources: # create value_fenics for all source objects - source.create_value_fenics( - mesh=self.mesh.mesh, - t=self.t, - ) + if isinstance( + source.value.input_value, + (int, float, fem.Constant, fem.Function, ufl.core.expr.Expr), + ): + source.value.convert_input_value( + mesh=self.mesh.mesh, + t=self.t, + ) + elif callable(source.value.input_value): + source.value.fenics_object = helpers.as_mapped_function( + value=source.value.input_value, + mesh=self.mesh.mesh, + t=self.t, + ) def create_flux_values_fenics(self): """For each heat flux create the value_fenics""" for bc in self.boundary_conditions: # create value_fenics for all F.HeatFluxBC objects if isinstance(bc, boundary_conditions.HeatFluxBC): - bc.create_value_fenics( + + bc.value.convert_input_value( + function_space=self.function_space, mesh=self.mesh.mesh, - temperature=self.u, t=self.t, ) @@ -164,13 +178,27 @@ def create_initial_conditions(self): # create value_fenics for condition - self.initial_condition.create_expr_fenics( - mesh=self.mesh.mesh, - function_space=self.function_space, - ) + # self.initial_condition.value.create_expr_fenics( + # mesh=self.mesh.mesh, + # function_space=self.function_space, + # ) + if isinstance(self.initial_condition.value.input_value, (int, float)): + self.initial_condition.value.fenics_interpolation_expression = ( + lambda x: np.full(x.shape[1], self.initial_condition.value.input_value) + ) + else: + self.initial_condition.value.fenics_interpolation_expression = ( + helpers.as_fenics_interpolation_expression( + value=self.initial_condition.value.input_value, + function_space=self.function_space, + mesh=self.mesh.mesh, + ) + ) # assign to previous solution of species - self.u_n.interpolate(self.initial_condition.expr_fenics) + self.u_n.interpolate( + self.initial_condition.value.fenics_interpolation_expression + ) def create_formulation(self): """Creates the formulation of the model""" @@ -205,14 +233,18 @@ def create_formulation(self): # add sources for source in self.sources: self.formulation -= ( - source.value_fenics * self.test_function * self.dx(source.volume.id) + source.value.fenics_object + * self.test_function + * self.dx(source.volume.id) ) # add fluxes for bc in self.boundary_conditions: if isinstance(bc, boundary_conditions.HeatFluxBC): self.formulation -= ( - bc.value_fenics * self.test_function * self.ds(bc.subdomain.id) + bc.value.fenics_object + * self.test_function + * self.ds(bc.subdomain.id) ) def initialise_exports(self): diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 2c90f6fb0..122edb7f7 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -5,6 +5,31 @@ import ufl +# def as_fenics_constant( +# value: float | int | fem.Constant, mesh: dolfinx.mesh.Mesh +# ) -> fem.Constant: +# """Converts a value to a dolfinx.Constant + +# Args: +# value: the value to convert +# mesh: the mesh of the domiain + +# Returns: +# The converted value + +# Raises: +# TypeError: if the value is not a float, an int or a dolfinx.Constant +# """ +# if isinstance(value, (float, int)): +# return fem.Constant(mesh, dolfinx.default_scalar_type(value)) +# elif isinstance(value, fem.Constant): +# return value +# else: +# raise TypeError( +# f"Value must be a float, an int or a dolfinx.Constant, not {type(value)}" +# ) + + def as_fenics_constant( value: float | int | fem.Constant, mesh: dolfinx.mesh.Mesh ) -> fem.Constant: @@ -15,13 +40,13 @@ def as_fenics_constant( mesh: the mesh of the domiain Returns: - The converted value + (fem.Constant) The converted value Raises: TypeError: if the value is not a float, an int or a dolfinx.Constant """ if isinstance(value, (float, int)): - return fem.Constant(mesh, dolfinx.default_scalar_type(value)) + return fem.Constant(mesh, dolfinx.default_scalar_type(float(value))) elif isinstance(value, fem.Constant): return value else: @@ -30,31 +55,99 @@ def as_fenics_constant( ) -class ConvertToFenicsObject: +def as_mapped_function( + value, mesh=None, t=None, temperature=None +) -> ufl.core.expr.Expr: + """Maps a user given function + + Args: + value (Callable): the callable to convert + mesh (dolfinx.mesh.Mesh): the mesh of the domain + t (fem.Constant): the time + temperature (fem.Function, fem.Constant or ufl.core.expr.Expr): the temperature + + Returns: + (ufl.core.expr.Expr) The mapped function """ - u = value + + arguments = value.__code__.co_varnames + + kwargs = {} + if "t" in arguments: + kwargs["t"] = t + if "x" in arguments: + x = ufl.SpatialCoordinate(mesh) + kwargs["x"] = x + if "T" in arguments: + kwargs["T"] = temperature + + return value(**kwargs) + + +def as_fenics_interpolation_expression( + value, function_space, mesh=None, t=None, temperature=None +) -> fem.Expression: + """Converts a callable input value to a fenics expression""" + + mapped_function = as_mapped_function( + value=value, mesh=mesh, t=t, temperature=temperature + ) + + return fem.Expression( + mapped_function, + function_space.element.interpolation_points(), + ) + + +def as_fenics_interp_expr_and_function( + value, function_space, mesh=None, t=None, temperature=None +) -> fem.Function: + + fenics_interpolation_expression = as_fenics_interpolation_expression( + value=value, + function_space=function_space, + mesh=mesh, + t=t, + temperature=temperature, + ) + fenics_object = fem.Function(function_space) + fenics_object.interpolate(fenics_interpolation_expression) + + return fenics_interpolation_expression, fenics_object + + +class Value: + """ + A class to handle input values from users and convert them to a relevent fenics object Args: - value: The value of the user input + input_value: The value of the user input Attributes: - value: The value of the user input - fenics_object: The value of the user input in fenics format - fenics_interpolation_expression: The expression of the user input that is used to + input_value : The value of the user input + fenics_interpolation_expression : The expression of the user input that is used to update the `fenics_object` + fenics_object : The value of the user input in fenics format """ input_value: ( - np.ndarray - | fem.Constant - | int + int | float + | np.ndarray | Callable[[np.ndarray], np.ndarray] | Callable[[np.ndarray, float], np.ndarray] | Callable[[float], float] + | fem.Constant + | fem.Expression + | ufl.core.expr.Expr + | fem.Function ) + mapped_function: ufl.core.expr.Expr + fenics_interpolation_expression: fem.Expression + fenics_object: fem.Function | fem.Constant | ufl.core.expr.Expr + def __init__(self, input_value): self, self.input_value = input_value @@ -63,6 +156,34 @@ def __init__(self, input_value): self.fenics_interpolation_expression = None self.fenics_object = None + @property + def input_value(self): + return self._input_value + + @input_value.setter + def input_value(self, value): + if value is None: + self._input_value = value + elif isinstance( + value, + ( + float, + int, + fem.Constant, + np.ndarray, + fem.Function, + ufl.core.expr.Expr, + fem.Function, + ), + ): + self._input_value = value + elif callable(value): + self._input_value = value + else: + raise TypeError( + "Value must be a float, int, fem.Constant, fem.Function, or callable" + ) + @property def time_dependent(self) -> bool: """Returns true if the value given is time dependent""" @@ -89,78 +210,30 @@ def temperature_dependent(self) -> bool: else: return False - def as_fenics_constant(self, value, mesh: dolfinx.mesh.Mesh) -> fem.Constant: - """Converts a value to a dolfinx.Constant - - Args: - value: the value to convert - mesh: the mesh of the domiain - - Raises: - TypeError: if the value is not a float, an int or a dolfinx.Constant - """ - if isinstance(value, (float, int)): - self.fenics_object = fem.Constant( - mesh, dolfinx.default_scalar_type(float(value)) - ) - else: - raise TypeError( - f"Value must be a float, an int or a dolfinx.Constant, not {type(value)}" - ) - - def as_mapped_function(self, mesh=None, t=None, temperature=None): - arguments = self.input_value.__code__.co_varnames - - kwargs = {} - if "t" in arguments: - kwargs["t"] = t - if "x" in arguments: - x = ufl.SpatialCoordinate(mesh) - kwargs["x"] = x - if "T" in arguments: - kwargs["T"] = temperature - - self.fenics_object = self.input_value(**kwargs) - - def as_fenics_interpolation_expression( - self, function_space, temperature=None, t=None, mesh=None + def convert_input_value( + self, mesh=None, function_space=None, t=None, temperature=None ): - - self.as_mapped_function(mesh=mesh, t=t, temperature=temperature) - # store the expression - self.fenics_interpolation_expression = fem.Expression( - self.fenics_object, - function_space.element.interpolation_points(), - ) - - def as_fenics_function(self, function_space, t, mesh=None, temperature=None): - - if self.fenics_interpolation_expression is None: - self.as_fenics_interpolation_expression( - mesh=mesh, function_space=function_space, temperature=temperature, t=t - ) - - self.fenics_object = fem.Function(function_space) - self.fenics_object.interpolate(self.fenics_interpolation_expression) - - def convert_value(self, mesh=None, function_space=None, t=None, temperature=None): - """Converts the value to a fenics object + """Converts a user given value to a relevent fenics object depending + on the type of the value provided Args: - function_space: the function space of the fenics object - t: the time - temperature: the temperature + mesh (dolfinx.mesh.Mesh): the mesh of the domain + function_space (dolfinx.fem.function.FunctionSpace): the function space of the fenics object + t (fem.Constant): the time + temperature (fem.Function, fem.Constant or ufl.core.expr.Expr): the temperature """ if isinstance(self.input_value, fem.Constant): self.fenics_object = self.input_value + elif isinstance(self.input_value, fem.Expression): self.fenics_interpolation_expression = self.input_value + elif isinstance(self.input_value, (fem.Function, ufl.core.expr.Expr)): self.fenics_object = self.input_value + elif isinstance(self.input_value, (float, int)): - if mesh is None: - raise ValueError("Mesh must be provided to create a constant") - self.as_fenics_constant(mesh=mesh, value=self.input_value) + self.fenics_object = as_fenics_constant(value=self.input_value, mesh=mesh) + elif callable(self.input_value): args = self.input_value.__code__.co_varnames # if only t is an argument, create constant object @@ -170,16 +243,32 @@ def convert_value(self, mesh=None, function_space=None, t=None, temperature=None "self.value should return a float or an int, not " + f"{type(self.input_value(t=float(t)))} " ) - self.as_fenics_constant(mesh=mesh, value=self.input_value(t=float(t))) + + self.fenics_object = as_fenics_constant( + value=self.input_value(t=float(t)), mesh=mesh + ) + elif self.temperature_dependent: - self.as_fenics_function( - mesh=mesh, - function_space=function_space, - t=t, - temperature=temperature, + self.fenics_interpolation_expression, self.fenics_object = ( + as_fenics_interp_expr_and_function( + value=self.input_value, + function_space=function_space, + mesh=mesh, + t=t, + temperature=temperature, + ) ) + else: - self.as_fenics_function(mesh=mesh, function_space=function_space, t=t) + self.fenics_interpolation_expression, self.fenics_object = ( + as_fenics_interp_expr_and_function( + value=self.input_value, + function_space=function_space, + mesh=mesh, + t=t, + ) + ) + else: raise TypeError( f"Value must be a float, an int or a callable, not {type(self.input_value)}" @@ -193,13 +282,10 @@ def update(self, t): """ if callable(self.input_value): arguments = self.input_value.__code__.co_varnames + if isinstance(self.fenics_object, fem.Constant) and "t" in arguments: self.fenics_object.value = float(self.input_value(t=t)) + elif isinstance(self.fenics_object, fem.Function): if self.fenics_interpolation_expression is not None: self.fenics_object.interpolate(self.fenics_interpolation_expression) - - # if isinstance(self.fenics_object, fem.Constant): - # self.fenics_object.value = self.input_value(t=t) - # elif isinstance(self.fenics_object, fem.Function): - # self.fenics_object.interpolate(self.fenics_interpolation_expression) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 158557748..b50575284 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -165,7 +165,6 @@ def __init__( self.reactions = reactions or [] self.initial_conditions = initial_conditions or [] self.traps = traps or [] - self.temperature_fenics = None self._vtxfiles: list[dolfinx.io.VTXWriter] = [] @property @@ -177,45 +176,23 @@ def temperature(self, value): if value is None: self._temperature = value elif isinstance(value, (float, int, fem.Constant, fem.Function)): - self._temperature = value + self._temperature = festim.Value(value) elif callable(value): arguments = value.__code__.co_varnames if "T" in arguments: raise ValueError("Temperature cannot be a function of temperature, T") - self._temperature = value + self._temperature = festim.Value(value) else: raise TypeError( "Value must be a float, int, fem.Constant, fem.Function, or callable" ) - @property - def temperature_fenics(self): - return self._temperature_fenics - - @temperature_fenics.setter - def temperature_fenics(self, value): - if value is None: - self._temperature_fenics = value - return - elif not isinstance( - value, - festim.ConvertToFenicsObject, - ): - raise TypeError("Value must be a fem.Constant or fem.Function") - - self._temperature_fenics = value - @property def temperature_time_dependent(self): if self.temperature is None: return False - if isinstance(self.temperature, fem.Constant): - return False - if callable(self.temperature): - arguments = self.temperature.__code__.co_varnames - return "t" in arguments else: - return False + return self.temperature.time_dependent @property def multispecies(self): @@ -317,9 +294,10 @@ def define_temperature(self): if self.temperature is None: raise ValueError("the temperature attribute needs to be defined") - self.temperature_fenics = festim.ConvertToFenicsObject( - input_value=self.temperature - ) + if self.temperature.temperature_dependent: + raise ValueError( + "the temperature input value cannot be dependent on temperature" + ) degree = 1 element_temperature = basix.ufl.element( @@ -332,7 +310,7 @@ def define_temperature(self): self.mesh.mesh, element_temperature ) - self.temperature_fenics.convert_value( + self.temperature.convert_input_value( mesh=self.mesh.mesh, function_space=function_space_temperature, t=self.t ) @@ -420,7 +398,7 @@ def define_D_global(self, species): expr = D_0 * ufl.exp( -E_D / as_fenics_constant(k_B, self.mesh.mesh) - / self.temperature_fenics.fenics_object + / self.temperature.fenics_object ) D_expr = fem.Expression(expr, self.V_DG_1.element.interpolation_points()) D.interpolate(D_expr) @@ -515,18 +493,21 @@ def define_boundary_conditions(self): # 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): - if isinstance(bc.value, (int, float)): - bc.value_fenics.as_fenics_constant( + if isinstance( + bc.value.input_value, + (int, float, fem.Constant, fem.Function, ufl.core.expr.Expr), + ): + bc.value.convert_input_value( mesh=self.mesh.mesh, - value=bc.value, + t=self.t, + temperature=self.temperature.fenics_object, ) - elif isinstance(bc.value, (fem.Function, ufl.core.expr.Expr)): - bc.value_fenics.fenics_object = bc.value else: - bc.value_fenics.as_mapped_function( + bc.value.fenics_object = festim.as_mapped_function( + value=bc.value.input_value, mesh=self.mesh.mesh, - temperature=self.temperature_fenics.fenics_object, + temperature=self.temperature.fenics_object, t=self.t, ) @@ -548,9 +529,9 @@ def create_dirichletbc_form(self, bc): else: function_space_value = bc.species.collapsed_function_space - bc.value_fenics.convert_value( + bc.value.convert_input_value( mesh=self.mesh.mesh, - temperature=self.temperature_fenics.fenics_object, + temperature=self.temperature.fenics_object, function_space=function_space_value, t=self.t, ) @@ -572,16 +553,14 @@ def create_dirichletbc_form(self, bc): ) # create form - if not self.multispecies and isinstance( - bc.value_fenics.fenics_object, (fem.Function) - ): + if not self.multispecies and isinstance(bc.value.fenics_object, (fem.Function)): # no need to pass the functionspace since value_fenics is already a Function function_space_form = None else: function_space_form = bc.species.sub_function_space form = fem.dirichletbc( - value=bc.value_fenics.fenics_object, + value=bc.value.fenics_object, dofs=bc_dofs, V=function_space_form, ) @@ -593,18 +572,21 @@ def create_source_values_fenics(self): for source in self.sources: # create value_fenics for all F.ParticleSource objects if isinstance(source, festim.source.ParticleSource): - if isinstance(source.value, (int, float)): - source.value_fenics.as_fenics_constant( + if isinstance( + source.value.input_value, + (int, float, fem.Constant, fem.Function, ufl.core.expr.Expr), + ): + source.value.convert_input_value( mesh=self.mesh.mesh, - value=source.value, + t=self.t, + temperature=self.temperature.fenics_object, ) - elif isinstance(source.value, (fem.Function, ufl.core.expr.Expr)): - source.value_fenics.fenics_object = self.value - elif callable(source.value): - source.value_fenics.as_mapped_function( + elif callable(source.value.input_value): + source.value.fenics_object = festim.as_mapped_function( + value=source.value.input_value, mesh=self.mesh.mesh, - temperature=self.temperature_fenics.fenics_object, t=self.t, + temperature=self.temperature.fenics_object, ) def create_flux_values_fenics(self): @@ -613,9 +595,9 @@ def create_flux_values_fenics(self): # create value_fenics for all F.ParticleFluxBC objects if isinstance(bc, boundary_conditions.ParticleFluxBC): - bc.value_fenics.convert_value( + bc.value.convert_input_value( mesh=self.mesh.mesh, - temperature=self.temperature_fenics, + temperature=self.temperature.fenics_object, t=self.t, ) @@ -631,35 +613,43 @@ def create_initial_conditions(self): function_space_value = None for condition in self.initial_conditions: + # create value_fenics for condition function_space_value = None - if callable(condition.value): + if callable(condition.value.input_value): + + if condition.value.time_dependent: + raise ValueError("Initial conditions cannot be time dependent") + # if bc.value is a callable then need to provide a functionspace if not self.multispecies: function_space_value = condition.species.sub_function_space else: function_space_value = condition.species.collapsed_function_space - if isinstance(condition.value, (int, float)): - condition.value_fenics.fenics_interpolation_expression = ( - lambda x: np.full(x.shape[1], condition.value) + if isinstance(condition.value.input_value, (int, float)): + condition.value.fenics_interpolation_expression = lambda x: np.full( + x.shape[1], condition.value.input_value ) else: - condition.value_fenics.as_fenics_interpolation_expression( - mesh=self.mesh.mesh, - temperature=self.temperature_fenics.fenics_object, - function_space=function_space_value, + condition.value.fenics_interpolation_expression = ( + festim.as_fenics_interpolation_expression( + value=condition.value.input_value, + function_space=function_space_value, + mesh=self.mesh.mesh, + temperature=self.temperature.fenics_object, + ) ) # assign to previous solution of species if not self.multispecies: condition.species.prev_solution.interpolate( - condition.value_fenics.fenics_interpolation_expression + condition.value.fenics_interpolation_expression ) else: idx = self.species.index(condition.species) self.u_n.sub(idx).interpolate( - condition.value_fenics.fenics_interpolation_expression + condition.value.fenics_interpolation_expression ) def create_formulation(self): @@ -675,7 +665,7 @@ def create_formulation(self): for vol in self.volume_subdomains: D = vol.material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature_fenics.fenics_object, spe + self.mesh.mesh, self.temperature.fenics_object, spe ) if spe.mobile: self.formulation += ufl.dot(D * ufl.grad(u), ufl.grad(v)) * self.dx( @@ -689,7 +679,7 @@ def create_formulation(self): for reactant in reaction.reactant: if isinstance(reactant, festim.species.Species): self.formulation += ( - reaction.reaction_term(self.temperature_fenics.fenics_object) + reaction.reaction_term(self.temperature.fenics_object) * reactant.test_function * self.dx(reaction.volume.id) ) @@ -701,14 +691,14 @@ def create_formulation(self): products = [reaction.product] for product in products: self.formulation += ( - -reaction.reaction_term(self.temperature_fenics.fenics_object) + -reaction.reaction_term(self.temperature.fenics_object) * product.test_function * self.dx(reaction.volume.id) ) # add sources for source in self.sources: self.formulation -= ( - source.value_fenics.fenics_object + source.value.fenics_object * source.species.test_function * self.dx(source.volume.id) ) @@ -717,7 +707,7 @@ def create_formulation(self): for bc in self.boundary_conditions: if isinstance(bc, boundary_conditions.ParticleFluxBC): self.formulation -= ( - bc.value_fenics.fenics_object + bc.value.fenics_object * bc.species.test_function * self.ds(bc.subdomain.id) ) @@ -760,7 +750,7 @@ def update_time_dependent_values(self): if not self.temperature_time_dependent: return - self.temperature_fenics.update(t=t) + self.temperature.update(t=t) for bc in self.boundary_conditions: if isinstance( @@ -770,12 +760,12 @@ def update_time_dependent_values(self): boundary_conditions.ParticleFluxBC, ), ): - if bc.value_fenics.temperature_dependent: - bc.value_fenics.update(t=t) + if bc.value.temperature_dependent: + bc.value.update(t=t) for source in self.sources: - if source.value_fenics.temperature_dependent: - source.value_fenics.update(t=t) + if source.value.temperature_dependent: + source.value.update(t=t) def post_processing(self): """Post processes the model""" diff --git a/src/festim/initial_condition.py b/src/festim/initial_condition.py index ff9c621fb..6177709a1 100644 --- a/src/festim/initial_condition.py +++ b/src/festim/initial_condition.py @@ -31,11 +31,41 @@ def __init__(self, value, species): self.value = value self.species = species - self.value_fenics = F.ConvertToFenicsObject(value) + @property + def value(self): + return self._value + + @value.setter + def value(self, value): + if value is None: + self._value = value + elif isinstance(value, (float, int, fem.Constant, fem.Function)): + self._value = F.Value(value) + elif callable(value): + self._value = F.Value(value) + else: + raise TypeError( + "Value must be a float, int, fem.Constant, fem.Function, or callable" + ) class InitialTemperature: def __init__(self, value) -> None: self.value = value - self.value_fenics = F.ConvertToFenicsObject(value) + @property + def value(self): + return self._value + + @value.setter + def value(self, value): + if value is None: + self._value = value + elif isinstance(value, (float, int, fem.Constant, fem.Function)): + self._value = F.Value(value) + elif callable(value): + self._value = F.Value(value) + else: + raise TypeError( + "Value must be a float, int, fem.Constant, fem.Function, or callable" + ) diff --git a/src/festim/problem.py b/src/festim/problem.py index cfff53269..60dd512b9 100644 --- a/src/festim/problem.py +++ b/src/festim/problem.py @@ -189,9 +189,9 @@ def iterate(self): def update_time_dependent_values(self): t = float(self.t) for bc in self.boundary_conditions: - if bc.value_fenics.time_dependent: - bc.value_fenics.update(t=t) + if bc.value.time_dependent: + bc.value.update(t=t) for source in self.sources: - if source.value_fenics.time_dependent: - source.value_fenics.update(t=t) + if source.value.time_dependent: + source.value.update(t=t) diff --git a/src/festim/source.py b/src/festim/source.py index 153273bd9..8028cc388 100644 --- a/src/festim/source.py +++ b/src/festim/source.py @@ -1,6 +1,3 @@ -import dolfinx -import numpy as np -import ufl from dolfinx import fem import festim as F @@ -40,7 +37,22 @@ def __init__(self, value, volume): self.value = value self.volume = volume - self.value_fenics = F.ConvertToFenicsObject(value) + @property + def value(self): + return self._value + + @value.setter + def value(self, value): + if value is None: + self._value = value + elif isinstance(value, (float, int, fem.Constant, fem.Function)): + self._value = F.Value(value) + elif callable(value): + self._value = F.Value(value) + else: + raise TypeError( + "Value must be a float, int, fem.Constant, fem.Function, or callable" + ) @property def volume(self): @@ -75,3 +87,6 @@ def species(self, value): class HeatSource(SourceBase): def __init__(self, value, volume): super().__init__(value, volume) + + if self.value.temperature_dependent: + raise ValueError("Heat source cannot be temperature dependent") diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 1c0887643..6323114d8 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -107,7 +107,7 @@ def test_define_temperature(input, expected_type): my_model.define_temperature() # TEST - assert isinstance(my_model.temperature_fenics.fenics_object, expected_type) + assert isinstance(my_model.temperature.fenics_object, expected_type) @pytest.mark.parametrize( @@ -214,8 +214,8 @@ def test_update_time_dependent_values_temperature(T_function, expected_values): my_model.update_time_dependent_values() # TEST - if isinstance(my_model.temperature_fenics, fem.Constant): - computed_value = float(my_model.temperature_fenics) + if isinstance(my_model.temperature.fenics_object, fem.Constant): + computed_value = float(my_model.temperature.fenics_object) assert np.isclose(computed_value, expected_values[i]) @@ -301,8 +301,12 @@ def test_define_D_global_different_materials(): computed_values = [D_computed.x.array[0], D_computed.x.array[-1]] - D_expected_left = D_0_left * np.exp(-E_D_left / (F.k_B * my_model.temperature)) - D_expected_right = D_0_right * np.exp(-E_D_right / (F.k_B * my_model.temperature)) + D_expected_left = D_0_left * np.exp( + -E_D_left / (F.k_B * my_model.temperature.input_value) + ) + D_expected_right = D_0_right * np.exp( + -E_D_right / (F.k_B * my_model.temperature.input_value) + ) expected_values = [D_expected_left, D_expected_right] @@ -413,8 +417,8 @@ def test_define_D_global_multispecies(): computed_values = [D_A_computed.x.array[-1], D_B_computed.x.array[-1]] - D_expected_A = D_0_A * np.exp(-E_D_A / (F.k_B * my_model.temperature)) - D_expected_B = D_0_B * np.exp(-E_D_B / (F.k_B * my_model.temperature)) + D_expected_A = D_0_A * np.exp(-E_D_A / (F.k_B * my_model.temperature.input_value)) + D_expected_B = D_0_B * np.exp(-E_D_B / (F.k_B * my_model.temperature.input_value)) expected_values = [D_expected_A, D_expected_B] @@ -636,8 +640,10 @@ def test_update_time_dependent_bcs_with_time_dependent_temperature( my_model.update_time_dependent_values() # TEST - if isinstance(my_model.boundary_conditions[0].value_fenics, fem.Constant): - computed_value = float(my_model.boundary_conditions[0].value_fenics) + if isinstance( + my_model.boundary_conditions[0].value.fenics_object, fem.Constant + ): + computed_value = float(my_model.boundary_conditions[0].value.fenics_object) assert np.isclose(computed_value, expected_values[i]) @@ -681,8 +687,8 @@ def test_update_time_dependent_values_source(source_value, expected_values): my_model.update_time_dependent_values() # TEST - if isinstance(my_model.sources[0].value_fenics.fenics_object, fem.Constant): - computed_value = float(my_model.sources[0].value_fenics.fenics_object) + if isinstance(my_model.sources[0].value.fenics_object, fem.Constant): + computed_value = float(my_model.sources[0].value.fenics_object) assert np.isclose(computed_value, expected_values[i]) @@ -735,8 +741,8 @@ def test_update_sources_with_time_dependent_temperature( my_model.update_time_dependent_values() # TEST - if isinstance(my_model.sources[0].value_fenics.fenics_object, fem.Constant): - computed_value = float(my_model.sources[0].value_fenics.fenics_object) + if isinstance(my_model.sources[0].value.fenics_object, fem.Constant): + computed_value = float(my_model.sources[0].value.fenics_object) assert np.isclose(computed_value, expected_values[i]) @@ -767,8 +773,8 @@ def test_create_source_values_fenics_multispecies(): my_model.create_source_values_fenics() # TEST - assert np.isclose(float(my_model.sources[0].value_fenics.fenics_object), 5) - assert np.isclose(float(my_model.sources[1].value_fenics.fenics_object), 11) + assert np.isclose(float(my_model.sources[0].value.fenics_object), 5) + assert np.isclose(float(my_model.sources[1].value.fenics_object), 11) # TODO replace this by a proper MMS test @@ -1146,8 +1152,10 @@ def test_update_time_dependent_values_flux(bc_value, expected_values): my_model.update_time_dependent_values() # TEST - if isinstance(my_model.boundary_conditions[0].value_fenics, fem.Constant): - computed_value = float(my_model.boundary_conditions[0].value_fenics) + if isinstance( + my_model.boundary_conditions[0].value.fenics_object, fem.Constant + ): + computed_value = float(my_model.boundary_conditions[0].value.fenics_object) assert np.isclose(computed_value, expected_values[i]) @@ -1200,8 +1208,10 @@ def test_update_fluxes_with_time_dependent_temperature( my_model.update_time_dependent_values() # TEST - if isinstance(my_model.boundary_conditions[0].value_fenics, fem.Constant): - computed_value = float(my_model.boundary_conditions[0].value_fenics) + if isinstance( + my_model.boundary_conditions[0].value.fenics_object, fem.Constant + ): + computed_value = float(my_model.boundary_conditions[0].value.fenics_object) assert np.isclose(computed_value, expected_values[i]) @@ -1233,9 +1243,5 @@ def test_create_flux_values_fenics_multispecies(): my_model.create_flux_values_fenics() # TEST - assert np.isclose( - float(my_model.boundary_conditions[0].value_fenics.fenics_object), 5 - ) - assert np.isclose( - float(my_model.boundary_conditions[1].value_fenics.fenics_object), 11 - ) + assert np.isclose(float(my_model.boundary_conditions[0].value.fenics_object), 5) + assert np.isclose(float(my_model.boundary_conditions[1].value.fenics_object), 11) diff --git a/test/test_heat_transfer_problem.py b/test/test_heat_transfer_problem.py index 224802b7e..336da634c 100644 --- a/test/test_heat_transfer_problem.py +++ b/test/test_heat_transfer_problem.py @@ -603,7 +603,6 @@ def test_adaptive_timestepping_shrinks(): (lambda t: t, [1.0, 2.0, 3.0]), (lambda t: 1.0 + t, [2.0, 3.0, 4.0]), (lambda x, t: 1.0 + x[0] + t, [6.0, 7.0, 8.0]), - (lambda T, t: T + 2 * t, [12.0, 14.0, 16.0]), ( lambda x, t: ufl.conditional(ufl.lt(t, 1.5), 100.0 + x[0], 0.0), [104.0, 0.0, 0.0], @@ -639,6 +638,8 @@ def test_update_time_dependent_values_HeatFluxBC(bc_value, expected_values): my_model.update_time_dependent_values() # TEST - if isinstance(my_model.boundary_conditions[0].value_fenics, fem.Constant): - computed_value = float(my_model.boundary_conditions[0].value_fenics) + if isinstance( + my_model.boundary_conditions[0].value.fenics_object, fem.Constant + ): + computed_value = float(my_model.boundary_conditions[0].value.fenics_object) assert np.isclose(computed_value, expected_values[i]) From 76ab530cc5eaac475d366d786f8b4c1f803fc1bb Mon Sep 17 00:00:00 2001 From: jhdark Date: Mon, 27 Jan 2025 11:22:34 -0500 Subject: [PATCH 03/40] addtional argument for convert value --- src/festim/helpers.py | 39 +++++++------------- src/festim/hydrogen_transport_problem.py | 47 +++++++----------------- 2 files changed, 26 insertions(+), 60 deletions(-) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 122edb7f7..07214ca30 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -5,31 +5,6 @@ import ufl -# def as_fenics_constant( -# value: float | int | fem.Constant, mesh: dolfinx.mesh.Mesh -# ) -> fem.Constant: -# """Converts a value to a dolfinx.Constant - -# Args: -# value: the value to convert -# mesh: the mesh of the domiain - -# Returns: -# The converted value - -# Raises: -# TypeError: if the value is not a float, an int or a dolfinx.Constant -# """ -# if isinstance(value, (float, int)): -# return fem.Constant(mesh, dolfinx.default_scalar_type(value)) -# elif isinstance(value, fem.Constant): -# return value -# else: -# raise TypeError( -# f"Value must be a float, an int or a dolfinx.Constant, not {type(value)}" -# ) - - def as_fenics_constant( value: float | int | fem.Constant, mesh: dolfinx.mesh.Mesh ) -> fem.Constant: @@ -211,7 +186,12 @@ def temperature_dependent(self) -> bool: return False def convert_input_value( - self, mesh=None, function_space=None, t=None, temperature=None + self, + mesh=None, + function_space=None, + t=None, + temperature=None, + up_to_mapping=False, ): """Converts a user given value to a relevent fenics object depending on the type of the value provided @@ -221,6 +201,8 @@ def convert_input_value( function_space (dolfinx.fem.function.FunctionSpace): the function space of the fenics object t (fem.Constant): the time temperature (fem.Function, fem.Constant or ufl.core.expr.Expr): the temperature + up_to_mapping (bool): if True, the value is only mapped to a function if the input is callable, + not interpolated or converted to a function """ if isinstance(self.input_value, fem.Constant): self.fenics_object = self.input_value @@ -248,6 +230,11 @@ def convert_input_value( value=self.input_value(t=float(t)), mesh=mesh ) + elif up_to_mapping: + self.fenics_object = as_mapped_function( + value=self.input_value, mesh=mesh, t=t, temperature=temperature + ) + elif self.temperature_dependent: self.fenics_interpolation_expression, self.fenics_object = ( as_fenics_interp_expr_and_function( diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index b50575284..46f995de5 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -493,23 +493,13 @@ def define_boundary_conditions(self): # 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): - if isinstance( - bc.value.input_value, - (int, float, fem.Constant, fem.Function, ufl.core.expr.Expr), - ): - bc.value.convert_input_value( - mesh=self.mesh.mesh, - t=self.t, - temperature=self.temperature.fenics_object, - ) - else: - bc.value.fenics_object = festim.as_mapped_function( - value=bc.value.input_value, - mesh=self.mesh.mesh, - temperature=self.temperature.fenics_object, - t=self.t, - ) + bc.value.convert_input_value( + mesh=self.mesh.mesh, + t=self.t, + temperature=self.temperature.fenics_object, + up_to_mapping=True, + ) super().define_boundary_conditions() @@ -571,34 +561,23 @@ def create_source_values_fenics(self): """For each source create the value_fenics""" for source in self.sources: # create value_fenics for all F.ParticleSource objects - if isinstance(source, festim.source.ParticleSource): - if isinstance( - source.value.input_value, - (int, float, fem.Constant, fem.Function, ufl.core.expr.Expr), - ): - source.value.convert_input_value( - mesh=self.mesh.mesh, - t=self.t, - temperature=self.temperature.fenics_object, - ) - elif callable(source.value.input_value): - source.value.fenics_object = festim.as_mapped_function( - value=source.value.input_value, - mesh=self.mesh.mesh, - t=self.t, - temperature=self.temperature.fenics_object, - ) + source.value.convert_input_value( + mesh=self.mesh.mesh, + t=self.t, + temperature=self.temperature.fenics_object, + up_to_mapping=True, + ) 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, boundary_conditions.ParticleFluxBC): - bc.value.convert_input_value( mesh=self.mesh.mesh, temperature=self.temperature.fenics_object, t=self.t, + up_to_mapping=True, ) def create_initial_conditions(self): From 19d6b28483c4bbda235a75f0ae2ee65aa92f98f2 Mon Sep 17 00:00:00 2001 From: jhdark Date: Mon, 27 Jan 2025 14:45:28 -0500 Subject: [PATCH 04/40] add repr dunder for when users print the value --- src/festim/helpers.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 07214ca30..69aa24fe5 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -131,6 +131,9 @@ def __init__(self, input_value): self.fenics_interpolation_expression = None self.fenics_object = None + def __repr__(self) -> str: + return str(self.input_value) + @property def input_value(self): return self._input_value From 0bf64c4d3da7d17823752d86a7dbccc1bc3777ef Mon Sep 17 00:00:00 2001 From: jhdark Date: Tue, 28 Jan 2025 16:25:32 -0500 Subject: [PATCH 05/40] better naming --- src/festim/__init__.py | 2 +- src/festim/heat_transfer_problem.py | 25 +++++------------------- src/festim/helpers.py | 18 ++++++++--------- src/festim/hydrogen_transport_problem.py | 6 +++--- 4 files changed, 17 insertions(+), 34 deletions(-) diff --git a/src/festim/__init__.py b/src/festim/__init__.py index d0f4c0f58..2a95882f1 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -38,7 +38,7 @@ from .heat_transfer_problem import HeatTransferProblem from .helpers import ( as_fenics_constant, - as_mapped_function, + as_ufl_expression, as_fenics_interpolation_expression, as_fenics_interp_expr_and_function, Value, diff --git a/src/festim/heat_transfer_problem.py b/src/festim/heat_transfer_problem.py index 90b1207c0..4c69db779 100644 --- a/src/festim/heat_transfer_problem.py +++ b/src/festim/heat_transfer_problem.py @@ -137,20 +137,11 @@ def create_source_values_fenics(self): """For each source create the value_fenics""" for source in self.sources: # create value_fenics for all source objects - if isinstance( - source.value.input_value, - (int, float, fem.Constant, fem.Function, ufl.core.expr.Expr), - ): - source.value.convert_input_value( - mesh=self.mesh.mesh, - t=self.t, - ) - elif callable(source.value.input_value): - source.value.fenics_object = helpers.as_mapped_function( - value=source.value.input_value, - mesh=self.mesh.mesh, - t=self.t, - ) + source.value.convert_input_value( + mesh=self.mesh.mesh, + t=self.t, + up_to_ufl_expr=True, + ) def create_flux_values_fenics(self): """For each heat flux create the value_fenics""" @@ -176,12 +167,6 @@ def create_initial_conditions(self): "Initial conditions can only be defined for transient simulations" ) - # create value_fenics for condition - - # self.initial_condition.value.create_expr_fenics( - # mesh=self.mesh.mesh, - # function_space=self.function_space, - # ) if isinstance(self.initial_condition.value.input_value, (int, float)): self.initial_condition.value.fenics_interpolation_expression = ( lambda x: np.full(x.shape[1], self.initial_condition.value.input_value) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 69aa24fe5..998306330 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -30,9 +30,7 @@ def as_fenics_constant( ) -def as_mapped_function( - value, mesh=None, t=None, temperature=None -) -> ufl.core.expr.Expr: +def as_ufl_expression(value, mesh=None, t=None, temperature=None) -> ufl.core.expr.Expr: """Maps a user given function Args: @@ -64,7 +62,7 @@ def as_fenics_interpolation_expression( ) -> fem.Expression: """Converts a callable input value to a fenics expression""" - mapped_function = as_mapped_function( + mapped_function = as_ufl_expression( value=value, mesh=mesh, t=t, temperature=temperature ) @@ -119,7 +117,7 @@ class Value: | fem.Function ) - mapped_function: ufl.core.expr.Expr + ufl_expression: ufl.core.expr.Expr fenics_interpolation_expression: fem.Expression fenics_object: fem.Function | fem.Constant | ufl.core.expr.Expr @@ -127,7 +125,7 @@ def __init__(self, input_value): self, self.input_value = input_value - self.mapped_function = None + self.ufl_expression = None self.fenics_interpolation_expression = None self.fenics_object = None @@ -194,7 +192,7 @@ def convert_input_value( function_space=None, t=None, temperature=None, - up_to_mapping=False, + up_to_ufl_expr=False, ): """Converts a user given value to a relevent fenics object depending on the type of the value provided @@ -204,7 +202,7 @@ def convert_input_value( function_space (dolfinx.fem.function.FunctionSpace): the function space of the fenics object t (fem.Constant): the time temperature (fem.Function, fem.Constant or ufl.core.expr.Expr): the temperature - up_to_mapping (bool): if True, the value is only mapped to a function if the input is callable, + up_to_ufl_expr (bool): if True, the value is only mapped to a function if the input is callable, not interpolated or converted to a function """ if isinstance(self.input_value, fem.Constant): @@ -233,8 +231,8 @@ def convert_input_value( value=self.input_value(t=float(t)), mesh=mesh ) - elif up_to_mapping: - self.fenics_object = as_mapped_function( + elif up_to_ufl_expr: + self.fenics_object = as_ufl_expression( value=self.input_value, mesh=mesh, t=t, temperature=temperature ) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 46f995de5..ead48bb24 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -498,7 +498,7 @@ def define_boundary_conditions(self): mesh=self.mesh.mesh, t=self.t, temperature=self.temperature.fenics_object, - up_to_mapping=True, + up_to_ufl_expr=True, ) super().define_boundary_conditions() @@ -565,7 +565,7 @@ def create_source_values_fenics(self): mesh=self.mesh.mesh, t=self.t, temperature=self.temperature.fenics_object, - up_to_mapping=True, + up_to_ufl_expr=True, ) def create_flux_values_fenics(self): @@ -577,7 +577,7 @@ def create_flux_values_fenics(self): mesh=self.mesh.mesh, temperature=self.temperature.fenics_object, t=self.t, - up_to_mapping=True, + up_to_ufl_expr=True, ) def create_initial_conditions(self): From 7ef47913339ee1539e442495c1daa68e8bc1404c Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 31 Jan 2025 16:43:45 -0500 Subject: [PATCH 06/40] back to as mapped function, typehinting --- src/festim/__init__.py | 3 +- src/festim/heat_transfer_problem.py | 4 +- src/festim/helpers.py | 86 ++++++++++++------------ src/festim/hydrogen_transport_problem.py | 4 +- 4 files changed, 48 insertions(+), 49 deletions(-) diff --git a/src/festim/__init__.py b/src/festim/__init__.py index 2a95882f1..78ae73722 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -38,8 +38,7 @@ from .heat_transfer_problem import HeatTransferProblem from .helpers import ( as_fenics_constant, - as_ufl_expression, - as_fenics_interpolation_expression, + as_mapped_function, as_fenics_interp_expr_and_function, Value, ) diff --git a/src/festim/heat_transfer_problem.py b/src/festim/heat_transfer_problem.py index 4c69db779..268bfca90 100644 --- a/src/festim/heat_transfer_problem.py +++ b/src/festim/heat_transfer_problem.py @@ -172,8 +172,8 @@ def create_initial_conditions(self): lambda x: np.full(x.shape[1], self.initial_condition.value.input_value) ) else: - self.initial_condition.value.fenics_interpolation_expression = ( - helpers.as_fenics_interpolation_expression( + self.initial_condition.value.fenics_interpolation_expression, _ = ( + helpers.as_fenics_interp_expr_and_function( value=self.initial_condition.value.input_value, function_space=self.function_space, mesh=self.mesh.mesh, diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 998306330..0c39df50d 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -15,7 +15,7 @@ def as_fenics_constant( mesh: the mesh of the domiain Returns: - (fem.Constant) The converted value + The converted value Raises: TypeError: if the value is not a float, an int or a dolfinx.Constant @@ -30,17 +30,22 @@ def as_fenics_constant( ) -def as_ufl_expression(value, mesh=None, t=None, temperature=None) -> ufl.core.expr.Expr: - """Maps a user given function +def as_mapped_function( + value: Callable, + mesh: dolfinx.mesh.Mesh = None, + t: fem.Constant = None, + temperature: fem.Function | fem.Constant | ufl.core.expr.Expr = None, +) -> ufl.core.expr.Expr: + """Maps a user given callable function to the mesh, time or temperature within festim as needed Args: - value (Callable): the callable to convert - mesh (dolfinx.mesh.Mesh): the mesh of the domain - t (fem.Constant): the time - temperature (fem.Function, fem.Constant or ufl.core.expr.Expr): the temperature + value: the callable to convert + mesh: the mesh of the domain + t: the time + temperature: the temperature Returns: - (ufl.core.expr.Expr) The mapped function + The mapped function """ arguments = value.__code__.co_varnames @@ -57,32 +62,37 @@ def as_ufl_expression(value, mesh=None, t=None, temperature=None) -> ufl.core.ex return value(**kwargs) -def as_fenics_interpolation_expression( - value, function_space, mesh=None, t=None, temperature=None -) -> fem.Expression: - """Converts a callable input value to a fenics expression""" +def as_fenics_interp_expr_and_function( + value: Callable, + function_space: dolfinx.fem.function.FunctionSpace, + mesh: dolfinx.mesh.Mesh = None, + t: fem.Constant = None, + temperature: fem.Function | fem.Constant | ufl.core.expr.Expr = None, +) -> tuple[fem.Expression, fem.Function]: + """Takes a user given callable function, maps the function to the mesh, time or + temperature within festim as needed. Then creates the fenics interpolation expression + and function objects + + Args: + value: the callable to convert + function_space: The function space to interpolate function over + mesh: the mesh of the domain + t: the time + temperature: the temperature + + Returns: + fenics interpolation expression, fenics function + """ - mapped_function = as_ufl_expression( + mapped_function = as_mapped_function( value=value, mesh=mesh, t=t, temperature=temperature ) - return fem.Expression( + fenics_interpolation_expression = fem.Expression( mapped_function, function_space.element.interpolation_points(), ) - -def as_fenics_interp_expr_and_function( - value, function_space, mesh=None, t=None, temperature=None -) -> fem.Function: - - fenics_interpolation_expression = as_fenics_interpolation_expression( - value=value, - function_space=function_space, - mesh=mesh, - t=t, - temperature=temperature, - ) fenics_object = fem.Function(function_space) fenics_object.interpolate(fenics_interpolation_expression) @@ -188,11 +198,11 @@ def temperature_dependent(self) -> bool: def convert_input_value( self, - mesh=None, - function_space=None, - t=None, - temperature=None, - up_to_ufl_expr=False, + function_space: dolfinx.fem.function.FunctionSpace = None, + mesh: dolfinx.mesh.Mesh = None, + t: fem.Constant = None, + temperature: fem.Function | fem.Constant | ufl.core.expr.Expr = None, + up_to_ufl_expr: bool = False, ): """Converts a user given value to a relevent fenics object depending on the type of the value provided @@ -232,21 +242,10 @@ def convert_input_value( ) elif up_to_ufl_expr: - self.fenics_object = as_ufl_expression( + self.fenics_object = as_mapped_function( value=self.input_value, mesh=mesh, t=t, temperature=temperature ) - elif self.temperature_dependent: - self.fenics_interpolation_expression, self.fenics_object = ( - as_fenics_interp_expr_and_function( - value=self.input_value, - function_space=function_space, - mesh=mesh, - t=t, - temperature=temperature, - ) - ) - else: self.fenics_interpolation_expression, self.fenics_object = ( as_fenics_interp_expr_and_function( @@ -254,6 +253,7 @@ def convert_input_value( function_space=function_space, mesh=mesh, t=t, + temperature=temperature, ) ) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index ead48bb24..879796e5a 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -611,8 +611,8 @@ def create_initial_conditions(self): x.shape[1], condition.value.input_value ) else: - condition.value.fenics_interpolation_expression = ( - festim.as_fenics_interpolation_expression( + condition.value.fenics_interpolation_expression, _ = ( + festim.as_fenics_interp_expr_and_function( value=condition.value.input_value, function_space=function_space_value, mesh=self.mesh.mesh, From 808286e4507213b88cb51ae2cc492458f5c3bdb3 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 10:46:26 -0500 Subject: [PATCH 07/40] revert dirichet bc file --- .../boundary_conditions/dirichlet_bc.py | 194 ++++++++++++++---- 1 file changed, 150 insertions(+), 44 deletions(-) diff --git a/src/festim/boundary_conditions/dirichlet_bc.py b/src/festim/boundary_conditions/dirichlet_bc.py index 80e834363..c9cfc37c0 100644 --- a/src/festim/boundary_conditions/dirichlet_bc.py +++ b/src/festim/boundary_conditions/dirichlet_bc.py @@ -49,22 +49,38 @@ def __init__( self.subdomain = subdomain self.value = value + self.value_fenics = None + self.bc_expr = None + @property - def value(self): - return self._value + def value_fenics(self): + return self._value_fenics - @value.setter - def value(self, value): + @value_fenics.setter + def value_fenics(self, value: None | fem.Function | fem.Constant | np.ndarray): if value is None: - self._value = value - elif isinstance(value, (float, int, fem.Constant, fem.Function)): - self._value = helpers.Value(value) - elif callable(value): - self._value = helpers.Value(value) - else: + self._value_fenics = value + return + if not isinstance(value, (fem.Function, fem.Constant, np.ndarray)): + # FIXME: Should we allow sending in a callable here? raise TypeError( - "Value must be a float, int, fem.Constant, fem.Function, or callable" + "Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, or a np.ndarray not" + + f"{type(value)}" ) + self._value_fenics = value + + @property + def time_dependent(self) -> bool: + """Returns true if the value of the boundary condition is time dependent""" + if self.value is None: + return False + if isinstance(self.value, fem.Constant): + return False + if callable(self.value): + arguments = self.value.__code__.co_varnames + return "t" in arguments + else: + return False def define_surface_subdomain_dofs( self, @@ -104,6 +120,18 @@ def define_surface_subdomain_dofs( return bc_dofs + def update(self, t: float): + """Updates the boundary condition value + Args: + t (float): the time + """ + if callable(self.value): + arguments = self.value.__code__.co_varnames + if isinstance(self.value_fenics, fem.Constant) and "t" in arguments: + self.value_fenics.value = self.value(t=t) + else: + self.value_fenics.interpolate(self.bc_expr) + class FixedConcentrationBC(DirichletBCBase): """ @@ -145,42 +173,120 @@ def __init__( self.species = species super().__init__(subdomain, value) + @property + def temperature_dependent(self): + if self.value is None: + return False + if isinstance(self.value, fem.Constant): + return False + if callable(self.value): + arguments = self.value.__code__.co_varnames + return "T" in arguments + else: + return False + + def create_value( + self, + function_space: fem.FunctionSpace, + temperature: float | fem.Constant, + t: float | 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 `dolfinx.fem.Constant`. + If the value is a function of t, it is converted to `dolfinx.fem.Constant`. + Otherwise, it is converted to a `dolfinx.fem.Function`.Function and the + expression of the function is stored in `bc_expr`. + Args: + function_space (dolfinx.fem.FunctionSpace): the function space + temperature: The temperature + t (dolfinx.fem.Constant): the time + """ + mesh = function_space.mesh + x = ufl.SpatialCoordinate(mesh) + + if isinstance(self.value, (int, float)): + self.value_fenics = helpers.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 + if not isinstance(self.value(t=float(t)), (float, int)): + raise ValueError( + "self.value should return a float or an int, not " + + f"{type(self.value(t=float(t)))} " + ) + self.value_fenics = helpers.as_fenics_constant( + mesh=mesh, value=self.value(t=float(t)) + ) + else: + self.value_fenics = fem.Function(function_space) + kwargs = {} + if "t" in arguments: + kwargs["t"] = t + if "x" in arguments: + 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(), + ) + self.value_fenics.interpolate(self.bc_expr) + # alias for FixedConcentrationBC DirichletBC = FixedConcentrationBC class FixedTemperatureBC(DirichletBCBase): - """ - 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 - - Examples: - - .. highlight:: python - .. code-block:: python - - from festim import FixedTemperatureBC - FixedTemperatureBC(subdomain=my_subdomain, value=1) - FixedTemperatureBC(subdomain=my_subdomain, - value=lambda x: 1 + x[0]) - FixedTemperatureBC(subdomain=my_subdomain, - value=lambda t: 1 + t) - FixedTemperatureBC(subdomain=my_subdomain, - value=lambda x, t: 1 + x[0] + t) - - """ - - def __init__( - self, - subdomain: _subdomain.SurfaceSubdomain, - value: np.ndarray | fem.Constant | int | float | Callable, - ): - super().__init__(subdomain, value) - - if self.value.temperature_dependent: - raise ValueError( - "Temperature dependent boundary conditions are not supported for FixedTemperatureBC" - ) + def create_value(self, function_space: fem.FunctionSpace, 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 `dolfinx.fem.Constant`. + If the value is a function of t, it is converted to a `dolfinx.fem.Constant`. + Otherwise, it is converted to a` dolfinx.fem.Function` and the + expression of the function is stored in `bc_expr`. + Args: + function_space: the function space + t: the time + """ + mesh = function_space.mesh + x = ufl.SpatialCoordinate(mesh) + + if isinstance(self.value, (int, float)): + self.value_fenics = helpers.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: + # only t is an argument + if not isinstance(self.value(t=float(t)), (float, int)): + raise ValueError( + "self.value should return a float or an int, not " + + f"{type(self.value(t=float(t)))} " + ) + self.value_fenics = helpers.as_fenics_constant( + mesh=mesh, value=self.value(t=float(t)) + ) + else: + self.value_fenics = fem.Function(function_space) + kwargs = {} + if "t" in arguments: + kwargs["t"] = t + if "x" in arguments: + kwargs["x"] = x + + # 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(), + ) + self.value_fenics.interpolate(self.bc_expr) From a8ef53652bf6c6eef82fda62d95d342bc3875dfc Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 10:47:28 -0500 Subject: [PATCH 08/40] add spaces --- src/festim/boundary_conditions/dirichlet_bc.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/festim/boundary_conditions/dirichlet_bc.py b/src/festim/boundary_conditions/dirichlet_bc.py index c9cfc37c0..0a5af01d6 100644 --- a/src/festim/boundary_conditions/dirichlet_bc.py +++ b/src/festim/boundary_conditions/dirichlet_bc.py @@ -122,6 +122,7 @@ def define_surface_subdomain_dofs( def update(self, t: float): """Updates the boundary condition value + Args: t (float): the time """ @@ -197,6 +198,7 @@ def create_value( If the value is a function of t, it is converted to `dolfinx.fem.Constant`. Otherwise, it is converted to a `dolfinx.fem.Function`.Function and the expression of the function is stored in `bc_expr`. + Args: function_space (dolfinx.fem.FunctionSpace): the function space temperature: The temperature @@ -252,6 +254,7 @@ def create_value(self, function_space: fem.FunctionSpace, t: fem.Constant): If the value is a function of t, it is converted to a `dolfinx.fem.Constant`. Otherwise, it is converted to a` dolfinx.fem.Function` and the expression of the function is stored in `bc_expr`. + Args: function_space: the function space t: the time From 61e2ab404bdc3244edc0ec9aef051bc0204b5393 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 10:49:32 -0500 Subject: [PATCH 09/40] revert flux bc --- src/festim/boundary_conditions/flux_bc.py | 145 +++++++++++++++++++--- 1 file changed, 131 insertions(+), 14 deletions(-) diff --git a/src/festim/boundary_conditions/flux_bc.py b/src/festim/boundary_conditions/flux_bc.py index 3a72c1fc0..b908b353e 100644 --- a/src/festim/boundary_conditions/flux_bc.py +++ b/src/festim/boundary_conditions/flux_bc.py @@ -1,3 +1,5 @@ +import numpy as np +import ufl from dolfinx import fem import festim as F @@ -32,22 +34,98 @@ def __init__(self, subdomain, value): self.subdomain = subdomain self.value = value + self.value_fenics = None + self.bc_expr = None + @property - def value(self): - return self._value + def value_fenics(self): + return self._value_fenics - @value.setter - def value(self, value): + @value_fenics.setter + def value_fenics(self, value): if value is None: - self._value = value - elif isinstance(value, (float, int, fem.Constant, fem.Function)): - self._value = F.Value(value) - elif callable(value): - self._value = F.Value(value) - else: + self._value_fenics = value + return + if not isinstance( + value, (fem.Function, fem.Constant, np.ndarray, ufl.core.expr.Expr) + ): raise TypeError( - "Value must be a float, int, fem.Constant, fem.Function, or callable" + f"Value must be a dolfinx.fem.Function, dolfinx.fem.Constant, np.ndarray or ufl.core.expr.Expr not {type(value)}" ) + self._value_fenics = value + + @property + def time_dependent(self): + if self.value is None: + raise TypeError("Value must be given to determine if its time dependent") + if isinstance(self.value, fem.Constant): + return False + if callable(self.value): + arguments = self.value.__code__.co_varnames + return "t" in arguments + else: + return False + + @property + def temperature_dependent(self): + if self.value is None: + return False + if isinstance(self.value, fem.Constant): + return False + if callable(self.value): + arguments = self.value.__code__.co_varnames + return "T" in arguments + else: + return False + + 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 + 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 ufl Expression + Args: + mesh (dolfinx.mesh.Mesh) : the mesh + 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 + if not isinstance(self.value(t=float(t)), (float, int)): + raise ValueError( + f"self.value should return a float or an int, not {type(self.value(t=float(t)))} " + ) + self.value_fenics = F.as_fenics_constant( + mesh=mesh, value=self.value(t=float(t)) + ) + else: + kwargs = {} + if "t" in arguments: + kwargs["t"] = t + if "x" in arguments: + kwargs["x"] = x + if "T" in arguments: + kwargs["T"] = temperature + + self.value_fenics = self.value(**kwargs) + + def update(self, t): + """Updates the flux bc value + Args: + t (float): the time + """ + if callable(self.value): + arguments = self.value.__code__.co_varnames + if isinstance(self.value_fenics, fem.Constant) and "t" in arguments: + self.value_fenics.value = self.value(t=t) class ParticleFluxBC(FluxBCBase): @@ -96,6 +174,48 @@ def __init__(self, subdomain, value, species, species_dependent_value={}): self.species = species self.species_dependent_value = species_dependent_value + 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 + 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 ufl Expression + Args: + mesh (dolfinx.mesh.Mesh) : the mesh + 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 + if not isinstance(self.value(t=float(t)), (float, int)): + raise ValueError( + f"self.value should return a float or an int, not {type(self.value(t=float(t)))} " + ) + self.value_fenics = F.as_fenics_constant( + mesh=mesh, value=self.value(t=float(t)) + ) + else: + kwargs = {} + if "t" in arguments: + kwargs["t"] = t + if "x" in arguments: + kwargs["x"] = x + if "T" in arguments: + kwargs["T"] = temperature + + for name, species in self.species_dependent_value.items(): + kwargs[name] = species.concentration + + self.value_fenics = self.value(**kwargs) + class HeatFluxBC(FluxBCBase): """ @@ -131,6 +251,3 @@ class HeatFluxBC(FluxBCBase): def __init__(self, subdomain, value): super().__init__(subdomain=subdomain, value=value) - - if self.value.temperature_dependent: - raise ValueError("Heat flux cannot be temperature dependent") From 595ec793b214a1c7346b75b524d3af578571c124 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 10:50:19 -0500 Subject: [PATCH 10/40] add spaces --- src/festim/boundary_conditions/flux_bc.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/festim/boundary_conditions/flux_bc.py b/src/festim/boundary_conditions/flux_bc.py index b908b353e..99349a0f1 100644 --- a/src/festim/boundary_conditions/flux_bc.py +++ b/src/festim/boundary_conditions/flux_bc.py @@ -84,6 +84,7 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant): 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 ufl Expression + Args: mesh (dolfinx.mesh.Mesh) : the mesh temperature (float): the temperature @@ -119,6 +120,7 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant): def update(self, t): """Updates the flux bc value + Args: t (float): the time """ @@ -180,6 +182,7 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant): 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 ufl Expression + Args: mesh (dolfinx.mesh.Mesh) : the mesh temperature (float): the temperature From fa2741157fe33cd7aa9c3fc7343e3c0c714147cb Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 10:52:11 -0500 Subject: [PATCH 11/40] revert initial condition --- src/festim/initial_condition.py | 88 ++++++++++++++++++++++----------- 1 file changed, 58 insertions(+), 30 deletions(-) diff --git a/src/festim/initial_condition.py b/src/festim/initial_condition.py index 6177709a1..0251bac21 100644 --- a/src/festim/initial_condition.py +++ b/src/festim/initial_condition.py @@ -31,41 +31,69 @@ def __init__(self, value, species): self.value = value self.species = species - @property - def value(self): - return self._value - - @value.setter - def value(self, value): - if value is None: - self._value = value - elif isinstance(value, (float, int, fem.Constant, fem.Function)): - self._value = F.Value(value) - elif callable(value): - self._value = F.Value(value) - else: - raise TypeError( - "Value must be a float, int, fem.Constant, fem.Function, or callable" + self.expr_fenics = None + + def create_expr_fenics(self, mesh, temperature, function_space): + """Creates the expr_fenics of the initial condition. + If the value is a float or int, a function is created with an array with + the shape of the mesh and all set to the value. + Otherwise, it is converted to a fem.Expression. + + Args: + mesh (dolfinx.mesh.Mesh) : the mesh + temperature (float): the temperature + function_space(dolfinx.fem.FunctionSpaceBase): the function space of the species + """ + x = ufl.SpatialCoordinate(mesh) + + if isinstance(self.value, (int, float)): + self.expr_fenics = lambda x: np.full(x.shape[1], self.value) + + elif callable(self.value): + arguments = self.value.__code__.co_varnames + kwargs = {} + if "t" in arguments: + raise ValueError("Initial condition cannot be a function of time.") + if "x" in arguments: + kwargs["x"] = x + if "T" in arguments: + kwargs["T"] = temperature + + self.expr_fenics = fem.Expression( + self.value(**kwargs), + function_space.element.interpolation_points(), ) class InitialTemperature: def __init__(self, value) -> None: self.value = value + self.expr_fenics = None + + def create_expr_fenics(self, mesh, function_space): + """Creates the expr_fenics of the initial condition. + If the value is a float or int, a function is created with an array with + the shape of the mesh and all set to the value. + Otherwise, it is converted to a fem.Expression. + + Args: + mesh (dolfinx.mesh.Mesh) : the mesh + function_space(dolfinx.fem.FunctionSpace): the function space of the species + """ + x = ufl.SpatialCoordinate(mesh) + + if isinstance(self.value, (int, float)): + self.expr_fenics = lambda x: np.full(x.shape[1], self.value) + + elif callable(self.value): + arguments = self.value.__code__.co_varnames + kwargs = {} + if "t" in arguments: + raise ValueError("Initial condition cannot be a function of time.") + if "x" in arguments: + kwargs["x"] = x - @property - def value(self): - return self._value - - @value.setter - def value(self, value): - if value is None: - self._value = value - elif isinstance(value, (float, int, fem.Constant, fem.Function)): - self._value = F.Value(value) - elif callable(value): - self._value = F.Value(value) - else: - raise TypeError( - "Value must be a float, int, fem.Constant, fem.Function, or callable" + self.expr_fenics = fem.Expression( + self.value(**kwargs), + function_space.element.interpolation_points(), ) From d57e3b0256d5ae2f5a70a84a6a098661665c7e89 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 10:53:50 -0500 Subject: [PATCH 12/40] remove import festim --- src/festim/initial_condition.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/festim/initial_condition.py b/src/festim/initial_condition.py index 0251bac21..b4947c521 100644 --- a/src/festim/initial_condition.py +++ b/src/festim/initial_condition.py @@ -1,7 +1,6 @@ import numpy as np import ufl from dolfinx import fem -import festim as F # TODO rename this to InitialConcentration and create a new base class From 4f52b044a3b8ffd56d08abf174f8417b7b6549e4 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 11:01:20 -0500 Subject: [PATCH 13/40] revert heat transfer problem --- src/festim/heat_transfer_problem.py | 41 ++++++++++------------------- src/festim/problem.py | 4 +-- test/test_heat_transfer_problem.py | 6 ++--- 3 files changed, 18 insertions(+), 33 deletions(-) diff --git a/src/festim/heat_transfer_problem.py b/src/festim/heat_transfer_problem.py index 268bfca90..8ad6ffdae 100644 --- a/src/festim/heat_transfer_problem.py +++ b/src/festim/heat_transfer_problem.py @@ -3,8 +3,6 @@ from dolfinx import fem from dolfinx.io import VTXWriter -import numpy as np - from festim import boundary_conditions, exports, helpers, problem from festim import source as _source @@ -113,9 +111,8 @@ def create_dirichletbc_form(self, bc): dolfinx.fem.bcs.DirichletBC: A representation of the boundary condition for modifying linear systems. """ - bc.value.convert_input_value( + bc.create_value( function_space=self.function_space, - mesh=self.mesh.mesh, t=self.t, ) @@ -124,11 +121,11 @@ def create_dirichletbc_form(self, bc): function_space=self.function_space, ) - if isinstance(bc.value.fenics_object, (fem.Function)): - form = fem.dirichletbc(value=bc.value.fenics_object, dofs=bc_dofs) + if isinstance(bc.value_fenics, (fem.Function)): + form = fem.dirichletbc(value=bc.value_fenics, dofs=bc_dofs) else: form = fem.dirichletbc( - value=bc.value.fenics_object, dofs=bc_dofs, V=self.function_space + value=bc.value_fenics, dofs=bc_dofs, V=self.function_space ) return form @@ -149,9 +146,9 @@ def create_flux_values_fenics(self): # create value_fenics for all F.HeatFluxBC objects if isinstance(bc, boundary_conditions.HeatFluxBC): - bc.value.convert_input_value( - function_space=self.function_space, + bc.create_value_fenics( mesh=self.mesh.mesh, + temperature=self.u, t=self.t, ) @@ -167,24 +164,16 @@ def create_initial_conditions(self): "Initial conditions can only be defined for transient simulations" ) - if isinstance(self.initial_condition.value.input_value, (int, float)): - self.initial_condition.value.fenics_interpolation_expression = ( - lambda x: np.full(x.shape[1], self.initial_condition.value.input_value) - ) - else: - self.initial_condition.value.fenics_interpolation_expression, _ = ( - helpers.as_fenics_interp_expr_and_function( - value=self.initial_condition.value.input_value, - function_space=self.function_space, - mesh=self.mesh.mesh, - ) - ) + # create value_fenics for condition - # assign to previous solution of species - self.u_n.interpolate( - self.initial_condition.value.fenics_interpolation_expression + self.initial_condition.create_expr_fenics( + mesh=self.mesh.mesh, + function_space=self.function_space, ) + # assign to previous solution of species + self.u_n.interpolate(self.initial_condition.expr_fenics) + def create_formulation(self): """Creates the formulation of the model""" @@ -227,9 +216,7 @@ def create_formulation(self): for bc in self.boundary_conditions: if isinstance(bc, boundary_conditions.HeatFluxBC): self.formulation -= ( - bc.value.fenics_object - * self.test_function - * self.ds(bc.subdomain.id) + bc.value_fenics * self.test_function * self.ds(bc.subdomain.id) ) def initialise_exports(self): diff --git a/src/festim/problem.py b/src/festim/problem.py index 60dd512b9..e40cc5638 100644 --- a/src/festim/problem.py +++ b/src/festim/problem.py @@ -189,8 +189,8 @@ def iterate(self): def update_time_dependent_values(self): t = float(self.t) for bc in self.boundary_conditions: - if bc.value.time_dependent: - bc.value.update(t=t) + if bc.time_dependent: + bc.update(t=t) for source in self.sources: if source.value.time_dependent: diff --git a/test/test_heat_transfer_problem.py b/test/test_heat_transfer_problem.py index 336da634c..885337f69 100644 --- a/test/test_heat_transfer_problem.py +++ b/test/test_heat_transfer_problem.py @@ -638,8 +638,6 @@ def test_update_time_dependent_values_HeatFluxBC(bc_value, expected_values): my_model.update_time_dependent_values() # TEST - if isinstance( - my_model.boundary_conditions[0].value.fenics_object, fem.Constant - ): - computed_value = float(my_model.boundary_conditions[0].value.fenics_object) + if isinstance(my_model.boundary_conditions[0].value_fenics, fem.Constant): + computed_value = float(my_model.boundary_conditions[0].value_fenics) assert np.isclose(computed_value, expected_values[i]) From abb58a45f82e7f83f758e71526cc55d9651834fa Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 11:04:33 -0500 Subject: [PATCH 14/40] revert temperature --- src/festim/hydrogen_transport_problem.py | 31 +++++++++++++++++++----- 1 file changed, 25 insertions(+), 6 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 879796e5a..0d6921ab3 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -165,6 +165,7 @@ def __init__( self.reactions = reactions or [] self.initial_conditions = initial_conditions or [] self.traps = traps or [] + self.temperature_fenics = None self._vtxfiles: list[dolfinx.io.VTXWriter] = [] @property @@ -176,23 +177,41 @@ def temperature(self, value): if value is None: self._temperature = value elif isinstance(value, (float, int, fem.Constant, fem.Function)): - self._temperature = festim.Value(value) + self._temperature = value elif callable(value): - arguments = value.__code__.co_varnames - if "T" in arguments: - raise ValueError("Temperature cannot be a function of temperature, T") - self._temperature = festim.Value(value) + self._temperature = value else: raise TypeError( "Value must be a float, int, fem.Constant, fem.Function, or callable" ) + @property + def temperature_fenics(self): + return self._temperature_fenics + + @temperature_fenics.setter + def temperature_fenics(self, value): + if value is None: + self._temperature_fenics = value + return + elif not isinstance( + value, + (fem.Constant, fem.Function), + ): + raise TypeError("Value must be a fem.Constant or fem.Function") + self._temperature_fenics = value + @property def temperature_time_dependent(self): if self.temperature is None: return False + if isinstance(self.temperature, fem.Constant): + return False + if callable(self.temperature): + arguments = self.temperature.__code__.co_varnames + return "t" in arguments else: - return self.temperature.time_dependent + return False @property def multispecies(self): From dac66e1c0941be146ee6d7abbb44b09cc9b3c081 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 11:17:45 -0500 Subject: [PATCH 15/40] revert hydrogen transport problem --- src/festim/hydrogen_transport_problem.py | 136 ++++++++++++----------- 1 file changed, 74 insertions(+), 62 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 0d6921ab3..01df1dd0e 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -313,25 +313,54 @@ def define_temperature(self): if self.temperature is None: raise ValueError("the temperature attribute needs to be defined") - if self.temperature.temperature_dependent: - raise ValueError( - "the temperature input value cannot be dependent on temperature" + # if temperature is a float or int, create a fem.Constant + elif isinstance(self.temperature, (float, int)): + self.temperature_fenics = as_fenics_constant( + self.temperature, self.mesh.mesh ) + # if temperature is a fem.Constant or function, pass it to temperature_fenics + elif isinstance(self.temperature, (fem.Constant, fem.Function)): + self.temperature_fenics = self.temperature - degree = 1 - element_temperature = basix.ufl.element( - basix.ElementFamily.P, - self.mesh.mesh.basix_cell(), - degree, - basix.LagrangeVariant.equispaced, - ) - function_space_temperature = fem.functionspace( - self.mesh.mesh, element_temperature - ) - - self.temperature.convert_input_value( - mesh=self.mesh.mesh, function_space=function_space_temperature, t=self.t - ) + # if temperature is callable, process accordingly + elif callable(self.temperature): + arguments = self.temperature.__code__.co_varnames + if "t" in arguments and "x" not in arguments: + if not isinstance(self.temperature(t=float(self.t)), (float, int)): + raise ValueError( + f"self.temperature should return a float or an int, not " + f"{type(self.temperature(t=float(self.t)))} " + ) + # only t is an argument + self.temperature_fenics = as_fenics_constant( + mesh=self.mesh.mesh, value=self.temperature(t=float(self.t)) + ) + else: + x = ufl.SpatialCoordinate(self.mesh.mesh) + degree = 1 + element_temperature = basix.ufl.element( + basix.ElementFamily.P, + self.mesh.mesh.basix_cell(), + degree, + basix.LagrangeVariant.equispaced, + ) + function_space_temperature = fem.functionspace( + self.mesh.mesh, element_temperature + ) + self.temperature_fenics = fem.Function(function_space_temperature) + kwargs = {} + if "t" in arguments: + kwargs["t"] = self.t + if "x" in arguments: + kwargs["x"] = x + + # store the expression of the temperature + # to update the temperature_fenics later + self.temperature_expr = fem.Expression( + self.temperature(**kwargs), + function_space_temperature.element.interpolation_points(), + ) + self.temperature_fenics.interpolate(self.temperature_expr) def initialise_exports(self): """Defines the export writers of the model, if field is given as @@ -415,9 +444,7 @@ def define_D_global(self, species): D = fem.Function(self.V_DG_1) expr = D_0 * ufl.exp( - -E_D - / as_fenics_constant(k_B, self.mesh.mesh) - / self.temperature.fenics_object + -E_D / as_fenics_constant(k_B, self.mesh.mesh) / self.temperature_fenics ) D_expr = fem.Expression(expr, self.V_DG_1.element.interpolation_points()) D.interpolate(D_expr) @@ -512,12 +539,10 @@ def define_boundary_conditions(self): # 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.value.convert_input_value( + bc.create_value_fenics( mesh=self.mesh.mesh, + temperature=self.temperature_fenics, t=self.t, - temperature=self.temperature.fenics_object, - up_to_ufl_expr=True, ) super().define_boundary_conditions() @@ -538,17 +563,14 @@ def create_dirichletbc_form(self, bc): else: function_space_value = bc.species.collapsed_function_space - bc.value.convert_input_value( - mesh=self.mesh.mesh, - temperature=self.temperature.fenics_object, + bc.create_value_fenics( + temperature=self.temperature_fenics, function_space=function_space_value, t=self.t, ) # get dofs - if self.multispecies and isinstance( - bc.value_fenics.fenics_object, (fem.Function) - ): + if self.multispecies and isinstance(bc.value_fenics, (fem.Function)): function_space_dofs = ( bc.species.sub_function_space, bc.species.collapsed_function_space, @@ -562,14 +584,14 @@ def create_dirichletbc_form(self, bc): ) # create form - if not self.multispecies and isinstance(bc.value.fenics_object, (fem.Function)): + if not self.multispecies and isinstance(bc.value_fenics, (fem.Function)): # no need to pass the functionspace since value_fenics is already a Function function_space_form = None else: function_space_form = bc.species.sub_function_space form = fem.dirichletbc( - value=bc.value.fenics_object, + value=bc.value_fenics, dofs=bc_dofs, V=function_space_form, ) @@ -583,7 +605,7 @@ def create_source_values_fenics(self): source.value.convert_input_value( mesh=self.mesh.mesh, t=self.t, - temperature=self.temperature.fenics_object, + temperature=self.temperature_fenics, up_to_ufl_expr=True, ) @@ -592,11 +614,10 @@ def create_flux_values_fenics(self): for bc in self.boundary_conditions: # create value_fenics for all F.ParticleFluxBC objects if isinstance(bc, boundary_conditions.ParticleFluxBC): - bc.value.convert_input_value( + bc.create_value_fenics( mesh=self.mesh.mesh, - temperature=self.temperature.fenics_object, + temperature=self.temperature_fenics, t=self.t, - up_to_ufl_expr=True, ) def create_initial_conditions(self): @@ -614,7 +635,7 @@ def create_initial_conditions(self): # create value_fenics for condition function_space_value = None - if callable(condition.value.input_value): + if callable(condition.value): if condition.value.time_dependent: raise ValueError("Initial conditions cannot be time dependent") @@ -625,30 +646,18 @@ def create_initial_conditions(self): else: function_space_value = condition.species.collapsed_function_space - if isinstance(condition.value.input_value, (int, float)): - condition.value.fenics_interpolation_expression = lambda x: np.full( - x.shape[1], condition.value.input_value - ) - else: - condition.value.fenics_interpolation_expression, _ = ( - festim.as_fenics_interp_expr_and_function( - value=condition.value.input_value, - function_space=function_space_value, - mesh=self.mesh.mesh, - temperature=self.temperature.fenics_object, - ) - ) + condition.create_expr_fenics( + mesh=self.mesh.mesh, + temperature=self.temperature_fenics, + function_space=function_space_value, + ) # assign to previous solution of species if not self.multispecies: - condition.species.prev_solution.interpolate( - condition.value.fenics_interpolation_expression - ) + condition.species.prev_solution.interpolate(condition.expr_fenics) else: idx = self.species.index(condition.species) - self.u_n.sub(idx).interpolate( - condition.value.fenics_interpolation_expression - ) + self.u_n.sub(idx).interpolate(condition.expr_fenics) def create_formulation(self): """Creates the formulation of the model""" @@ -663,7 +672,7 @@ def create_formulation(self): for vol in self.volume_subdomains: D = vol.material.get_diffusion_coefficient( - self.mesh.mesh, self.temperature.fenics_object, spe + self.mesh.mesh, self.temperature_fenics, spe ) if spe.mobile: self.formulation += ufl.dot(D * ufl.grad(u), ufl.grad(v)) * self.dx( @@ -677,7 +686,7 @@ def create_formulation(self): for reactant in reaction.reactant: if isinstance(reactant, festim.species.Species): self.formulation += ( - reaction.reaction_term(self.temperature.fenics_object) + reaction.reaction_term(self.temperature_fenics) * reactant.test_function * self.dx(reaction.volume.id) ) @@ -689,7 +698,7 @@ def create_formulation(self): products = [reaction.product] for product in products: self.formulation += ( - -reaction.reaction_term(self.temperature.fenics_object) + -reaction.reaction_term(self.temperature_fenics) * product.test_function * self.dx(reaction.volume.id) ) @@ -748,7 +757,10 @@ def update_time_dependent_values(self): if not self.temperature_time_dependent: return - self.temperature.update(t=t) + if isinstance(self.temperature_fenics, fem.Constant): + self.temperature_fenics.value = self.temperature(t=t) + elif isinstance(self.temperature_fenics, fem.Function): + self.temperature_fenics.interpolate(self.temperature_expr) for bc in self.boundary_conditions: if isinstance( @@ -758,8 +770,8 @@ def update_time_dependent_values(self): boundary_conditions.ParticleFluxBC, ), ): - if bc.value.temperature_dependent: - bc.value.update(t=t) + if bc.temperature_dependent: + bc.update(t=t) for source in self.sources: if source.value.temperature_dependent: From eb7d1a0863769d72282daf3bfb65c873640ec438 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 11:28:05 -0500 Subject: [PATCH 16/40] revert h transport problem and tests --- src/festim/hydrogen_transport_problem.py | 6 ++-- test/test_h_transport_problem.py | 38 +++++++++--------------- 2 files changed, 17 insertions(+), 27 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 01df1dd0e..6b3087db9 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -563,7 +563,7 @@ def create_dirichletbc_form(self, bc): else: function_space_value = bc.species.collapsed_function_space - bc.create_value_fenics( + bc.create_value( temperature=self.temperature_fenics, function_space=function_space_value, t=self.t, @@ -714,14 +714,14 @@ def create_formulation(self): for bc in self.boundary_conditions: if isinstance(bc, boundary_conditions.ParticleFluxBC): self.formulation -= ( - bc.value.fenics_object + 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.fenics_object + flux_bc.value_fenics * flux_bc.species.test_function * self.ds(flux_bc.subdomain.id) ) diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 6323114d8..51eb6f253 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -107,7 +107,7 @@ def test_define_temperature(input, expected_type): my_model.define_temperature() # TEST - assert isinstance(my_model.temperature.fenics_object, expected_type) + assert isinstance(my_model.temperature_fenics, expected_type) @pytest.mark.parametrize( @@ -129,7 +129,7 @@ def test_define_temperature_error_if_ufl_conditional_t_only(input): with pytest.raises( ValueError, - match="self.value should return a float or an int, not", + match="self.temperature should return a float or an int, not", ): my_model.define_temperature() @@ -214,8 +214,8 @@ def test_update_time_dependent_values_temperature(T_function, expected_values): my_model.update_time_dependent_values() # TEST - if isinstance(my_model.temperature.fenics_object, fem.Constant): - computed_value = float(my_model.temperature.fenics_object) + if isinstance(my_model.temperature_fenics, fem.Constant): + computed_value = float(my_model.temperature_fenics) assert np.isclose(computed_value, expected_values[i]) @@ -301,12 +301,8 @@ def test_define_D_global_different_materials(): computed_values = [D_computed.x.array[0], D_computed.x.array[-1]] - D_expected_left = D_0_left * np.exp( - -E_D_left / (F.k_B * my_model.temperature.input_value) - ) - D_expected_right = D_0_right * np.exp( - -E_D_right / (F.k_B * my_model.temperature.input_value) - ) + D_expected_left = D_0_left * np.exp(-E_D_left / (F.k_B * my_model.temperature)) + D_expected_right = D_0_right * np.exp(-E_D_right / (F.k_B * my_model.temperature)) expected_values = [D_expected_left, D_expected_right] @@ -417,8 +413,8 @@ def test_define_D_global_multispecies(): computed_values = [D_A_computed.x.array[-1], D_B_computed.x.array[-1]] - D_expected_A = D_0_A * np.exp(-E_D_A / (F.k_B * my_model.temperature.input_value)) - D_expected_B = D_0_B * np.exp(-E_D_B / (F.k_B * my_model.temperature.input_value)) + D_expected_A = D_0_A * np.exp(-E_D_A / (F.k_B * my_model.temperature)) + D_expected_B = D_0_B * np.exp(-E_D_B / (F.k_B * my_model.temperature)) expected_values = [D_expected_A, D_expected_B] @@ -640,10 +636,8 @@ def test_update_time_dependent_bcs_with_time_dependent_temperature( my_model.update_time_dependent_values() # TEST - if isinstance( - my_model.boundary_conditions[0].value.fenics_object, fem.Constant - ): - computed_value = float(my_model.boundary_conditions[0].value.fenics_object) + if isinstance(my_model.boundary_conditions[0].value_fenics, fem.Constant): + computed_value = float(my_model.boundary_conditions[0].value_fenics) assert np.isclose(computed_value, expected_values[i]) @@ -1152,10 +1146,8 @@ def test_update_time_dependent_values_flux(bc_value, expected_values): my_model.update_time_dependent_values() # TEST - if isinstance( - my_model.boundary_conditions[0].value.fenics_object, fem.Constant - ): - computed_value = float(my_model.boundary_conditions[0].value.fenics_object) + if isinstance(my_model.boundary_conditions[0].value_fenics, fem.Constant): + computed_value = float(my_model.boundary_conditions[0].value_fenics) assert np.isclose(computed_value, expected_values[i]) @@ -1208,10 +1200,8 @@ def test_update_fluxes_with_time_dependent_temperature( my_model.update_time_dependent_values() # TEST - if isinstance( - my_model.boundary_conditions[0].value.fenics_object, fem.Constant - ): - computed_value = float(my_model.boundary_conditions[0].value.fenics_object) + if isinstance(my_model.boundary_conditions[0].value_fenics, fem.Constant): + computed_value = float(my_model.boundary_conditions[0].value_fenics) assert np.isclose(computed_value, expected_values[i]) From cf7cecb141920e74fd37d3c415bbdf69d2d1d5b2 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 11:30:47 -0500 Subject: [PATCH 17/40] revert changes to just source --- src/festim/hydrogen_transport_problem.py | 5 ----- test/test_h_transport_problem.py | 4 ++-- 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 6b3087db9..da5a13a57 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -632,14 +632,9 @@ def create_initial_conditions(self): function_space_value = None for condition in self.initial_conditions: - # create value_fenics for condition function_space_value = None if callable(condition.value): - - if condition.value.time_dependent: - raise ValueError("Initial conditions cannot be time dependent") - # if bc.value is a callable then need to provide a functionspace if not self.multispecies: function_space_value = condition.species.sub_function_space diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 51eb6f253..200420221 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -1233,5 +1233,5 @@ def test_create_flux_values_fenics_multispecies(): my_model.create_flux_values_fenics() # TEST - assert np.isclose(float(my_model.boundary_conditions[0].value.fenics_object), 5) - assert np.isclose(float(my_model.boundary_conditions[1].value.fenics_object), 11) + assert np.isclose(float(my_model.boundary_conditions[0].value_fenics), 5) + assert np.isclose(float(my_model.boundary_conditions[1].value_fenics), 11) From 41396f379b2e88e98e5463b1933964419e8dd842 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 13:00:56 -0500 Subject: [PATCH 18/40] use fenics object --- src/festim/hydrogen_transport_problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index da5a13a57..b4c3fd31b 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1109,7 +1109,7 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain): for source in self.sources: v = source.species.subdomain_to_test_function[subdomain] if source.volume == subdomain: - form -= source.value_fenics * v * self.dx(subdomain.id) + form -= source.value.fenics_object * v * self.dx(subdomain.id) # store the form in the subdomain object subdomain.F = form From 77a8422ac71009e85033f14ed14f1afb7fd1819f Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 13:01:15 -0500 Subject: [PATCH 19/40] accept None as input for value --- src/festim/source.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/festim/source.py b/src/festim/source.py index 8028cc388..0abe02196 100644 --- a/src/festim/source.py +++ b/src/festim/source.py @@ -1,5 +1,5 @@ from dolfinx import fem - +import ufl import festim as F @@ -44,8 +44,11 @@ def value(self): @value.setter def value(self, value): if value is None: - self._value = value - elif isinstance(value, (float, int, fem.Constant, fem.Function)): + self._value = F.Value(value) + elif isinstance( + value, + (float, int, fem.Constant, fem.Function, ufl.core.expr.Expr), + ): self._value = F.Value(value) elif callable(value): self._value = F.Value(value) From 83dfca51a49317e3ecd8702c949e2a88ad7d9afb Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 13:01:30 -0500 Subject: [PATCH 20/40] fix tests --- src/festim/helpers.py | 8 ++++++-- test/test_source.py | 46 ++++++++++++------------------------------- 2 files changed, 19 insertions(+), 35 deletions(-) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 0c39df50d..598d0543f 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -177,6 +177,8 @@ def time_dependent(self) -> bool: return False if isinstance(self.input_value, fem.Constant): return False + if isinstance(self.input_value, ufl.core.expr.Expr): + return False if callable(self.input_value): arguments = self.input_value.__code__.co_varnames return "t" in arguments @@ -190,6 +192,8 @@ def temperature_dependent(self) -> bool: return False if isinstance(self.input_value, fem.Constant): return False + if isinstance(self.input_value, ufl.core.expr.Expr): + return False if callable(self.input_value): arguments = self.input_value.__code__.co_varnames return "T" in arguments @@ -262,11 +266,11 @@ def convert_input_value( f"Value must be a float, an int or a callable, not {type(self.input_value)}" ) - def update(self, t): + def update(self, t: float): """Updates the value Args: - t (float): the time + t: the time """ if callable(self.input_value): arguments = self.input_value.__code__.co_varnames diff --git a/test/test_source.py b/test/test_source.py index 9bddca94f..411d87eb5 100644 --- a/test/test_source.py +++ b/test/test_source.py @@ -22,32 +22,11 @@ def test_init(): # check that the attributes are set correctly assert source.volume == volume - assert source.value == value assert source.species == species - assert source.value_fenics is None - assert source.source_expr is None - -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 Source object - volume = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=dummy_mat) - value = 1.0 - species = F.Species("test") - source = F.ParticleSource(volume=volume, value=value, species=species) - - # set the value_fenics attribute to a valid value - value_fenics = fem.Constant(mesh, 2.0) - source.value_fenics = value_fenics - - # check that the value_fenics attribute is set correctly - assert source.value_fenics == value_fenics - - # set the value_fenics attribute to an invalid value - with pytest.raises(TypeError): - source.value_fenics = "invalid" + # check value is processed correctly + assert source.value.input_value == value + assert isinstance(source.value, F.Value) @pytest.mark.parametrize( @@ -66,9 +45,8 @@ def test_value_fenics(): ), ], ) -def test_create_value_fenics(value, expected_type): - """Test that the create value method produces either a fem.Constant or - fem.Function depending on the value input""" +def test_create_fenics_object(value, expected_type): + """Test that the correct fenics object is created depending on the value input""" # BUILD vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) @@ -80,11 +58,11 @@ def test_create_value_fenics(value, expected_type): t = fem.Constant(mesh, 0.0) # RUN - source.create_value_fenics(mesh, T, t) + source.value.convert_input_value(mesh=mesh, temperature=T, t=t, up_to_ufl_expr=True) # TEST # check that the value_fenics attribute is set correctly - assert isinstance(source.value_fenics, expected_type) + assert isinstance(source.value.fenics_object, expected_type) @pytest.mark.parametrize( @@ -107,7 +85,7 @@ def test_source_time_dependent_attribute(input, expected_value): species = F.Species("test") my_source = F.ParticleSource(input, volume, species) - assert my_source.time_dependent is expected_value + assert my_source.value.time_dependent is expected_value @pytest.mark.parametrize( @@ -132,7 +110,7 @@ def test_source_temperature_dependent_attribute(input, expected_value): species = F.Species("test") my_source = F.ParticleSource(input, volume, species) - assert my_source.temperature_dependent is expected_value + assert my_source.value.temperature_dependent is expected_value def test_ValueError_raised_when_callable_returns_wrong_type(): @@ -154,7 +132,9 @@ def my_value(t): ValueError, match="self.value should return a float or an int, not Date: Thu, 20 Feb 2025 14:30:29 -0500 Subject: [PATCH 21/40] format ruff --- src/festim/heat_transfer_problem.py | 1 - src/festim/helpers.py | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/festim/heat_transfer_problem.py b/src/festim/heat_transfer_problem.py index 8ad6ffdae..477ef815f 100644 --- a/src/festim/heat_transfer_problem.py +++ b/src/festim/heat_transfer_problem.py @@ -145,7 +145,6 @@ def create_flux_values_fenics(self): for bc in self.boundary_conditions: # create value_fenics for all F.HeatFluxBC objects if isinstance(bc, boundary_conditions.HeatFluxBC): - bc.create_value_fenics( mesh=self.mesh.mesh, temperature=self.u, diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 598d0543f..2a2220ebe 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -132,7 +132,7 @@ class Value: fenics_object: fem.Function | fem.Constant | ufl.core.expr.Expr def __init__(self, input_value): - self, + (self,) self.input_value = input_value self.ufl_expression = None From 5af11cc0eeb16654f62c6b361ea7bce27bf34fa2 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 14:35:44 -0500 Subject: [PATCH 22/40] use fenics object --- src/festim/problem_change_of_var.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/festim/problem_change_of_var.py b/src/festim/problem_change_of_var.py index c0771c2c2..4f39285ec 100644 --- a/src/festim/problem_change_of_var.py +++ b/src/festim/problem_change_of_var.py @@ -104,7 +104,7 @@ def create_formulation(self): # add sources for source in self.sources: self.formulation -= ( - source.value_fenics + source.value.fenics_object * source.species.test_function * self.dx(source.volume.id) ) From 2cca0b1676f12fcedb07a790206c677199cf8c1c Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 14:56:26 -0500 Subject: [PATCH 23/40] format ruff --- src/festim/helpers.py | 50 +++++++++++++++++++++---------------------- 1 file changed, 24 insertions(+), 26 deletions(-) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 2a2220ebe..1178037df 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -1,8 +1,9 @@ +from collections.abc import Callable + import dolfinx -from dolfinx import fem import numpy as np -from collections.abc import Callable import ufl +from dolfinx import fem def as_fenics_constant( @@ -20,7 +21,7 @@ def as_fenics_constant( Raises: TypeError: if the value is not a float, an int or a dolfinx.Constant """ - if isinstance(value, (float, int)): + if isinstance(value, "float | int"): return fem.Constant(mesh, dolfinx.default_scalar_type(float(value))) elif isinstance(value, fem.Constant): return value @@ -36,7 +37,8 @@ def as_mapped_function( t: fem.Constant = None, temperature: fem.Function | fem.Constant | ufl.core.expr.Expr = None, ) -> ufl.core.expr.Expr: - """Maps a user given callable function to the mesh, time or temperature within festim as needed + """Maps a user given callable function to the mesh, time or temperature within + festim as needed Args: value: the callable to convert @@ -70,8 +72,8 @@ def as_fenics_interp_expr_and_function( temperature: fem.Function | fem.Constant | ufl.core.expr.Expr = None, ) -> tuple[fem.Expression, fem.Function]: """Takes a user given callable function, maps the function to the mesh, time or - temperature within festim as needed. Then creates the fenics interpolation expression - and function objects + temperature within festim as needed. Then creates the fenics interpolation + expression and function objects Args: value: the callable to convert @@ -101,15 +103,16 @@ def as_fenics_interp_expr_and_function( class Value: """ - A class to handle input values from users and convert them to a relevent fenics object + A class to handle input values from users and convert them to a relevent fenics + object Args: input_value: The value of the user input Attributes: input_value : The value of the user input - fenics_interpolation_expression : The expression of the user input that is used to - update the `fenics_object` + fenics_interpolation_expression : The expression of the user input that is used + to update the `fenics_object` fenics_object : The value of the user input in fenics format """ @@ -152,15 +155,7 @@ def input_value(self, value): self._input_value = value elif isinstance( value, - ( - float, - int, - fem.Constant, - np.ndarray, - fem.Function, - ufl.core.expr.Expr, - fem.Function, - ), + "float | int | fem.Constant | np.ndarray | fem.Function | ufl.core.expr.Expr | fem.Function", # noqa: E501 ): self._input_value = value elif callable(value): @@ -213,11 +208,13 @@ def convert_input_value( Args: mesh (dolfinx.mesh.Mesh): the mesh of the domain - function_space (dolfinx.fem.function.FunctionSpace): the function space of the fenics object + function_space (dolfinx.fem.function.FunctionSpace): the function space of + the fenics object t (fem.Constant): the time - temperature (fem.Function, fem.Constant or ufl.core.expr.Expr): the temperature - up_to_ufl_expr (bool): if True, the value is only mapped to a function if the input is callable, - not interpolated or converted to a function + temperature (fem.Function, fem.Constant or ufl.core.expr.Expr): the + temperature + up_to_ufl_expr (bool): if True, the value is only mapped to a function if + the input is callable, not interpolated or converted to a function """ if isinstance(self.input_value, fem.Constant): self.fenics_object = self.input_value @@ -225,17 +222,17 @@ def convert_input_value( elif isinstance(self.input_value, fem.Expression): self.fenics_interpolation_expression = self.input_value - elif isinstance(self.input_value, (fem.Function, ufl.core.expr.Expr)): + elif isinstance(self.input_value, fem.Function | ufl.core.expr.Expr): self.fenics_object = self.input_value - elif isinstance(self.input_value, (float, int)): + elif isinstance(self.input_value, float | int): self.fenics_object = as_fenics_constant(value=self.input_value, mesh=mesh) elif callable(self.input_value): args = self.input_value.__code__.co_varnames # if only t is an argument, create constant object if "t" in args and "x" not in args and "T" not in args: - if not isinstance(self.input_value(t=float(t)), (float, int)): + if not isinstance(self.input_value(t=float(t)), float | int): raise ValueError( "self.value should return a float or an int, not " + f"{type(self.input_value(t=float(t)))} " @@ -263,7 +260,8 @@ def convert_input_value( else: raise TypeError( - f"Value must be a float, an int or a callable, not {type(self.input_value)}" + f"Value must be a float, an int or a callable, not " + f"{type(self.input_value)}" ) def update(self, t: float): From 9e4e8db6869e7309bda19fc49356189ebadd24a5 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 16:03:28 -0500 Subject: [PATCH 24/40] formatted ruff --- src/festim/helpers.py | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 1178037df..d26001590 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -21,7 +21,7 @@ def as_fenics_constant( Raises: TypeError: if the value is not a float, an int or a dolfinx.Constant """ - if isinstance(value, "float | int"): + if isinstance(value, float | int): return fem.Constant(mesh, dolfinx.default_scalar_type(float(value))) elif isinstance(value, fem.Constant): return value @@ -155,7 +155,13 @@ def input_value(self, value): self._input_value = value elif isinstance( value, - "float | int | fem.Constant | np.ndarray | fem.Function | ufl.core.expr.Expr | fem.Function", # noqa: E501 + float + | int + | fem.Constant + | np.ndarray + | fem.Function + | ufl.core.expr.Expr + | fem.Function, ): self._input_value = value elif callable(value): @@ -170,9 +176,7 @@ def time_dependent(self) -> bool: """Returns true if the value given is time dependent""" if self.input_value is None: return False - if isinstance(self.input_value, fem.Constant): - return False - if isinstance(self.input_value, ufl.core.expr.Expr): + if isinstance(self.input_value, fem.Constant | ufl.core.expr.Expr): return False if callable(self.input_value): arguments = self.input_value.__code__.co_varnames @@ -185,9 +189,7 @@ def temperature_dependent(self) -> bool: """Returns true if the value given is temperature dependent""" if self.input_value is None: return False - if isinstance(self.input_value, fem.Constant): - return False - if isinstance(self.input_value, ufl.core.expr.Expr): + if isinstance(self.input_value, fem.Constant | ufl.core.expr.Expr): return False if callable(self.input_value): arguments = self.input_value.__code__.co_varnames @@ -216,15 +218,14 @@ def convert_input_value( up_to_ufl_expr (bool): if True, the value is only mapped to a function if the input is callable, not interpolated or converted to a function """ - if isinstance(self.input_value, fem.Constant): + if isinstance( + self.input_value, fem.Constant | fem.Function | ufl.core.expr.Expr + ): self.fenics_object = self.input_value elif isinstance(self.input_value, fem.Expression): self.fenics_interpolation_expression = self.input_value - elif isinstance(self.input_value, fem.Function | ufl.core.expr.Expr): - self.fenics_object = self.input_value - elif isinstance(self.input_value, float | int): self.fenics_object = as_fenics_constant(value=self.input_value, mesh=mesh) From aac0a525a89cd7c339479165a864c9ff7e21a44d Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 16:03:36 -0500 Subject: [PATCH 25/40] add tests for coverage --- test/test_helpers.py | 200 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 199 insertions(+), 1 deletion(-) diff --git a/test/test_helpers.py b/test/test_helpers.py index 9220b1451..4d21467e5 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -1,3 +1,7 @@ +from mpi4py import MPI + +import basix +import dolfinx.mesh import numpy as np import pytest import ufl @@ -8,6 +12,15 @@ test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) x = ufl.SpatialCoordinate(test_mesh.mesh) +element_CG = basix.ufl.element( + basix.ElementFamily.P, + test_mesh.mesh.basix_cell(), + 1, + basix.LagrangeVariant.equispaced, +) +test_function_space = fem.functionspace(test_mesh.mesh, element_CG) +test_function = fem.Function(test_function_space) + @pytest.mark.parametrize( "value", @@ -22,8 +35,193 @@ def test_temperature_type_and_processing(value): """Test that the temperature type is correctly set""" - if not isinstance(value, (fem.Constant, int, float)): + if not isinstance(value, fem.Constant | int | float): with pytest.raises(TypeError): F.as_fenics_constant(value, test_mesh.mesh) else: assert isinstance(F.as_fenics_constant(value, test_mesh.mesh), fem.Constant) + + +@pytest.mark.parametrize( + "input_value, expected_output_type", [(1.0, fem.Constant), (3, fem.Constant)] +) +def test_value_convert_float_int_inputs(input_value, expected_output_type): + """Test that float and value is correctly converted""" + + mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10) + + test_value = F.Value(input_value) + + test_value.convert_input_value(mesh=mesh) + + assert isinstance(test_value.fenics_object, expected_output_type) + + +@pytest.mark.parametrize( + "input_value, expected_output_type", + [ + (lambda t: t, fem.Constant), + (lambda t: 1.0 + t, fem.Constant), + (lambda x: 1.0 + x[0], ufl.core.expr.Expr), + (lambda x, t: 1.0 + x[0] + t, ufl.core.expr.Expr), + (lambda x, t, T: 1.0 + x[0] + t + T, ufl.core.expr.Expr), + ( + lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), + ufl.core.expr.Expr, + ), + ], +) +def test_value_convert_up_to_ufl_inputs(input_value, expected_output_type): + """Test that float and value is correctly converted""" + + my_mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10) + my_t = fem.Constant(my_mesh, default_scalar_type(10)) + my_T = fem.Constant(my_mesh, default_scalar_type(3)) + + test_value = F.Value(input_value) + + test_value.convert_input_value( + mesh=my_mesh, + t=my_t, + temperature=my_T, + up_to_ufl_expr=True, + ) + + assert isinstance(test_value.fenics_object, expected_output_type) + + +@pytest.mark.parametrize( + "input_value, expected_output_type", + [ + (lambda t: t, fem.Constant), + (lambda t: 1.0 + t, fem.Constant), + (lambda x: 1.0 + x[0], fem.Function), + (lambda x, t: 1.0 + x[0] + t, fem.Function), + (lambda x, t, T: 1.0 + x[0] + t + T, fem.Function), + ( + lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), + ufl.core.expr.Expr, + ), + ], +) +def test_value_convert_callable_inputs(input_value, expected_output_type): + """Test that float and value is correctly converted""" + + my_mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 12) + my_t = fem.Constant(my_mesh, default_scalar_type(8)) + my_T = fem.Constant(my_mesh, default_scalar_type(5)) + + element_CG = basix.ufl.element( + basix.ElementFamily.P, + my_mesh.basix_cell(), + 1, + basix.LagrangeVariant.equispaced, + ) + my_function_space = fem.functionspace(my_mesh, element_CG) + + test_value = F.Value(input_value) + + test_value.convert_input_value( + function_space=my_function_space, + mesh=my_mesh, + t=my_t, + temperature=my_T, + ) + + assert isinstance(test_value.fenics_object, expected_output_type) + + +def test_error_raised_wehn_input_value_is_not_accepted(): + """Test that an error is raised when the input value is not accepted""" + + with pytest.raises( + TypeError, + match="Value must be a float, int, fem.Constant, fem.Function, or callable", + ): + F.Value("coucou") + + +@pytest.mark.parametrize( + "input_value, expected_output", + [ + (1.0, False), + (None, False), + (fem.Constant(test_mesh.mesh, default_scalar_type(1.0)), False), + (lambda t: t, True), + (lambda t: 1.0 + t, True), + (lambda x: 1.0 + x[0], False), + (lambda x, t: 1.0 + x[0] + t, True), + (lambda x, t, T: 1.0 + x[0] + t + T, True), + ( + lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), + True, + ), + ], +) +def test_time_dependent_values(input_value, expected_output): + """Test that the time_dependent attribute is correctly set""" + + test_value = F.Value(input_value) + + assert test_value.time_dependent == expected_output + + +@pytest.mark.parametrize( + "input_value, expected_output", + [ + (1.0, False), + (None, False), + (fem.Constant(test_mesh.mesh, default_scalar_type(1.0)), False), + (lambda t: t, False), + (lambda T: 1.0 + T, True), + (lambda x: 1.0 + x[0], False), + (lambda x, t: 1.0 + x[0] + t, False), + (lambda x, t, T: 1.0 + x[0] + t + T, True), + ( + lambda T, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + T[0], 0.0), + True, + ), + ], +) +def test_temperature_dependent_values(input_value, expected_output): + """Test that the time_dependent attribute is correctly set""" + + test_value = F.Value(input_value) + + assert test_value.temperature_dependent == expected_output + + +@pytest.mark.parametrize( + "value", + [ + fem.Constant(test_mesh.mesh, default_scalar_type(1.0)), + test_function, + ], +) +def test_input_values_of_constants_and_functions_are_accepted(value): + """Test that the input values of constants and functions are accepted""" + + test_value = F.Value(value) + + test_value.convert_input_value() + + assert test_value.fenics_object == value + + +def test_ValueError_raised_when_callable_returns_wrong_type(): + """The create_value_fenics method should raise a ValueError when the callable + returns an object which is not a float or int""" + + def my_value(t): + return ufl.conditional(ufl.lt(t, 0.5), 100, 0) + + test_value = F.Value(my_value) + + T = fem.Constant(test_mesh.mesh, 550.0) + t = fem.Constant(test_mesh.mesh, 0.0) + + with pytest.raises( + ValueError, + match="self.value should return a float or an int, not Date: Thu, 20 Feb 2025 16:59:12 -0500 Subject: [PATCH 26/40] setter already in Value --- src/festim/source.py | 16 +--------------- 1 file changed, 1 insertion(+), 15 deletions(-) diff --git a/src/festim/source.py b/src/festim/source.py index 0abe02196..393d5aede 100644 --- a/src/festim/source.py +++ b/src/festim/source.py @@ -1,5 +1,3 @@ -from dolfinx import fem -import ufl import festim as F @@ -43,19 +41,7 @@ def value(self): @value.setter def value(self, value): - if value is None: - self._value = F.Value(value) - elif isinstance( - value, - (float, int, fem.Constant, fem.Function, ufl.core.expr.Expr), - ): - self._value = F.Value(value) - elif callable(value): - self._value = F.Value(value) - else: - raise TypeError( - "Value must be a float, int, fem.Constant, fem.Function, or callable" - ) + self._value = F.Value(value) @property def volume(self): From 6bc191b8fc1f5d76334fabe019eeab950d2a98dc Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 16:59:30 -0500 Subject: [PATCH 27/40] remove bloat, accept expression --- src/festim/helpers.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index d26001590..761a914d2 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -159,7 +159,7 @@ def input_value(self, value): | int | fem.Constant | np.ndarray - | fem.Function + | fem.Expression | ufl.core.expr.Expr | fem.Function, ): @@ -168,7 +168,8 @@ def input_value(self, value): self._input_value = value else: raise TypeError( - "Value must be a float, int, fem.Constant, fem.Function, or callable" + "Value must be a float, int, fem.Constant, np.ndarray, fem.Expression," + f" ufl.core.expr.Expr, fem.Function, or callable not {value}" ) @property @@ -259,12 +260,6 @@ def convert_input_value( ) ) - else: - raise TypeError( - f"Value must be a float, an int or a callable, not " - f"{type(self.input_value)}" - ) - def update(self, t: float): """Updates the value From f4f86d60aa6c732ae50ae35fc8c48089b504ae93 Mon Sep 17 00:00:00 2001 From: jhdark Date: Thu, 20 Feb 2025 16:59:40 -0500 Subject: [PATCH 28/40] more tests for coverage --- test/test_helpers.py | 46 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) diff --git a/test/test_helpers.py b/test/test_helpers.py index 4d21467e5..2419683da 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -136,7 +136,10 @@ def test_error_raised_wehn_input_value_is_not_accepted(): with pytest.raises( TypeError, - match="Value must be a float, int, fem.Constant, fem.Function, or callable", + match=( + "Value must be a float, int, fem.Constant, np.ndarray, fem.Expression, " + "ufl.core.expr.Expr, fem.Function, or callable not coucou" + ), ): F.Value("coucou") @@ -208,6 +211,25 @@ def test_input_values_of_constants_and_functions_are_accepted(value): assert test_value.fenics_object == value +def test_input_values_of_expressions_are_accepted(): + """Test that the input values of constants and functions are accepted""" + + my_func = lambda x: 1.0 + x[0] + kwargs = {} + kwargs["x"] = x + mapped_func = my_func(**kwargs) + + test_expression = fem.Expression( + mapped_func, + test_function_space.element.interpolation_points(), + ) + test_value = F.Value(input_value=test_expression) + + test_value.convert_input_value() + + assert test_value.fenics_interpolation_expression == test_expression + + def test_ValueError_raised_when_callable_returns_wrong_type(): """The create_value_fenics method should raise a ValueError when the callable returns an object which is not a float or int""" @@ -225,3 +247,25 @@ def my_value(t): match="self.value should return a float or an int, not Date: Fri, 21 Feb 2025 12:20:09 -0500 Subject: [PATCH 29/40] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Jørgen Schartum Dokken --- src/festim/helpers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 761a914d2..9618fddf4 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -50,6 +50,7 @@ def as_mapped_function( The mapped function """ + # Extract the input variable names in the callable function `value` arguments = value.__code__.co_varnames kwargs = {} @@ -135,7 +136,6 @@ class Value: fenics_object: fem.Function | fem.Constant | ufl.core.expr.Expr def __init__(self, input_value): - (self,) self.input_value = input_value self.ufl_expression = None From ae2763c72afd82d9696092d3b9e5a695790dac77 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 21 Feb 2025 14:23:45 -0500 Subject: [PATCH 30/40] use mesh from functionspace --- src/festim/heat_transfer_problem.py | 1 + src/festim/helpers.py | 71 ++++++++++-------------- src/festim/hydrogen_transport_problem.py | 30 ++++++++-- test/test_helpers.py | 2 + test/test_source.py | 5 +- 5 files changed, 61 insertions(+), 48 deletions(-) diff --git a/src/festim/heat_transfer_problem.py b/src/festim/heat_transfer_problem.py index 477ef815f..19cd6205e 100644 --- a/src/festim/heat_transfer_problem.py +++ b/src/festim/heat_transfer_problem.py @@ -136,6 +136,7 @@ def create_source_values_fenics(self): # create value_fenics for all source objects source.value.convert_input_value( mesh=self.mesh.mesh, + function_space=self.function_space, t=self.t, up_to_ufl_expr=True, ) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 9618fddf4..3bfec168e 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -5,6 +5,8 @@ import ufl from dolfinx import fem +from typing import Optional + def as_fenics_constant( value: float | int | fem.Constant, mesh: dolfinx.mesh.Mesh @@ -33,18 +35,18 @@ def as_fenics_constant( def as_mapped_function( value: Callable, - mesh: dolfinx.mesh.Mesh = None, - t: fem.Constant = None, - temperature: fem.Function | fem.Constant | ufl.core.expr.Expr = None, + function_space: Optional[fem.functionspace] = None, + t: Optional[fem.Constant] = None, + temperature: Optional[fem.Function | fem.Constant | ufl.core.expr.Expr] = None, ) -> ufl.core.expr.Expr: """Maps a user given callable function to the mesh, time or temperature within festim as needed Args: value: the callable to convert - mesh: the mesh of the domain - t: the time - temperature: the temperature + function_space: the function space of the domain, optional + t: the time, optional + temperature: the temperature, optional Returns: The mapped function @@ -57,7 +59,7 @@ def as_mapped_function( if "t" in arguments: kwargs["t"] = t if "x" in arguments: - x = ufl.SpatialCoordinate(mesh) + x = ufl.SpatialCoordinate(function_space.mesh) kwargs["x"] = x if "T" in arguments: kwargs["T"] = temperature @@ -68,9 +70,8 @@ def as_mapped_function( def as_fenics_interp_expr_and_function( value: Callable, function_space: dolfinx.fem.function.FunctionSpace, - mesh: dolfinx.mesh.Mesh = None, - t: fem.Constant = None, - temperature: fem.Function | fem.Constant | ufl.core.expr.Expr = None, + t: Optional[fem.Constant] = None, + temperature: Optional[fem.Function | fem.Constant | ufl.core.expr.Expr] = None, ) -> tuple[fem.Expression, fem.Function]: """Takes a user given callable function, maps the function to the mesh, time or temperature within festim as needed. Then creates the fenics interpolation @@ -79,16 +80,15 @@ def as_fenics_interp_expr_and_function( Args: value: the callable to convert function_space: The function space to interpolate function over - mesh: the mesh of the domain - t: the time - temperature: the temperature + t: the time, optional + temperature: the temperature, optional Returns: fenics interpolation expression, fenics function """ mapped_function = as_mapped_function( - value=value, mesh=mesh, t=t, temperature=temperature + value=value, function_space=function_space, t=t, temperature=temperature ) fenics_interpolation_expression = fem.Expression( @@ -118,18 +118,7 @@ class Value: """ - input_value: ( - int - | float - | np.ndarray - | Callable[[np.ndarray], np.ndarray] - | Callable[[np.ndarray, float], np.ndarray] - | Callable[[float], float] - | fem.Constant - | fem.Expression - | ufl.core.expr.Expr - | fem.Function - ) + input_value: float | int | fem.Constant | np.ndarray | fem.Expression | ufl.core.expr.Expr | fem.Function ufl_expression: ufl.core.expr.Expr fenics_interpolation_expression: fem.Expression @@ -200,24 +189,22 @@ def temperature_dependent(self) -> bool: def convert_input_value( self, - function_space: dolfinx.fem.function.FunctionSpace = None, - mesh: dolfinx.mesh.Mesh = None, - t: fem.Constant = None, - temperature: fem.Function | fem.Constant | ufl.core.expr.Expr = None, - up_to_ufl_expr: bool = False, + mesh: Optional[dolfinx.mesh.Mesh] = None, + function_space: Optional[dolfinx.fem.function.FunctionSpace] = None, + t: Optional[fem.Constant] = None, + temperature: Optional[fem.Function | fem.Constant | ufl.core.expr.Expr] = None, + up_to_ufl_expr: Optional[bool] = False, ): """Converts a user given value to a relevent fenics object depending on the type of the value provided Args: - mesh (dolfinx.mesh.Mesh): the mesh of the domain - function_space (dolfinx.fem.function.FunctionSpace): the function space of - the fenics object - t (fem.Constant): the time - temperature (fem.Function, fem.Constant or ufl.core.expr.Expr): the - temperature - up_to_ufl_expr (bool): if True, the value is only mapped to a function if - the input is callable, not interpolated or converted to a function + mesh: the mesh of the domain, optional + function_space: the function space of the fenics object, optional + t: the time, optional + temperature: the temperature, optional + up_to_ufl_expr: if True, the value is only mapped to a function if the input + is callable, not interpolated or converted to a function, optional """ if isinstance( self.input_value, fem.Constant | fem.Function | ufl.core.expr.Expr @@ -246,7 +233,10 @@ def convert_input_value( elif up_to_ufl_expr: self.fenics_object = as_mapped_function( - value=self.input_value, mesh=mesh, t=t, temperature=temperature + value=self.input_value, + function_space=function_space, + t=t, + temperature=temperature, ) else: @@ -254,7 +244,6 @@ def convert_input_value( as_fenics_interp_expr_and_function( value=self.input_value, function_space=function_space, - mesh=mesh, t=t, temperature=temperature, ) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index dc7eaed23..0627304c6 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -29,6 +29,7 @@ ) from festim.helpers import as_fenics_constant from festim.mesh import Mesh +from festim import source as _source __all__ = ["HydrogenTransportProblem", "HTransportProblemDiscontinuous"] @@ -615,12 +616,14 @@ def create_source_values_fenics(self): """For each source create the value_fenics""" for source in self.sources: # create value_fenics for all F.ParticleSource objects - source.value.convert_input_value( - mesh=self.mesh.mesh, - t=self.t, - temperature=self.temperature_fenics, - up_to_ufl_expr=True, - ) + if isinstance(source, _source.ParticleSource): + source.value.convert_input_value( + mesh=self.mesh.mesh, + function_space=self.function_space, + t=self.t, + temperature=self.temperature_fenics, + up_to_ufl_expr=True, + ) def create_flux_values_fenics(self): """For each particle flux create the value_fenics""" @@ -1054,6 +1057,21 @@ def define_function_spaces(self, subdomain: _subdomain.VolumeSubdomain): name = f"{species.name}_{subdomain.id}" species.subdomain_to_post_processing_solution[subdomain].name = name + def create_source_values_fenics(self): + """For each source create the value_fenics""" + for source in self.sources: + # create value_fenics for all F.ParticleSource objects + if isinstance(source, _source.ParticleSource): + V = dolfinx.fem.functionspace(self.mesh.mesh, ("Lagrange", 1)) + + source.value.convert_input_value( + mesh=self.mesh.mesh, + function_space=V, + t=self.t, + temperature=self.temperature_fenics, + up_to_ufl_expr=True, + ) + def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain): """ Creates the variational formulation for each subdomain and stores it in ``subdomain.F`` diff --git a/test/test_helpers.py b/test/test_helpers.py index 2419683da..ab6697f7a 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -75,6 +75,7 @@ def test_value_convert_up_to_ufl_inputs(input_value, expected_output_type): """Test that float and value is correctly converted""" my_mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10) + V = fem.functionspace(my_mesh, ("Lagrange", 1)) my_t = fem.Constant(my_mesh, default_scalar_type(10)) my_T = fem.Constant(my_mesh, default_scalar_type(3)) @@ -82,6 +83,7 @@ def test_value_convert_up_to_ufl_inputs(input_value, expected_output_type): test_value.convert_input_value( mesh=my_mesh, + function_space=V, t=my_t, temperature=my_T, up_to_ufl_expr=True, diff --git a/test/test_source.py b/test/test_source.py index 411d87eb5..9ec7ab5ee 100644 --- a/test/test_source.py +++ b/test/test_source.py @@ -51,6 +51,7 @@ def test_create_fenics_object(value, expected_type): # BUILD vol_subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) species = F.Species("test") + V = fem.functionspace(mesh, ("Lagrange", 1)) source = F.ParticleSource(volume=vol_subdomain, value=value, species=species) @@ -58,7 +59,9 @@ def test_create_fenics_object(value, expected_type): t = fem.Constant(mesh, 0.0) # RUN - source.value.convert_input_value(mesh=mesh, temperature=T, t=t, up_to_ufl_expr=True) + source.value.convert_input_value( + function_space=V, mesh=mesh, temperature=T, t=t, up_to_ufl_expr=True + ) # TEST # check that the value_fenics attribute is set correctly From a0a10a6274291d6e516210bedf1725dd8ee436fc Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 21 Feb 2025 14:35:13 -0500 Subject: [PATCH 31/40] format ruff --- src/festim/helpers.py | 3 +-- src/festim/hydrogen_transport_problem.py | 6 +++--- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 3bfec168e..99b6d8b2c 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -1,12 +1,11 @@ from collections.abc import Callable +from typing import Optional import dolfinx import numpy as np import ufl from dolfinx import fem -from typing import Optional - def as_fenics_constant( value: float | int | fem.Constant, mesh: dolfinx.mesh.Mesh diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 0627304c6..923de05ee 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -21,6 +21,7 @@ from festim import ( reaction as _reaction, ) +from festim import source as _source from festim import ( species as _species, ) @@ -29,9 +30,8 @@ ) from festim.helpers import as_fenics_constant from festim.mesh import Mesh -from festim import source as _source -__all__ = ["HydrogenTransportProblem", "HTransportProblemDiscontinuous"] +__all__ = ["HTransportProblemDiscontinuous", "HydrogenTransportProblem"] class HydrogenTransportProblem(problem.ProblemBase): @@ -1063,7 +1063,7 @@ def create_source_values_fenics(self): # create value_fenics for all F.ParticleSource objects if isinstance(source, _source.ParticleSource): V = dolfinx.fem.functionspace(self.mesh.mesh, ("Lagrange", 1)) - + source.value.convert_input_value( mesh=self.mesh.mesh, function_space=V, From 1022bef64455c7a6fd471b820dadf92f3598f3c0 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 21 Feb 2025 15:10:45 -0500 Subject: [PATCH 32/40] refactor source: better docs, typehinting --- src/festim/source.py | 107 ++++++++++++++++++++++++++++++++----------- 1 file changed, 79 insertions(+), 28 deletions(-) diff --git a/src/festim/source.py b/src/festim/source.py index 393d5aede..525062c67 100644 --- a/src/festim/source.py +++ b/src/festim/source.py @@ -1,37 +1,36 @@ -import festim as F +import numpy as np +import ufl +from dolfinx import fem + +from festim.helpers import Value +from festim.species import Species +from festim.subdomain.volume_subdomain import VolumeSubdomain class SourceBase: """ - Source class + Source base class Args: - volume (festim.VolumeSubdomain1D): the volume subdomains where the source is applied - value (float, int, fem.Constant or callable): the value of the soure - species (festim.Species): the species to which the source is applied + value: the value of the source + volume: the volume subdomains where the source is applied Attributes: - volume (festim.VolumeSubdomain1D): the volume subdomains where the source is applied - value (float, int, fem.Constant or callable): the value of the soure - species (festim.Species): the species to which the source is applied - value_fenics (fem.Function or fem.Constant): the value of the source in - fenics format - source_expr (fem.Expression): the expression of the source term that is - used to update the value_fenics - time_dependent (bool): True if the value of the source is time dependent - temperature_dependent (bool): True if the value of the source is temperature - dependent - - Usage: - >>> from festim import Source - >>> Source(volume=my_vol, value=1, species="H") - >>> Source(volume=my_vol, value=lambda x: 1 + x[0], species="H") - >>> Source(volume=my_vol, value=lambda t: 1 + t, species="H") - >>> Source(volume=my_vol, value=lambda T: 1 + T, species="H") - >>> Source(volume=my_vol, value=lambda x, t: 1 + x[0] + t, species="H") + value: the value of the source + volume: the volume subdomains where the source is applied """ - def __init__(self, value, volume): + value: float | int | fem.Constant | np.ndarray | fem.Expression | \ + ufl.core.expr.Expr | fem.Function + volume: VolumeSubdomain + + def __init__( + self, + value: float | int | fem.Constant | np.ndarray | fem.Expression | + ufl.core.expr.Expr | fem.Function, + volume: VolumeSubdomain + ): + self.value = value self.volume = volume @@ -41,7 +40,7 @@ def value(self): @value.setter def value(self, value): - self._value = F.Value(value) + self._value = Value(value) @property def volume(self): @@ -50,13 +49,42 @@ def volume(self): @volume.setter def volume(self, value): # check that volume is festim.VolumeSubdomain - if not isinstance(value, F.VolumeSubdomain): + if not isinstance(value, VolumeSubdomain): raise TypeError("volume must be of type festim.VolumeSubdomain") self._volume = value class ParticleSource(SourceBase): - def __init__(self, value, volume, species): + """ + Particle source class + + Args: + value: the value of the source + volume: the volume subdomains where the source is applied + species: the species to which the source is applied + + Attributes: + value: the value of the source + volume: the volume subdomains where the source is applied + species: the species to which the source is applied + + Examples: + + .. highlight:: python + .. code-block:: python + + from festim import ParticleSource + + ParticleSource(volume=my_vol, value=1, species="H") + ParticleSource(volume=my_vol, value=lambda x: 1 + x[0], species="H") + ParticleSource(volume=my_vol, value=lambda t: 1 + t, species="H") + ParticleSource(volume=my_vol, value=lambda T: 1 + T, species="H") + ParticleSource(volume=my_vol, value=lambda x, t: 1 + x[0] + t, species="H") + """ + + species: Species + + def __init__(self, value, volume, species: Species): self.species = species super().__init__(value, volume) @@ -67,13 +95,36 @@ def species(self): @species.setter def species(self, value): # check that species is festim.Species or list of festim.Species - if not isinstance(value, F.Species): + if not isinstance(value, Species): raise TypeError("species must be of type festim.Species") self._species = value class HeatSource(SourceBase): + """ + Heat source class + + Args: + value: the value of the source + volume: the volume subdomains where the source is applied + + Attributes: + value: the value of the source + volume: the volume subdomains where the source is applied + + Examples: + + .. highlight:: python + .. code-block:: python + + from festim import HeatSource + + HeatSource(volume=my_vol, value=1) + HeatSource(volume=my_vol, value=lambda x: 1 + x[0]) + HeatSource(volume=my_vol, value=lambda t: 1 + t) + HeatSource(volume=my_vol, value=lambda x, t: 1 + x[0] + t) + """ def __init__(self, value, volume): super().__init__(value, volume) From d800685667ad20c14a5fadcc00142257a9317a83 Mon Sep 17 00:00:00 2001 From: jhdark Date: Mon, 24 Feb 2025 10:43:09 -0500 Subject: [PATCH 33/40] format ruff --- src/festim/helpers.py | 12 ++++++++++-- src/festim/source.py | 28 ++++++++++++++++++++-------- 2 files changed, 30 insertions(+), 10 deletions(-) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 99b6d8b2c..bb8c7e923 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -117,7 +117,15 @@ class Value: """ - input_value: float | int | fem.Constant | np.ndarray | fem.Expression | ufl.core.expr.Expr | fem.Function + input_value: ( + float + | int + | fem.Constant + | np.ndarray + | fem.Expression + | ufl.core.expr.Expr + | fem.Function + ) ufl_expression: ufl.core.expr.Expr fenics_interpolation_expression: fem.Expression @@ -202,7 +210,7 @@ def convert_input_value( function_space: the function space of the fenics object, optional t: the time, optional temperature: the temperature, optional - up_to_ufl_expr: if True, the value is only mapped to a function if the input + up_to_ufl_expr: if True, the value is only mapped to a function if the input is callable, not interpolated or converted to a function, optional """ if isinstance( diff --git a/src/festim/source.py b/src/festim/source.py index 525062c67..63842feca 100644 --- a/src/festim/source.py +++ b/src/festim/source.py @@ -20,17 +20,28 @@ class SourceBase: volume: the volume subdomains where the source is applied """ - value: float | int | fem.Constant | np.ndarray | fem.Expression | \ - ufl.core.expr.Expr | fem.Function + value: ( + float + | int + | fem.Constant + | np.ndarray + | fem.Expression + | ufl.core.expr.Expr + | fem.Function + ) volume: VolumeSubdomain def __init__( - self, - value: float | int | fem.Constant | np.ndarray | fem.Expression | - ufl.core.expr.Expr | fem.Function, - volume: VolumeSubdomain - ): - + self, + value: float + | int + | fem.Constant + | np.ndarray + | fem.Expression + | ufl.core.expr.Expr + | fem.Function, + volume: VolumeSubdomain, + ): self.value = value self.volume = volume @@ -125,6 +136,7 @@ class HeatSource(SourceBase): HeatSource(volume=my_vol, value=lambda t: 1 + t) HeatSource(volume=my_vol, value=lambda x, t: 1 + x[0] + t) """ + def __init__(self, value, volume): super().__init__(value, volume) From f0a640be649df6a20bb87b58d91581493786f65e Mon Sep 17 00:00:00 2001 From: jhdark Date: Mon, 24 Feb 2025 11:23:14 -0500 Subject: [PATCH 34/40] use functionspace to get mesh --- src/festim/heat_transfer_problem.py | 1 - src/festim/helpers.py | 8 ++++---- src/festim/hydrogen_transport_problem.py | 2 -- test/test_helpers.py | 26 ++++++------------------ test/test_source.py | 10 ++++++--- 5 files changed, 17 insertions(+), 30 deletions(-) diff --git a/src/festim/heat_transfer_problem.py b/src/festim/heat_transfer_problem.py index 19cd6205e..b02b5d4f3 100644 --- a/src/festim/heat_transfer_problem.py +++ b/src/festim/heat_transfer_problem.py @@ -135,7 +135,6 @@ def create_source_values_fenics(self): for source in self.sources: # create value_fenics for all source objects source.value.convert_input_value( - mesh=self.mesh.mesh, function_space=self.function_space, t=self.t, up_to_ufl_expr=True, diff --git a/src/festim/helpers.py b/src/festim/helpers.py index bb8c7e923..7f518bda6 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -196,7 +196,6 @@ def temperature_dependent(self) -> bool: def convert_input_value( self, - mesh: Optional[dolfinx.mesh.Mesh] = None, function_space: Optional[dolfinx.fem.function.FunctionSpace] = None, t: Optional[fem.Constant] = None, temperature: Optional[fem.Function | fem.Constant | ufl.core.expr.Expr] = None, @@ -206,7 +205,6 @@ def convert_input_value( on the type of the value provided Args: - mesh: the mesh of the domain, optional function_space: the function space of the fenics object, optional t: the time, optional temperature: the temperature, optional @@ -222,7 +220,9 @@ def convert_input_value( self.fenics_interpolation_expression = self.input_value elif isinstance(self.input_value, float | int): - self.fenics_object = as_fenics_constant(value=self.input_value, mesh=mesh) + self.fenics_object = as_fenics_constant( + value=self.input_value, mesh=function_space.mesh + ) elif callable(self.input_value): args = self.input_value.__code__.co_varnames @@ -235,7 +235,7 @@ def convert_input_value( ) self.fenics_object = as_fenics_constant( - value=self.input_value(t=float(t)), mesh=mesh + value=self.input_value(t=float(t)), mesh=function_space.mesh ) elif up_to_ufl_expr: diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 923de05ee..dd14a763b 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -618,7 +618,6 @@ def create_source_values_fenics(self): # create value_fenics for all F.ParticleSource objects if isinstance(source, _source.ParticleSource): source.value.convert_input_value( - mesh=self.mesh.mesh, function_space=self.function_space, t=self.t, temperature=self.temperature_fenics, @@ -1065,7 +1064,6 @@ def create_source_values_fenics(self): V = dolfinx.fem.functionspace(self.mesh.mesh, ("Lagrange", 1)) source.value.convert_input_value( - mesh=self.mesh.mesh, function_space=V, t=self.t, temperature=self.temperature_fenics, diff --git a/test/test_helpers.py b/test/test_helpers.py index ab6697f7a..fe0b99993 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -12,13 +12,7 @@ test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) x = ufl.SpatialCoordinate(test_mesh.mesh) -element_CG = basix.ufl.element( - basix.ElementFamily.P, - test_mesh.mesh.basix_cell(), - 1, - basix.LagrangeVariant.equispaced, -) -test_function_space = fem.functionspace(test_mesh.mesh, element_CG) +test_function_space = fem.functionspace(test_mesh.mesh, ("Lagrange", 1)) test_function = fem.Function(test_function_space) @@ -48,11 +42,9 @@ def test_temperature_type_and_processing(value): def test_value_convert_float_int_inputs(input_value, expected_output_type): """Test that float and value is correctly converted""" - mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10) - test_value = F.Value(input_value) - test_value.convert_input_value(mesh=mesh) + test_value.convert_input_value(function_space=test_function_space) assert isinstance(test_value.fenics_object, expected_output_type) @@ -82,7 +74,6 @@ def test_value_convert_up_to_ufl_inputs(input_value, expected_output_type): test_value = F.Value(input_value) test_value.convert_input_value( - mesh=my_mesh, function_space=V, t=my_t, temperature=my_T, @@ -113,19 +104,12 @@ def test_value_convert_callable_inputs(input_value, expected_output_type): my_t = fem.Constant(my_mesh, default_scalar_type(8)) my_T = fem.Constant(my_mesh, default_scalar_type(5)) - element_CG = basix.ufl.element( - basix.ElementFamily.P, - my_mesh.basix_cell(), - 1, - basix.LagrangeVariant.equispaced, - ) - my_function_space = fem.functionspace(my_mesh, element_CG) + my_function_space = fem.functionspace(my_mesh, ("Lagrange", 1)) test_value = F.Value(input_value) test_value.convert_input_value( function_space=my_function_space, - mesh=my_mesh, t=my_t, temperature=my_T, ) @@ -248,7 +232,9 @@ def my_value(t): ValueError, match="self.value should return a float or an int, not Date: Tue, 25 Feb 2025 09:51:29 -0500 Subject: [PATCH 35/40] Update src/festim/helpers.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Rémi Delaporte-Mathurin <40028739+RemDelaporteMathurin@users.noreply.github.com> --- src/festim/helpers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 7f518bda6..16b114db7 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -169,7 +169,7 @@ def input_value(self, value): ) @property - def time_dependent(self) -> bool: + def explicit_time_dependent(self) -> bool: """Returns true if the value given is time dependent""" if self.input_value is None: return False From dc76b643833e7398d6d82c3397c367086d45891f Mon Sep 17 00:00:00 2001 From: jhdark Date: Tue, 25 Feb 2025 10:23:52 -0500 Subject: [PATCH 36/40] create fucntion spaces first, use function space from subdomain --- src/festim/hydrogen_transport_problem.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index dd14a763b..b58492d4a 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -934,12 +934,15 @@ def initialise(self): self.create_implicit_species_value_fenics() + for subdomain in self.volume_subdomains: + self.define_function_spaces(subdomain) + self.define_temperature() self.create_source_values_fenics() self.create_flux_values_fenics() self.create_initial_conditions() + for subdomain in self.volume_subdomains: - self.define_function_spaces(subdomain) self.create_subdomain_formulation(subdomain) subdomain.u.name = f"u_{subdomain.id}" @@ -1061,14 +1064,15 @@ def create_source_values_fenics(self): for source in self.sources: # create value_fenics for all F.ParticleSource objects if isinstance(source, _source.ParticleSource): - V = dolfinx.fem.functionspace(self.mesh.mesh, ("Lagrange", 1)) + for subdomain in source.species.subdomains: + V = source.species.subdomain_to_function_space[subdomain] - source.value.convert_input_value( - function_space=V, - t=self.t, - temperature=self.temperature_fenics, - up_to_ufl_expr=True, - ) + source.value.convert_input_value( + function_space=V, + t=self.t, + temperature=self.temperature_fenics, + up_to_ufl_expr=True, + ) def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain): """ From 422e4f3b0a0a837baf6e223b9c6672c64bf4670e Mon Sep 17 00:00:00 2001 From: jhdark Date: Tue, 25 Feb 2025 10:26:46 -0500 Subject: [PATCH 37/40] time dependent propety changed to explicitly time dependent --- src/festim/helpers.py | 5 +++++ src/festim/problem.py | 2 +- test/test_helpers.py | 2 +- test/test_source.py | 4 ++-- 4 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 16b114db7..d0e4d91c5 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -114,6 +114,9 @@ class Value: fenics_interpolation_expression : The expression of the user input that is used to update the `fenics_object` fenics_object : The value of the user input in fenics format + explicit_time_dependent : True if the user input value is explicitly time + dependent + temperature_dependent : True if the user input value is temperature dependent """ @@ -130,6 +133,8 @@ class Value: ufl_expression: ufl.core.expr.Expr fenics_interpolation_expression: fem.Expression fenics_object: fem.Function | fem.Constant | ufl.core.expr.Expr + explicit_time_dependent: bool + temperature_dependent: bool def __init__(self, input_value): self.input_value = input_value diff --git a/src/festim/problem.py b/src/festim/problem.py index e40cc5638..5463d88fd 100644 --- a/src/festim/problem.py +++ b/src/festim/problem.py @@ -193,5 +193,5 @@ def update_time_dependent_values(self): bc.update(t=t) for source in self.sources: - if source.value.time_dependent: + if source.value.explicit_time_dependent: source.value.update(t=t) diff --git a/test/test_helpers.py b/test/test_helpers.py index fe0b99993..b282dac8d 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -152,7 +152,7 @@ def test_time_dependent_values(input_value, expected_output): test_value = F.Value(input_value) - assert test_value.time_dependent == expected_output + assert test_value.explicit_time_dependent == expected_output @pytest.mark.parametrize( diff --git a/test/test_source.py b/test/test_source.py index 32576cccb..676d346a9 100644 --- a/test/test_source.py +++ b/test/test_source.py @@ -82,13 +82,13 @@ def test_create_fenics_object(value, expected_type): (lambda x, t: ufl.conditional(ufl.lt(t, 1.0), 100.0 + x[0], 0.0), True), ], ) -def test_source_time_dependent_attribute(input, expected_value): +def test_source_explicit_time_dependent_attribute(input, expected_value): """Test that the time_dependent attribute is correctly set""" volume = F.VolumeSubdomain1D(1, borders=[0, 1], material=dummy_mat) species = F.Species("test") my_source = F.ParticleSource(input, volume, species) - assert my_source.value.time_dependent is expected_value + assert my_source.value.explicit_time_dependent is expected_value @pytest.mark.parametrize( From d004c62025f3e89f6507bf3d3646f4dd88fc74ea Mon Sep 17 00:00:00 2001 From: jhdark Date: Tue, 25 Feb 2025 10:37:03 -0500 Subject: [PATCH 38/40] rename function processing sources --- src/festim/hydrogen_transport_problem.py | 8 ++++---- src/festim/problem_change_of_var.py | 2 +- test/test_h_transport_problem.py | 8 ++++---- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index b58492d4a..134b5d1d3 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -279,7 +279,7 @@ def initialise(self): self.define_temperature() self.define_boundary_conditions() - self.create_source_values_fenics() + self.convert_source_input_values_to_fenics_objects() self.create_flux_values_fenics() self.create_initial_conditions() self.create_formulation() @@ -612,7 +612,7 @@ def create_dirichletbc_form(self, bc): return form - def create_source_values_fenics(self): + def convert_source_input_values_to_fenics_objects(self): """For each source create the value_fenics""" for source in self.sources: # create value_fenics for all F.ParticleSource objects @@ -938,7 +938,7 @@ def initialise(self): self.define_function_spaces(subdomain) self.define_temperature() - self.create_source_values_fenics() + self.convert_source_input_values_to_fenics_objects() self.create_flux_values_fenics() self.create_initial_conditions() @@ -1059,7 +1059,7 @@ def define_function_spaces(self, subdomain: _subdomain.VolumeSubdomain): name = f"{species.name}_{subdomain.id}" species.subdomain_to_post_processing_solution[subdomain].name = name - def create_source_values_fenics(self): + def convert_source_input_values_to_fenics_objects(self): """For each source create the value_fenics""" for source in self.sources: # create value_fenics for all F.ParticleSource objects diff --git a/src/festim/problem_change_of_var.py b/src/festim/problem_change_of_var.py index 4f39285ec..273691cf3 100644 --- a/src/festim/problem_change_of_var.py +++ b/src/festim/problem_change_of_var.py @@ -30,7 +30,7 @@ def initialise(self): self.define_temperature() self.define_boundary_conditions() - self.create_source_values_fenics() + self.convert_source_input_values_to_fenics_objects() self.create_flux_values_fenics() self.create_initial_conditions() self.create_formulation() diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py index 200420221..4f3e2702a 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -673,7 +673,7 @@ def test_update_time_dependent_values_source(source_value, expected_values): my_model.define_meshtags_and_measures() my_model.assign_functions_to_species() my_model.define_temperature() - my_model.create_source_values_fenics() + my_model.convert_source_input_values_to_fenics_objects() for i in range(3): # RUN @@ -727,7 +727,7 @@ def test_update_sources_with_time_dependent_temperature( my_model.define_function_spaces() my_model.assign_functions_to_species() my_model.define_meshtags_and_measures() - my_model.create_source_values_fenics() + my_model.convert_source_input_values_to_fenics_objects() for i in range(3): # RUN @@ -740,7 +740,7 @@ def test_update_sources_with_time_dependent_temperature( assert np.isclose(computed_value, expected_values[i]) -def test_create_source_values_fenics_multispecies(): +def test_convert_source_input_values_to_fenics_objects_multispecies(): """Test that the define_sources method correctly sets the value_fenics attribute in a multispecies case""" # BUILD @@ -764,7 +764,7 @@ def test_create_source_values_fenics_multispecies(): my_model.define_temperature() # RUN - my_model.create_source_values_fenics() + my_model.convert_source_input_values_to_fenics_objects() # TEST assert np.isclose(float(my_model.sources[0].value.fenics_object), 5) From 9c9ef4db2547a7b7b692acc225ec7fcad53e7068 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 28 Feb 2025 11:24:28 -0500 Subject: [PATCH 39/40] support nightly --- src/festim/helpers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 29a32b8bd..956d42ab0 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -93,7 +93,7 @@ def as_fenics_interp_expr_and_function( fenics_interpolation_expression = fem.Expression( mapped_function, - function_space.element.interpolation_points(), + get_interpolation_points(function_space.element), ) fenics_object = fem.Function(function_space) From 4e04da8ea6c6e2a66760ee9ae11858c88914ea96 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 28 Feb 2025 11:32:25 -0500 Subject: [PATCH 40/40] new method for interpolation points --- test/test_helpers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_helpers.py b/test/test_helpers.py index b282dac8d..41bcb6e90 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -207,7 +207,7 @@ def test_input_values_of_expressions_are_accepted(): test_expression = fem.Expression( mapped_func, - test_function_space.element.interpolation_points(), + F.get_interpolation_points(test_function_space.element), ) test_value = F.Value(input_value=test_expression)