Simplify the Simulation interface, improve MMS verification, and clea…

…n up validation tests (#76)

* Increment version number

Simulation class:
* Simplify interface for parameters / keyword arguments, including solver parameters. Before, there was a pattern where many parameters were set with init while others were assigned after instantiation. Parameters can still be assigned after instantiation, but the default pattern now is to include all parameters also as arguments to init. This cleans up a lot of the test/user code.

MMS:
* Report mesh cell count and function space dof count
* Simplify interface for passing simulation parameters

Convection-coupled phase-change:
* Make interfaces more flexible, allow for different element degrees, choose norms in MMS verification, use correct norms for velocity and temperature throughout verification tests
* MMS: Shift pressure solution to have zero gradient at boundaries; doesn't actually change the pressure error; but this will be useful when trying different boundary conditions.
* Remove option for top wall heat flux which is no longer relevant
* Use temperature function_space for the postprocessed porosity

sapphire/benchmarks/freeze_water_in_cavity.py

@@ -44,8 +44,8 @@ def heat_driven_cavity_variational_form_residual(sim, solution):
     mass = sapphire.simulations.convection_coupled_phasechange.\
         mass(sim, solution)
 
-    stabilization = sapphire.simulations.convection_coupled_phasechange.\
-        stabilization(sim, solution)
+    pressure_penalty = sapphire.simulations.convection_coupled_phasechange.\
+        pressure_penalty(sim, solution)
 
     p, u, T = fe.split(solution)
 
@@ -63,7 +63,7 @@ def heat_driven_cavity_variational_form_residual(sim, solution):
 
     energy = psi_T*dot(u, grad(T)) + dot(grad(psi_T), 1./Pr*grad(T))
 
-    return mass + momentum + energy + stabilization
+    return mass + momentum + energy + pressure_penalty
 
 
 def dirichlet_boundary_conditions(sim):
@@ -114,6 +114,7 @@ def initial_values(sim):
         problem = problem,
         solver_parameters = {
                 "snes_type": "newtonls",
+                "snes_max_it": sim.solver_parameters["snes_max_it"],
                 "snes_monitor": None,
                 "ksp_type": "preonly", 
                 "pc_type": "lu", 
@@ -150,40 +151,46 @@ def variational_form_residual(sim, solution):
                     solution = solution,
                     buoyancy = water_buoyancy),
             sapphire.simulations.convection_coupled_phasechange.energy,
-            sapphire.simulations.convection_coupled_phasechange.stabilization)])\
+            sapphire.simulations.convection_coupled_phasechange.pressure_penalty)])\
         *fe.dx(degree = sim.quadrature_degree)
 
 
 class Simulation(sapphire.simulations.convection_coupled_phasechange.Simulation):
 
-    def __init__(self, *args, meshsize, **kwargs):
-
-        self.reference_temperature_range__degC = fe.Constant(10.)
+    def __init__(self, *args, 
+            meshsize,
+            reference_temperature_range__degC = 10.,
+            cold_wall_temperature_before_freezing = 0.,
+            cold_wall_temperature_during_freezing = -1.,
+            stefan_number = 0.125,
+            liquidus_temperature = 0.,
+            density_solid_to_liquid_ratio = 916.70/999.84,
+            heat_capacity_solid_to_liquid_ratio = 0.500,
+            thermal_conductivity_solid_to_liquid_ratio = 2.14/0.561,
+            **kwargs):
+
+        self.reference_temperature_range__degC = fe.Constant(
+            reference_temperature_range__degC)
 
         self.hot_wall_temperature = fe.Constant(1.)
 
-        self.cold_wall_temperature = fe.Constant(0.)
+        self.cold_wall_temperature = fe.Constant(cold_wall_temperature_before_freezing)
 
         super().__init__(
             *args,
             mesh = fe.UnitSquareMesh(meshsize, meshsize),
             variational_form_residual = variational_form_residual,
             initial_values = initial_values,
             dirichlet_boundary_conditions = dirichlet_boundary_conditions,
+            stefan_number = stefan_number,
+            liquidus_temperature = liquidus_temperature,
+            density_solid_to_liquid_ratio = density_solid_to_liquid_ratio,
+            heat_capacity_solid_to_liquid_ratio = \
+                heat_capacity_solid_to_liquid_ratio,
+            thermal_conductivity_solid_to_liquid_ratio = \
+                thermal_conductivity_solid_to_liquid_ratio,
             **kwargs)
 
-        self.stefan_number = self.stefan_number.assign(0.125)
-
-        self.liquidus_temperature = self.liquidus_temperature.assign(0.)
-
-        self.density_solid_to_liquid_ratio = \
-            self.density_solid_to_liquid_ratio.assign(916.70/999.84)
-
-        self.heat_capacity_solid_to_liquid_ratio = \
-            self.heat_capacity_solid_to_liquid_ratio.assign(0.500)
-
-        self.thermal_conductivity_solid_to_liquid_ratio = \
-            self.thermal_conductivity_solid_to_liquid_ratio.assign(2.14/0.561)
-
-        self.cold_wall_temperature = self.cold_wall_temperature.assign(-1.)
+        self.cold_wall_temperature = self.cold_wall_temperature.assign(
+            cold_wall_temperature_during_freezing)
 

sapphire/benchmarks/heat_driven_cavity.py

@@ -16,27 +16,33 @@ def dirichlet_boundary_conditions(sim):
         fe.DirichletBC(W.sub(2), sim.hot_wall_temperature, 1),
         fe.DirichletBC(W.sub(2), sim.cold_wall_temperature, 2)]
 
-
+
+default_Ra = 1.e6
+
+default_Pr = 0.71
+
+default_Gr = default_Ra/default_Pr
+
 class Simulation(sapphire.simulations.navier_stokes_boussinesq.Simulation):
 
-    def __init__(self, *args, meshsize, **kwargs):
-
-        self.hot_wall_temperature = fe.Constant(0.5)
+    def __init__(self, *args, 
+            meshsize, 
+            hot_wall_temperature = 0.5,
+            cold_wall_temperature = -0.5,
+            grashof_number = default_Gr,
+            prandtl_number = default_Pr,
+            **kwargs):
+
+        self.hot_wall_temperature = fe.Constant(hot_wall_temperature)
 
-        self.cold_wall_temperature = fe.Constant(-0.5)
+        self.cold_wall_temperature = fe.Constant(cold_wall_temperature)
 
         super().__init__(
             *args,
             mesh = fe.UnitSquareMesh(meshsize, meshsize),
             initial_values = initial_values,
             dirichlet_boundary_conditions = dirichlet_boundary_conditions,
+            grashof_number = grashof_number,
+            prandtl_number = prandtl_number,
             **kwargs)
-
-        Ra = 1.e6
-
-        Pr = 0.71
-
-        self.grashof_number = self.grashof_number.assign(Ra/Pr)
-
-        self.prandtl_number = self.prandtl_number.assign(Pr)
 

sapphire/benchmarks/melt_octadecane_in_cavity.py

@@ -25,78 +25,37 @@ def dirichlet_boundary_conditions(sim):
 
     return [
         fe.DirichletBC(W.sub(1), (0., 0.), "on_boundary"),
-        fe.DirichletBC(W.sub(2), sim.hot_wall_temperature, 1),
+        fe.DirichletBC(W.sub(2), sim.hotwall_temperature, 1),
         fe.DirichletBC(W.sub(2), sim.initial_temperature, 2)]
 
 
 class Simulation(sapphire.simulations.\
         convection_coupled_phasechange.Simulation):
 
-    def __init__(self, *args, meshsize, initial_temperature = -0.01, **kwargs):
+    def __init__(self, *args, 
+            meshsize,
+            hotwall_temperature = 1.,
+            initial_temperature = -0.01, 
+            stefan_number = 0.045,
+            rayleigh_number = 3.27e5,
+            prandtl_number = 56.2,
+            liquidus_temperature = 0.,
+            **kwargs):
 
-        self.hot_wall_temperature = fe.Constant(1.)
+        self.hotwall_temperature = fe.Constant(hotwall_temperature)
 
         self.initial_temperature = fe.Constant(initial_temperature)
 
-        self.topwall_heatflux = fe.Constant(0.)
+        grashof_number = rayleigh_number/prandtl_number
 
         super().__init__(
             *args,
+            liquidus_temperature = liquidus_temperature,
+            stefan_number = stefan_number,
+            grashof_number = grashof_number,
+            prandtl_number = prandtl_number,
             mesh = fe.UnitSquareMesh(meshsize, meshsize),
             initial_values = initial_values,
             dirichlet_boundary_conditions = dirichlet_boundary_conditions,
             **kwargs)
-
-        q = self.topwall_heatflux
-
-        _, _, psi_T = fe.TestFunctions(self.function_space)
-
-        ds = fe.ds(domain = self.mesh, subdomain_id = 4)
-
-        self.variational_form_residual += psi_T*q*ds
-
-        Ra = 3.27e5
-
-        Pr = 56.2
-
-        self.grashof_number = self.grashof_number.assign(Ra/Pr)
-
-        self.prandtl_number = self.prandtl_number.assign(Pr)
-
-        self.stefan_number = self.stefan_number.assign(0.045)
-
-        self.liquidus_temperature = self.liquidus_temperature.assign(0.)
-
-    def run(self, *args,
-            endtime,
-            topwall_heatflux_poststart = -0.02,
-            topwall_heatflux_starttime = 40.,
-            **kwargs):
-
-        final_endtime = endtime
-
-        original_topwall_heatflux = self.topwall_heatflux.__float__()
-
-        if final_endtime < topwall_heatflux_starttime:
-
-            self.solutions, self.time = super().run(*args,
-                endtime = final_endtime,
-                **kwargs)
-
-            return self.solutions, self.time
-
-        self.solutions, self.time = super().run(*args,
-            endtime = topwall_heatflux_starttime,
-            write_initial_outputs = True,
-            **kwargs)
-
-        self.topwall_heatflux = self.topwall_heatflux.assign(
-            topwall_heatflux_poststart)
-
-        self.solutions, self.time = super().run(*args,
-            endtime = final_endtime,
-            write_initial_outputs = False,
-            **kwargs)
-
-        return self.solutions, self.time