Verify second order accurate enthalpy-porosity method via MMS (#43)

* Test second order accuracy (in time and space) for enthalpy-porosity method. This shows spatial order of 1.7 and temporal order of 1.8. Need more computing time to do better.

* Correct initialization of solution for subsequent time step sizes during MMS temporal accuracy verification

* Test MMS plotting with pytest

* Create subdirectories with h and Delta_t values when plotting MMS solutions

* Have unsteady_model plot initial values before running time loop

fempy/mms.py

@@ -96,6 +96,18 @@ def verify_spatial_order_of_accuracy(
 
         model = MMSVerificationModel(meshsize = meshsize)
 
+        model.output_directory_path = model.output_directory_path.joinpath(
+            "mms_space")
+
+        if timestep_size is not None:
+
+            model.output_directory_path = \
+                model.output_directory_path.joinpath(
+                "Deltat" + str(timestep_size))
+
+        model.output_directory_path = model.output_directory_path.joinpath(
+            "m" + str(meshsize) + "/")
+
         model.assign_parameters(parameters)
 
         if hasattr(model, "time"):
@@ -104,7 +116,7 @@ def verify_spatial_order_of_accuracy(
 
             model.timestep_size.assign(timestep_size)
 
-            model.run(endtime = endtime)
+            model.run(endtime = endtime, plot = plot_solution)
 
         else:
 
@@ -156,11 +168,6 @@ def verify_spatial_order_of_accuracy(
 
         plt.close()
 
-
-    if plot_solution:
-
-        model.plot()
-
     assert(abs(order - expected_order) < tolerance)
 
 
@@ -182,6 +189,8 @@ def verify_temporal_order_of_accuracy(
 
     model = MMSVerificationModel(meshsize = meshsize)
 
+    basepath = model.output_directory_path
+
     model.assign_parameters(parameters)
 
     u_m = model.initial_values
@@ -198,6 +207,14 @@ def verify_temporal_order_of_accuracy(
 
         model.timestep_size.assign(timestep_size)
 
+        model.output_directory_path = basepath.joinpath("mms_time")
+
+        model.output_directory_path = model.output_directory_path.joinpath(
+            "m" + str(meshsize))
+
+        model.output_directory_path = model.output_directory_path.joinpath(
+            "Deltat" + str(timestep_size))
+
         model.time.assign(starttime)
 
         if (type(model.initial_values) == type((0,)) or
@@ -210,7 +227,9 @@ def verify_temporal_order_of_accuracy(
 
             model.initial_values.assign(initial_values[0])
 
-        model.run(endtime = endtime)
+        model.solution.assign(model.initial_values[0])
+
+        model.run(endtime = endtime, plot = plot_solution)
 
         table.append({
             "Delta_t": timestep_size,
@@ -257,11 +276,7 @@ def verify_temporal_order_of_accuracy(
         plt.savefig(str(filepath))
 
         plt.close()
-
-    if plot_solution:
-
-        model.plot()
-
+
     assert(abs(order - expected_order) < tolerance)
 
 

fempy/models/enthalpy_porosity.py

@@ -233,6 +233,11 @@ def plot(self):
 
 class ModelWithBDF2(Model):
 
+    def init_initial_values(self):
+
+        self.initial_values = [fe.Function(self.function_space)
+            for i in range(2)]
+
     def init_time_discrete_terms(self):
 
         Delta_t = self.timestep_size

fempy/unsteady_model.py

@@ -53,6 +53,10 @@ def push_back_initial_values(self):
 
     def run(self, endtime, plot = False):
 
+        if plot:
+
+            self.plot()
+
         while self.time.__float__() < (endtime - self.time_tolerance):
 
             self.time.assign(self.time + self.timestep_size)

tests/test__enthalpy_porosity.py

@@ -8,14 +8,21 @@ def test__melting_octadecane_benchmark__regression():
 
     endtime, expected_liquid_area, tolerance = 30., 0.24, 0.01
 
-    model = fempy.benchmarks.melting_octadecane.Model(meshsize = 32)
+    nx = 32
 
-    model.timestep_size.assign(10.)
+    Delta_t = 10.
 
-    model.smoothing.assign(1./256.)
+    model = fempy.benchmarks.melting_octadecane.Model(meshsize = nx)
+
+    model.timestep_size.assign(Delta_t)
+
+    s = 1./256.
+
+    model.smoothing.assign(s)
 
     model.output_directory_path = model.output_directory_path.joinpath(
-        "melting_octadecane/")
+        "melting_octadecane/" + "nx" + str(nx) + "_Deltat" + str(Delta_t) 
+        + "_s" + str(s) + "_tf" + str(endtime) + "/")
 
     model.run(endtime = endtime, plot = False)
 
@@ -76,7 +83,7 @@ def test__melting_octadecane_benchmark_with_darcy_resistance__regression():
     assert(abs(max_phil - 1.) < tolerance)
 
 
-class VerifiableModel(fempy.models.enthalpy_porosity.Model):
+class VerifiableModel(fempy.models.enthalpy_porosity.ModelWithBDF2):
 
     def __init__(self, meshsize):
 
@@ -116,11 +123,11 @@ def strong_form_residual(self, solution):
 
         mu = mu_s + (mu_l - mu_s)*phil
 
-        r_p = div(u) + gamma*p
+        r_p = div(u)
 
         r_u = diff(u, t) + grad(u)*u + grad(p) - 2.*div(mu*sym(grad(u))) + b
 
-        r_T = diff(T, t) + 1./Ste*diff(phil, t) + div(T*u) - 1./Pr*div(grad(T))
+        r_T = diff(T + 1./Ste*phil, t) + div(T*u) - 1./Pr*div(grad(T))
 
         return r_p, r_u, r_T
 
@@ -139,9 +146,9 @@ def init_manufactured_solution(self):
         u = exp(t)*sin(2.*pi*x[0])*sin(pi*x[1])*ihat + \
             exp(t)*sin(pi*x[0])*sin(2.*pi*x[1])*jhat
 
-        p = -0.5*(u[0]**2 + u[1]**2)
+        p = -sin(pi*x[0])*sin(2.*pi*x[1])
 
-        T = 0.5*sin(2.*pi*x[0])*sin(pi*x[1])*(1. - 2*exp(-3.*t**2))
+        T = 0.5*sin(2.*pi*x[0])*sin(pi*x[1])*(1. - exp(-t**2))
 
         self.manufactured_solution = p, u, T
 
@@ -161,15 +168,15 @@ def solve(self):
         print("Solved at time t = " + str(self.time.__float__()))
 
 
-def fails__test__verify_spatial_convergence_order_via_mms(
+def test__verify_spatial_convergence_order_via_mms(
         parameters = {
             "grashof_number": 2.,
             "prandtl_number": 5.,
             "stefan_number": 0.2,
             "smoothing": 1./16.},
-        mesh_sizes = (4, 8, 16),
-        timestep_size = 1./64.,
-        tolerance = 0.2):
+        mesh_sizes = (8, 16, 32),
+        timestep_size = 1./256.,
+        tolerance = 0.4):
 
     fempy.mms.verify_spatial_order_of_accuracy(
         Model = VerifiableModel,
@@ -178,5 +185,27 @@ def fails__test__verify_spatial_convergence_order_via_mms(
         mesh_sizes = mesh_sizes,
         tolerance = tolerance,
         timestep_size = timestep_size,
-        endtime = 1.)
-
+        endtime = 0.5,
+        plot_solution = False)
+
+
+def test__verify_temporal_convergence_order_via_mms(
+        parameters = {
+            "grashof_number": 2.,
+            "prandtl_number": 5.,
+            "stefan_number": 0.2,
+            "smoothing": 1./16.},
+        meshsize = 32,
+        timestep_sizes = (1./4., 1./8., 1./16., 1./32.),
+        tolerance = 0.2):
+
+    fempy.mms.verify_temporal_order_of_accuracy(
+        Model = VerifiableModel,
+        parameters = parameters,
+        expected_order = 2,
+        meshsize = meshsize,
+        tolerance = tolerance,
+        timestep_sizes = timestep_sizes,
+        endtime = 0.5,
+        plot_solution = False)
+

tests/test__heat.py

@@ -128,7 +128,8 @@ def test__verify_temporal_convergence_order_via_mms__bdf2(
         meshsize = meshsize,
         endtime = 1.,
         timestep_sizes = timestep_sizes,
-        tolerance = tolerance)
+        tolerance = tolerance,
+        plot_solution = True)
 
 
 class ModelWithWave(Model):

tests/test__laplace.py

@@ -43,5 +43,6 @@ def test__verify_convergence_order_via_mms(
         Model = OneDimMMSModel,
         expected_order = 2,
         mesh_sizes = mesh_sizes,
-        tolerance = tolerance)
+        tolerance = tolerance,
+        plot_solution = True)