diff --git a/src/festim/__init__.py b/src/festim/__init__.py index 8fa79b968..ac904f92a 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, get_interpolation_points +from .helpers import ( + as_fenics_constant, + as_mapped_function, + as_fenics_interp_expr_and_function, + Value, + get_interpolation_points, +) from .hydrogen_transport_problem import ( HTransportProblemDiscontinuous, HydrogenTransportProblem, diff --git a/src/festim/heat_transfer_problem.py b/src/festim/heat_transfer_problem.py index 9c49fce6c..b02b5d4f3 100644 --- a/src/festim/heat_transfer_problem.py +++ b/src/festim/heat_transfer_problem.py @@ -134,9 +134,10 @@ 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, + source.value.convert_input_value( + function_space=self.function_space, t=self.t, + up_to_ufl_expr=True, ) def create_flux_values_fenics(self): @@ -205,7 +206,9 @@ 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 diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 7afcfa7a0..956d42ab0 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -1,5 +1,11 @@ +from collections.abc import Callable +from typing import Optional + import dolfinx +import numpy as np +import ufl from dolfinx import fem +from packaging import version def as_fenics_constant( @@ -17,8 +23,8 @@ def as_fenics_constant( 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)) + if isinstance(value, float | int): + return fem.Constant(mesh, dolfinx.default_scalar_type(float(value))) elif isinstance(value, fem.Constant): return value else: @@ -27,7 +33,251 @@ def as_fenics_constant( ) -from packaging import version +def as_mapped_function( + value: Callable, + 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 + function_space: the function space of the domain, optional + t: the time, optional + temperature: the temperature, optional + + Returns: + The mapped function + """ + + # Extract the input variable names in the callable function `value` + arguments = value.__code__.co_varnames + + kwargs = {} + if "t" in arguments: + kwargs["t"] = t + if "x" in arguments: + x = ufl.SpatialCoordinate(function_space.mesh) + kwargs["x"] = x + if "T" in arguments: + kwargs["T"] = temperature + + return value(**kwargs) + + +def as_fenics_interp_expr_and_function( + value: Callable, + function_space: dolfinx.fem.function.FunctionSpace, + 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 + expression and function objects + + Args: + value: the callable to convert + function_space: The function space to interpolate function over + t: the time, optional + temperature: the temperature, optional + + Returns: + fenics interpolation expression, fenics function + """ + + mapped_function = as_mapped_function( + value=value, function_space=function_space, t=t, temperature=temperature + ) + + fenics_interpolation_expression = fem.Expression( + mapped_function, + get_interpolation_points(function_space.element), + ) + + 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: + 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_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 + + """ + + 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 + 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 + + self.ufl_expression = None + 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 + + @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.Expression + | 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, np.ndarray, fem.Expression," + f" ufl.core.expr.Expr, fem.Function, or callable not {value}" + ) + + @property + def explicit_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 | ufl.core.expr.Expr): + 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 | ufl.core.expr.Expr): + return False + if callable(self.input_value): + arguments = self.input_value.__code__.co_varnames + return "T" in arguments + else: + return False + + def convert_input_value( + self, + 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: + 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 + ): + 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, float | int): + 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 + # 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.fenics_object = as_fenics_constant( + value=self.input_value(t=float(t)), mesh=function_space.mesh + ) + + elif up_to_ufl_expr: + self.fenics_object = as_mapped_function( + value=self.input_value, + function_space=function_space, + t=t, + temperature=temperature, + ) + + else: + self.fenics_interpolation_expression, self.fenics_object = ( + as_fenics_interp_expr_and_function( + value=self.input_value, + function_space=function_space, + t=t, + temperature=temperature, + ) + ) + + def update(self, t: float): + """Updates the value + + Args: + t: 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) + # Check the version of dolfinx dolfinx_version = dolfinx.__version__ diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index ef7b45f51..1898392f0 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, ) @@ -30,7 +31,7 @@ from festim.helpers import as_fenics_constant, get_interpolation_points from festim.mesh import Mesh -__all__ = ["HydrogenTransportProblem", "HTransportProblemDiscontinuous"] +__all__ = ["HTransportProblemDiscontinuous", "HydrogenTransportProblem"] class HydrogenTransportProblem(problem.ProblemBase): @@ -278,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() @@ -611,15 +612,16 @@ 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 - if isinstance(source, festim.source.ParticleSource): - source.create_value_fenics( - mesh=self.mesh.mesh, - temperature=self.temperature_fenics, + if isinstance(source, _source.ParticleSource): + source.value.convert_input_value( + function_space=self.function_space, t=self.t, + temperature=self.temperature_fenics, + up_to_ufl_expr=True, ) def create_flux_values_fenics(self): @@ -713,7 +715,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) ) @@ -783,8 +785,8 @@ def update_time_dependent_values(self): bc.update(t=t) for source in self.sources: - if source.temperature_dependent: - source.update(t=t) + if source.value.temperature_dependent: + source.value.update(t=t) def post_processing(self): """Post processes the model""" @@ -932,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.convert_source_input_values_to_fenics_objects() 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}" @@ -1054,6 +1059,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 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 + if isinstance(source, _source.ParticleSource): + 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, + ) + def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain): """ Creates the variational formulation for each subdomain and stores it in ``subdomain.F`` @@ -1123,7 +1143,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 diff --git a/src/festim/problem.py b/src/festim/problem.py index d1066ef79..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.time_dependent: - source.update(t=t) + if source.value.explicit_time_dependent: + source.value.update(t=t) diff --git a/src/festim/problem_change_of_var.py b/src/festim/problem_change_of_var.py index 76cc905a3..643164af2 100644 --- a/src/festim/problem_change_of_var.py +++ b/src/festim/problem_change_of_var.py @@ -31,7 +31,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() @@ -105,7 +105,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) ) diff --git a/src/festim/source.py b/src/festim/source.py index a792dda49..63842feca 100644 --- a/src/festim/source.py +++ b/src/festim/source.py @@ -1,47 +1,57 @@ -import dolfinx import numpy as np import ufl from dolfinx import fem -import festim as F +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 - self.value_fenics = None - self.source_expr = None + @property + def value(self): + return self._value + + @value.setter + def value(self, value): + self._value = Value(value) @property def volume(self): @@ -50,53 +60,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 - @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): + """ + 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 -class ParticleSource(SourceBase): - def __init__(self, value, volume, species): + 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) @@ -107,112 +106,39 @@ 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 - @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): + """ + 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) - 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) + 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 b6639d8a4..4f3e2702a 100644 --- a/test/test_h_transport_problem.py +++ b/test/test_h_transport_problem.py @@ -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.temperature 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) @@ -669,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 @@ -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_object, fem.Constant): + computed_value = float(my_model.sources[0].value.fenics_object) assert np.isclose(computed_value, expected_values[i]) @@ -723,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 @@ -731,12 +735,12 @@ 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_object, fem.Constant): + computed_value = float(my_model.sources[0].value.fenics_object) 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 @@ -760,11 +764,11 @@ 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(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_object), 5) + assert np.isclose(float(my_model.sources[1].value.fenics_object), 11) # TODO replace this by a proper MMS test @@ -1229,5 +1233,5 @@ 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), 5) + assert np.isclose(float(my_model.boundary_conditions[1].value_fenics), 11) diff --git a/test/test_heat_transfer_problem.py b/test/test_heat_transfer_problem.py index 224802b7e..885337f69 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], diff --git a/test/test_helpers.py b/test/test_helpers.py index 9220b1451..41bcb6e90 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,9 @@ test_mesh = F.Mesh1D(vertices=np.array([0.0, 1.0, 2.0, 3.0, 4.0])) x = ufl.SpatialCoordinate(test_mesh.mesh) +test_function_space = fem.functionspace(test_mesh.mesh, ("Lagrange", 1)) +test_function = fem.Function(test_function_space) + @pytest.mark.parametrize( "value", @@ -22,8 +29,231 @@ 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""" + + test_value = F.Value(input_value) + + test_value.convert_input_value(function_space=test_function_space) + + 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) + 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)) + + test_value = F.Value(input_value) + + test_value.convert_input_value( + function_space=V, + 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)) + + 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, + 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, np.ndarray, fem.Expression, " + "ufl.core.expr.Expr, fem.Function, or callable not coucou" + ), + ): + 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.explicit_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_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, + F.get_interpolation_points(test_function_space.element), + ) + 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""" + + 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