Validate enthalpy-porosity for octadecane melting, extend with phase-…

…dependent material properties, verify extended model, validate for water freezing (#49)

Models:
* Primarily use solid velocity relaxation for enthalpy porosity model; also test general model with solid viscosity
* Extend enthalpy porosity model for phase-dependent material properties
* Extend enthalpy porosity implementation to 3D

Numerical Methods:
* Extend for BDF1-6
* Remove use of Taylor-Hood element for enthalpy method (and now only stabilize with pressure penalty)
* Generalize the continuation algorithm, which we now also use for high Grashof number in natural convection of water (which initializes the water freezing benchmark)

Verification:
* Include both second order and third order convergence tests

Validation:
* Add octadecane melting benchmark test
* Add heat flux to top wall to reproduce experimental result
* Validate for water freezing

Code Design:
* Stop using abc abstract class; just assert False for methods that must be redefined
* Always call fempy.model.Model.solve instead of fempy.model.Model.solver.solve
* Simplify handling of initial values and time discretization
* Add quadrature_degree and spatial_order parameters to model and temporal_order parameter to unsteady_model. 

Quality-of-Life:
* Add reporting feature, including counting SNES iterations and other auxiliary info
* Fix bug for Python 3.5 compatibility with PathLib

Other
* Increment version number

README.md

@@ -1,8 +1,9 @@
 # FemPy
-FemPy is a Python package using Firedrake to construct finite element models for monolithic multi-physics simulations.
+FemPy is an engine for constructing simulations
+based on mixed finite element models using Firedrake.
 
-## Setup
 
+## Setup
 [Install Firedrake](https://www.firedrakeproject.org/download.html).
 
 Download FemPy with 
@@ -11,15 +12,19 @@ Download FemPy with
 
 Test with
 
-    python3 -m pytest fempy/
+    python3 -m pytest -v -s -k "not long" fempy/
 
 Install with
 
     python3 fempy/setup.py install
 
+
 ## Development
+
+
 ### Project structure
 This project mostly follows the structure suggested by [The Hitchhiker's Guide to Python](http://docs.python-guide.org/en/latest/).
 
+
 ### Guidelines
 Mostly we try to follow PEP proposed guidelines, e.g. [The Zen of Python (PEP 20)](https://www.python.org/dev/peps/pep-0020/), and do not ever `from firedrake import *` ([PEP 8](https://www.python.org/dev/peps/pep-0008/)).

fempy/benchmarks/freeze_water.py

@@ -0,0 +1,183 @@
+import firedrake as fe
+import fempy.models.enthalpy_porosity
+import numpy
+
+
+class Model(fempy.models.enthalpy_porosity.Model):
+
+    def __init__(self, *args, spatial_dimensions, meshsize, **kwargs):
+
+        self.spatial_dimensions = spatial_dimensions
+
+        self.meshsize = meshsize
+
+        self.reference_temperature_range__degC = fe.Constant(10.)  # [deg C]
+
+        self.hot_wall_temperature = fe.Constant(1.)
+
+        self.cold_wall_temperature_before_freezing = fe.Constant(0.)
+
+        self.cold_wall_temperature_during_freezing = fe.Constant(-1.)
+
+        self.cold_wall_temperature = fe.Constant(0.)
+
+        self.cold_wall_temperature.assign(
+            self.cold_wall_temperature_before_freezing)
+
+        super().__init__(*args, **kwargs)
+
+        Ra = 2.518084e6
+
+        Pr = 6.99
+
+        self.grashof_number.assign(Ra/Pr)
+
+        self.prandtl_number.assign(Pr)
+
+        self.stefan_number.assign(0.125)
+
+        self.liquidus_temperature.assign(0.)
+
+        self.heat_capacity_solid_to_liquid_ratio.assign(0.500)
+
+        self.thermal_conductivity_solid_to_liquid_ratio.assign(2.14/0.561)
+
+    def buoyancy(self, T):
+        """ Eq. (25) from @cite{danaila2014newton} """
+        T_anomaly_degC = fe.Constant(4.0293)  # [deg C]
+
+        rho_anomaly_SI = fe.Constant(999.972)  # [kg/m^3]
+
+        w = fe.Constant(9.2793e-6)  # [(deg C)^(-q)]
+
+        q = fe.Constant(1.894816)
+
+        M = self.reference_temperature_range__degC
+
+        T_L = self.liquidus_temperature
+
+        def T_degC(T):
+            """ T = (T_degC - T_L)/M """
+            return M*T + T_L
+
+        def rho_of_T_degC(T_degC):
+            """ Eq. (24) from @cite{danaila2014newton} """
+            return rho_anomaly_SI*(1. - w*abs(T_degC - T_anomaly_degC)**q)
+
+        def rho(T):
+
+            return rho_of_T_degC(T_degC(T))
+
+        beta = fe.Constant(6.91e-5)  # [K^-1]
+
+        Gr = self.grashof_number
+
+        _, _, T = fe.split(self.solution)
+
+        ghat = fe.Constant(-self.unit_vectors()[1])
+
+        return Gr/(beta*M)*(rho(T_L) - rho(T))/rho(T_L)*ghat
+
+    def init_mesh(self):
+
+        if self.spatial_dimensions == 2:
+
+            Mesh = fe.UnitSquareMesh
+
+        elif self.spatial_dimensions == 3:
+
+            Mesh = fe.UnitCubeMesh
+
+        self.mesh = Mesh(*(self.meshsize,)*self.spatial_dimensions)
+
+    def initial_values(self):
+
+        print("Solving steady heat driven cavity to obtain initial values")
+
+        r = self.heat_driven_cavity_weak_form_residual()\
+            *self.integration_measure
+
+        u = self.solution
+
+        problem = fe.NonlinearVariationalProblem(
+            r, u, self.dirichlet_boundary_conditions, fe.derivative(r, u))
+
+        solver = fe.NonlinearVariationalSolver(
+            problem, 
+            solver_parameters = {
+                "snes_type": "newtonls",
+                "snes_max_it": 50,
+                "snes_monitor": True,
+                "ksp_type": "preonly", 
+                "pc_type": "lu", 
+                "mat_type": "aij",
+                "pc_factor_mat_solver_type": "mumps"})
+
+        x = fe.SpatialCoordinate(self.mesh)
+
+        T_c = self.cold_wall_temperature_before_freezing.__float__()
+
+        dim = self.mesh.geometric_dimension()
+
+        u.assign(fe.interpolate(
+            fe.Expression(
+                (0.,) + (0.,)*dim + (T_c,),
+                element = self.element),
+            self.function_space))
+
+        self.cold_wall_temperature.assign(
+            self.cold_wall_temperature_before_freezing)
+
+        T_h = self.hot_wall_temperature.__float__()
+
+        fempy.continuation.solve(
+            model = self,
+            solver = solver,
+            continuation_parameter = self.grashof_number,
+            continuation_sequence = None,
+            leftval = 0.,
+            rightval = self.grashof_number.__float__(),
+            startleft = True,
+            maxcount = 16)
+
+        self.cold_wall_temperature.assign(
+            self.cold_wall_temperature_during_freezing)
+
+        return self.solution
+
+    def init_dirichlet_boundary_conditions(self):
+
+        W = self.function_space
+
+        dim = self.mesh.geometric_dimension()
+
+        self.dirichlet_boundary_conditions = [
+            fe.DirichletBC(
+                W.sub(1), (0.,)*dim, "on_boundary"),
+            fe.DirichletBC(W.sub(2), self.hot_wall_temperature, 1),
+            fe.DirichletBC(W.sub(2), self.cold_wall_temperature, 2)]
+
+    def heat_driven_cavity_weak_form_residual(self):
+
+        mass = self.mass()
+
+        stabilization = self.stabilization()
+
+        p, u, T = fe.split(self.solution)
+
+        b = self.buoyancy(T)
+
+        Pr = self.prandtl_number
+
+        _, psi_u, psi_T = fe.TestFunctions(self.function_space)
+
+        inner, dot, grad, div, sym = \
+            fe.inner, fe.dot, fe.grad, fe.div, fe.sym
+
+        momentum = dot(psi_u, grad(u)*u + b) \
+            - div(psi_u)*p + 2.*inner(sym(grad(psi_u)), sym(grad(u)))
+
+        energy = psi_T*dot(u, grad(T)) + dot(grad(psi_T), 1./Pr*grad(T))
+
+        return mass + momentum + energy + stabilization
+

fempy/benchmarks/heat_driven_cavity.py

@@ -4,15 +4,17 @@
 
 class Model(fempy.models.navier_stokes_boussinesq.Model):
 
-    def __init__(self, meshsize):
+    def __init__(self, quadrature_degree, spatial_order, meshsize):
 
         self.meshsize = meshsize
 
         self.hot_wall_temperature = fe.Constant(0.5)
 
         self.cold_wall_temperature = fe.Constant(-0.5)
 
-        super().__init__()
+        super().__init__(
+            quadrature_degree = quadrature_degree,
+            spatial_order = spatial_order)
 
         Ra = 1.e6