Merge pull request #662 from KulaginVladimir/main

InitialCondition for HeatTransferProblem

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,15 +36,15 @@ 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,
         linear_solver=None,
     ) -> 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)

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(

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