diff --git a/festim/temperature/temperature_solver.py b/festim/temperature/temperature_solver.py index 92297c97d..b3111b398 100644 --- a/festim/temperature/temperature_solver.py +++ b/festim/temperature/temperature_solver.py @@ -8,8 +8,8 @@ class HeatTransferProblem(festim.Temperature): Args: transient (bool, optional): If True, a transient simulation will be run. Defaults to True. - initial_value (sp.Add, float, optional): The initial value. - Only needed if transient is True. Defaults to 0. + initial_condition (festim.InitialCondition, optional): The initial condition. + Only needed if transient is True. absolute_tolerance (float, optional): the absolute tolerance of the newton solver. Defaults to 1e-03 relative_tolerance (float, optional): the relative tolerance of the newton @@ -25,7 +25,7 @@ class HeatTransferProblem(festim.Temperature): Attributes: F (fenics.Form): the variational form of the heat transfer problem v_T (fenics.TestFunction): the test function - initial_value (sp.Add, int, float): the initial value + initial_condition (festim.InitialCondition): the initial condition sub_expressions (list): contains time dependent fenics.Expression to be updated sources (list): contains festim.Source objects for volumetric heat @@ -36,7 +36,7 @@ class HeatTransferProblem(festim.Temperature): def __init__( self, transient=True, - initial_value=0.0, + initial_condition=None, absolute_tolerance=1e-3, relative_tolerance=1e-10, maximum_iterations=30, @@ -44,7 +44,7 @@ def __init__( ) -> None: super().__init__() self.transient = transient - self.initial_value = initial_value + self.initial_condition = initial_condition self.absolute_tolerance = absolute_tolerance self.relative_tolerance = relative_tolerance self.maximum_iterations = maximum_iterations @@ -74,10 +74,19 @@ def create_functions(self, materials, mesh, dt=None): self.T_n = f.Function(V, name="T_n") self.v_T = f.TestFunction(V) - if self.transient: - ccode_T_ini = sp.printing.ccode(self.initial_value) - self.initial_value = f.Expression(ccode_T_ini, degree=2, t=0) - self.T_n.assign(f.interpolate(self.initial_value, V)) + if self.transient and self.initial_condition: + if isinstance(self.initial_condition.value, str): + if self.initial_condition.value.endswith(".xdmf"): + with f.XDMFFile(self.initial_condition.value) as file: + file.read_checkpoint( + self.T_n, + self.initial_condition.label, + self.initial_condition.time_step, + ) + else: + ccode_T_ini = sp.printing.ccode(self.initial_condition.value) + self.initial_condition.value = f.Expression(ccode_T_ini, degree=2, t=0) + self.T_n.assign(f.interpolate(self.initial_condition.value, V)) self.define_variational_problem(materials, mesh, dt) self.create_dirichlet_bcs(mesh.surface_markers) diff --git a/test/system/test_system.py b/test/system/test_system.py index 8bdfe6e0a..f30d5fceb 100644 --- a/test/system/test_system.py +++ b/test/system/test_system.py @@ -107,7 +107,9 @@ def test_run_temperature_transient(tmpdir): festim.DirichletBC(surfaces=[1, 2], value=u, field="T"), ] - my_temp = festim.HeatTransferProblem(transient=True, initial_value=u) + my_temp = festim.HeatTransferProblem( + transient=True, initial_condition=festim.InitialCondition(field="T", value=u) + ) my_sources = [ festim.Source( diff --git a/test/unit/test_temperature.py b/test/unit/test_temperature.py index 1891972cb..93140ff1e 100644 --- a/test/unit/test_temperature.py +++ b/test/unit/test_temperature.py @@ -69,7 +69,9 @@ def thermal_cond(a): bc1 = festim.DirichletBC(surfaces=[1], value=u, field="T") bc2 = festim.FluxBC(surfaces=[2], value=2, field="T") - my_temp = festim.HeatTransferProblem(transient=True, initial_value=0) + my_temp = festim.HeatTransferProblem( + transient=True, initial_condition=festim.InitialCondition(field="T", value=0) + ) my_temp.boundary_conditions = [bc1, bc2] my_temp.sources = [festim.Source(-4, volume=[1], field="T")] my_temp.create_functions(my_mats, my_mesh, dt=dt) @@ -93,6 +95,41 @@ def thermal_cond(a): assert expected_form.equals(F) +def test_heat_transfer_create_functions_transient(tmpdir): + """Test for the HeatTransferProblem.create_functions(). + Creates a function, writes it to an XDMF file, then a HeatTransferProblem + class is created from this file and the error norm between the written and + read fuctions is computed to ensure they are the same. + """ + # create function to be comapared + mesh = fenics.UnitIntervalMesh(10) + V = fenics.FunctionSpace(mesh, "CG", 1) + expr = fenics.Expression("1 + x[0]", degree=2) + T = fenics.interpolate(expr, V) + # write function to temporary file + T_file = tmpdir.join("T.xdmf") + fenics.XDMFFile(str(Path(T_file))).write_checkpoint( + T, "T", 0, fenics.XDMFFile.Encoding.HDF5, append=False + ) + # HeatTransferProblem needs a festim mesh and a festim material + my_mesh = festim.Mesh(mesh) + my_mesh.dx = fenics.dx + my_mesh.ds = fenics.ds + my_mats = festim.Materials( + [festim.Material(id=1, D_0=1, E_D=0, thermal_cond=1, heat_capacity=1, rho=1)] + ) + my_temp = festim.HeatTransferProblem( + transient=True, + initial_condition=festim.InitialCondition( + field="T", value=str(Path(T_file)), label="T", time_step=0 + ), + ) + my_temp.create_functions(my_mats, my_mesh, dt=festim.Stepsize(initial_value=2)) + # evaluate error between original and read function + error_L2 = fenics.errornorm(T, my_temp.T_n, "L2") + assert error_L2 < 1e-9 + + def test_temperature_from_xdmf_create_functions(tmpdir): """Test for the TemperatureFromXDMF.create_functions(). Creates a function, writes it to an XDMF file, then a TemperatureFromXDMF