diff --git a/README.md b/README.md index b0b41fc..664187f 100644 --- a/README.md +++ b/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/)). diff --git a/fempy/autosmooth.py b/fempy/autosmooth.py deleted file mode 100644 index c616c5d..0000000 --- a/fempy/autosmooth.py +++ /dev/null @@ -1,113 +0,0 @@ -""" Solve a strongly nonlinear problem by solving a sequence of over-regularized -problems with successively reduced regularization. -Here we use the word 'smooth' as a shorter synonym for 'regularize'. -""" -import firedrake as fe - - -def solve( - model, - maxval = 4., - firstval = None, - maxcount = 32): - - if model.smoothing_sequence is None: - - smoothing_sequence = (model.smoothing.__float__(),) - - if firstval is None: - - firstval = model.smoothing.__float__() - - while smoothing_sequence[0] < firstval: - - smoothing_sequence = (2.*smoothing_sequence[0],) + \ - smoothing_sequence - - else: - - smoothing_sequence = model.smoothing_sequence - - first_s_to_solve = smoothing_sequence[0] - - attempts = range(maxcount - len(smoothing_sequence)) - - solved = False - - for attempt in attempts: - - s_start_index = smoothing_sequence.index(first_s_to_solve) - - try: - - for s in smoothing_sequence[s_start_index:]: - - model.smoothing.assign(s) - - model.backup_solution.assign(model.solution) - - model.solver.solve() - - if not model.quiet: - - print("Solved with s = " + str(s)) - - solved = True - - break - - except fe.exceptions.ConvergenceError: - - current_s = model.smoothing.__float__() - - ss = smoothing_sequence - - if not model.quiet: - - print("Failed to solve with s = " + str(current_s) + - " from the sequence " + str(ss)) - - if attempt == attempts[-1]: - - break - - if current_s >= maxval: - - print("Exceeded maximum regularization (s_max = " + - str(maxval) + ")") - - break - - index = ss.index(current_s) - - if index == 0: - - s_to_insert = 2.*ss[0] - - new_ss = (s_to_insert,) + ss - - model.solution.assign(model.initial_values[0]) - - else: - - s_to_insert = (current_s + ss[index - 1])/2. - - new_ss = ss[:index] + (s_to_insert,) + ss[index:] - - model.solution.assign(model.backup_solution) - - smoothing_sequence = new_ss - - if not model.quiet: - - print("Inserted new value of " + str(s_to_insert)) - - first_s_to_solve = s_to_insert - - assert(solved) - - assert(model.smoothing.__float__() == - smoothing_sequence[-1]) - - model.smoothing_sequence = smoothing_sequence - \ No newline at end of file diff --git a/fempy/benchmarks/freeze_water.py b/fempy/benchmarks/freeze_water.py new file mode 100644 index 0000000..2976087 --- /dev/null +++ b/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 + \ No newline at end of file diff --git a/fempy/benchmarks/heat_driven_cavity.py b/fempy/benchmarks/heat_driven_cavity.py index 5f4980c..b635e67 100644 --- a/fempy/benchmarks/heat_driven_cavity.py +++ b/fempy/benchmarks/heat_driven_cavity.py @@ -4,7 +4,7 @@ class Model(fempy.models.navier_stokes_boussinesq.Model): - def __init__(self, meshsize): + def __init__(self, quadrature_degree, spatial_order, meshsize): self.meshsize = meshsize @@ -12,7 +12,9 @@ def __init__(self, meshsize): self.cold_wall_temperature = fe.Constant(-0.5) - super().__init__() + super().__init__( + quadrature_degree = quadrature_degree, + spatial_order = spatial_order) Ra = 1.e6 diff --git a/fempy/benchmarks/melt_octadecane.py b/fempy/benchmarks/melt_octadecane.py new file mode 100644 index 0000000..c71cbeb --- /dev/null +++ b/fempy/benchmarks/melt_octadecane.py @@ -0,0 +1,80 @@ +import firedrake as fe +import fempy.models.enthalpy_porosity + + +class Model(fempy.models.enthalpy_porosity.Model): + + def __init__(self, *args, meshsize, **kwargs): + + self.meshsize = meshsize + + self.hot_wall_temperature = fe.Constant(1.) + + self.cold_wall_temperature = fe.Constant(-0.01) + + self.topwall_heatflux = fe.Constant(0.) + + super().__init__(*args, **kwargs) + + Ra = 3.27e5 + + Pr = 56.2 + + self.grashof_number.assign(Ra/Pr) + + self.prandtl_number.assign(Pr) + + self.stefan_number.assign(0.045) + + self.liquidus_temperature.assign(0.) + + self.topwall_heatflux_postswitch = 0. + + self.topwall_heatflux_switchtime = 40. + 2.*self.time_tolerance + + def init_mesh(self): + + self.mesh = fe.UnitSquareMesh(self.meshsize, self.meshsize) + + def init_dirichlet_boundary_conditions(self): + + W = self.function_space + + self.dirichlet_boundary_conditions = [ + fe.DirichletBC(W.sub(1), (0., 0.), "on_boundary"), + fe.DirichletBC(W.sub(2), self.hot_wall_temperature, 1), + fe.DirichletBC(W.sub(2), self.cold_wall_temperature, 2)] + + def init_problem(self): + + q = self.topwall_heatflux + + _, _, psi_T = fe.TestFunctions(self.function_space) + + ds = fe.ds(domain = self.mesh, subdomain_id = 4) + + r = self.weak_form_residual*self.integration_measure + psi_T*q*ds + + u = self.solution + + self.problem = fe.NonlinearVariationalProblem( + r, u, self.dirichlet_boundary_conditions, fe.derivative(r, u)) + + def initial_values(self): + + return fe.interpolate( + fe.Expression( + (0., 0., 0., self.cold_wall_temperature.__float__()), + element = self.element), + self.function_space) + + def solve(self): + + if self.time.__float__() > \ + (self.topwall_heatflux_switchtime - self.time_tolerance): + + self.topwall_heatflux.assign( + self.topwall_heatflux_postswitch) + + super().solve() + \ No newline at end of file diff --git a/fempy/benchmarks/melting_octadecane.py b/fempy/benchmarks/melting_octadecane.py deleted file mode 100644 index 3623f34..0000000 --- a/fempy/benchmarks/melting_octadecane.py +++ /dev/null @@ -1,115 +0,0 @@ -import firedrake as fe -import fempy.models.enthalpy_porosity - - -class Model(fempy.models.enthalpy_porosity.Model): - - def __init__(self, meshsize): - - self.meshsize = meshsize - - self.hot_wall_temperature = fe.Constant(1.) - - self.cold_wall_temperature = fe.Constant(-0.01) - - super().__init__() - - Ra = 3.27e5 - - Pr = 56.2 - - self.grashof_number.assign(Ra/Pr) - - self.prandtl_number.assign(Pr) - - self.stefan_number.assign(0.045) - - self.liquidus_temperature.assign(0.) - - self.update_initial_values() - - def init_mesh(self): - - self.mesh = fe.UnitSquareMesh(self.meshsize, self.meshsize) - - def init_dirichlet_boundary_conditions(self): - - W = self.function_space - - self.dirichlet_boundary_conditions = [ - fe.DirichletBC(W.sub(1), (0., 0.), "on_boundary"), - fe.DirichletBC(W.sub(2), self.hot_wall_temperature, 1), - fe.DirichletBC(W.sub(2), self.cold_wall_temperature, 2)] - - def update_initial_values(self): - - initial_values = fe.interpolate( - fe.Expression( - (0., 0., 0., self.cold_wall_temperature.__float__()), - element = self.element), - self.function_space) - - self.initial_values.assign(initial_values) - - -class ModelWithDarcyResistance( - fempy.models.enthalpy_porosity.ModelWithDarcyResistance): - - def __init__(self, meshsize): - - self.meshsize = meshsize - - self.hot_wall_temperature = fe.Constant(1.) - - self.cold_wall_temperature = fe.Constant(-0.01) - - super().__init__() - - Ra = 3.27e5 - - Pr = 56.2 - - self.grashof_number.assign(Ra/Pr) - - self.prandtl_number.assign(Pr) - - self.stefan_number.assign(0.045) - - self.liquidus_temperature.assign(0.) - - self.update_initial_values() - - def init_mesh(self): - - self.mesh = fe.UnitSquareMesh(self.meshsize, self.meshsize) - - def init_dirichlet_boundary_conditions(self): - - W = self.function_space - - self.dirichlet_boundary_conditions = [ - fe.DirichletBC(W.sub(1), (0., 0.), "on_boundary"), - fe.DirichletBC(W.sub(2), self.hot_wall_temperature, 1), - fe.DirichletBC(W.sub(2), self.cold_wall_temperature, 2)] - - def update_initial_values(self): - - initial_values = fe.interpolate( - fe.Expression( - (0., 0., 0., self.cold_wall_temperature.__float__()), - element = self.element), - self.function_space) - - self.initial_values.assign(initial_values) - - -if __name__ == "__main__": - - model = Model(meshsize = 32) - - model.timestep_size.assign(10.) - - model.latent_heat_smoothing.assign(1./256.) - - model.run(endtime = 70., plot = True) - \ No newline at end of file diff --git a/fempy/continuation.py b/fempy/continuation.py new file mode 100644 index 0000000..82bc5f0 --- /dev/null +++ b/fempy/continuation.py @@ -0,0 +1,115 @@ +import firedrake as fe +import fempy.model + + +def solve( + model, + solver, + continuation_parameter, + continuation_sequence, + leftval, + rightval, + startleft = False, + maxcount = 32): + """ Solve a strongly nonlinear problem + by solving a sequence of over-regularized problems + with successively reduced regularization. + Here we use the word 'smooth' as a short synonym for 'regularize'. + + Always continue from left to right. + """ + if continuation_sequence is None: + + my_continuation_sequence = (leftval, rightval) + + else: + + my_continuation_sequence = continuation_sequence + + if startleft: + + first_s_to_solve = my_continuation_sequence[0] + + else: + + first_s_to_solve = my_continuation_sequence[-1] + + attempts = range(maxcount - len(my_continuation_sequence)) + + solved = False + + def bounded(val): + + if leftval < rightval: + + return leftval <= val and val <= rightval + + if leftval > rightval: + + return leftval >= val and val >= rightval + + for attempt in attempts: + + s_start_index = my_continuation_sequence.index(first_s_to_solve) + + try: + + for s in my_continuation_sequence[s_start_index:]: + + continuation_parameter.assign(s) + + model.backup_solution.assign(model.solution) + + solver.solve() + + if not model.quiet: + + print("Solved with continuation parameter = " + str(s)) + + solved = True + + break + + except fe.exceptions.ConvergenceError: + + current_s = continuation_parameter.__float__() + + ss = my_continuation_sequence + + if not model.quiet: + + print("Failed to solve with continuation paramter = " + + str(current_s) + + " from the sequence " + str(ss)) + + if attempt == attempts[-1]: + + break + + index = ss.index(current_s) + + assert(index > 0) + + s_to_insert = (current_s + ss[index - 1])/2. + + assert(bounded(s_to_insert)) + + new_ss = ss[:index] + (s_to_insert,) + ss[index:] + + model.solution.assign(model.backup_solution) + + my_continuation_sequence = new_ss + + if not model.quiet: + + print("Inserted new value of " + str(s_to_insert)) + + first_s_to_solve = s_to_insert + + assert(solved) + + assert(continuation_parameter.__float__() == + my_continuation_sequence[-1]) + + return my_continuation_sequence + \ No newline at end of file diff --git a/fempy/mms.py b/fempy/mms.py index c5c0361..08ad764 100644 --- a/fempy/mms.py +++ b/fempy/mms.py @@ -12,12 +12,18 @@ def make_mms_verification_model_class(Model): class MMSVerificationModel(Model): - def init_mesh(self): + def __init__(self, *args, **kwargs): - super().init_mesh() + super().__init__(*args, **kwargs) + self._L2_error = None + + def init_solutions(self): + self.init_manufactured_solution() + super().init_solutions() + def init_weak_form_residual(self): super().init_weak_form_residual() @@ -48,6 +54,26 @@ def init_dirichlet_boundary_conditions(self): fe.DirichletBC(V, g, "on_boundary") for V, g in zip(W, u_m)] + def initial_values(self): + + initial_values = fe.Function(self.function_space) + + if (type(self.manufactured_solution) is not type([0,])) \ + and (type(self.manufactured_solution) is not type((0,))): + + w_m = (self.manufactured_solution,) + + else: + + w_m = self.manufactured_solution + + for u_m, V in zip( + w_m, self.function_space): + + initial_values.assign(fe.interpolate(u_m, V)) + + return initial_values + def L2_error(self): dx = self.integration_measure @@ -71,6 +97,12 @@ def L2_error(self): return e + def report(self, write_header = True): + + self._L2_error = self.L2_error() + + super().report(write_header = write_header) + return MMSVerificationModel @@ -79,12 +111,14 @@ def verify_spatial_order_of_accuracy( expected_order, mesh_sizes, tolerance, + constructor_kwargs = {}, parameters = {}, timestep_size = None, endtime = None, starttime = 0., plot_errors = False, - plot_solution = False): + plot_solution = False, + report = False): MMSVerificationModel = make_mms_verification_model_class(Model) @@ -94,10 +128,10 @@ def verify_spatial_order_of_accuracy( for meshsize in mesh_sizes: - model = MMSVerificationModel(meshsize = meshsize) + model = MMSVerificationModel(**constructor_kwargs, meshsize = meshsize) model.output_directory_path = model.output_directory_path.joinpath( - "mms_space") + "mms_space_p" + str(expected_order) + "/") if timestep_size is not None: @@ -116,7 +150,7 @@ def verify_spatial_order_of_accuracy( model.timestep_size.assign(timestep_size) - model.run(endtime = endtime, plot = plot_solution) + model.run(endtime = endtime, plot = plot_solution, report = report) else: @@ -156,11 +190,8 @@ def verify_spatial_order_of_accuracy( plt.grid(True) - outdir_path = pathlib.Path("output/mms/") - - outdir_path.mkdir(parents = True, exist_ok = True) - - filepath = outdir_path.joinpath("e_vs_h").with_suffix(".png") + filepath = model.output_directory_path.joinpath( + "e_vs_h").with_suffix(".png") print("Writing plot to " + str(filepath)) @@ -178,36 +209,31 @@ def verify_temporal_order_of_accuracy( timestep_sizes, endtime, tolerance, + constructor_kwargs = {}, parameters = {}, starttime = 0., plot_errors = False, - plot_solution = False): + plot_solution = False, + report = False): MMSVerificationModel = make_mms_verification_model_class(Model) table = fempy.table.Table(("Delta_t", "L2_error", "temporal_order")) - model = MMSVerificationModel(meshsize = meshsize) + model = MMSVerificationModel(**constructor_kwargs, meshsize = meshsize) basepath = model.output_directory_path model.assign_parameters(parameters) - u_m = model.initial_values - - if not ((type(u_m) == type((0,))) or (type(u_m) == type([0,]))): - - u_m = (u_m), - - initial_values = [fe.Function(iv) for iv in u_m] - print("") for timestep_size in timestep_sizes: model.timestep_size.assign(timestep_size) - model.output_directory_path = basepath.joinpath("mms_time") + model.output_directory_path = basepath.joinpath( + "mms_time_q" + str(expected_order) + "/") model.output_directory_path = model.output_directory_path.joinpath( "m" + str(meshsize)) @@ -217,19 +243,9 @@ def verify_temporal_order_of_accuracy( model.time.assign(starttime) - if (type(model.initial_values) == type((0,)) or - (type(model.initial_values) == type([0,]))): + model.update_initial_values() - for i, iv in enumerate(model.initial_values): - - iv.assign(initial_values[i]) - else: - - model.initial_values.assign(initial_values[0]) - - model.solution.assign(model.initial_values[0]) - - model.run(endtime = endtime, plot = plot_solution) + model.run(endtime = endtime, plot = plot_solution, report = report) table.append({ "Delta_t": timestep_size, @@ -265,11 +281,8 @@ def verify_temporal_order_of_accuracy( plt.grid(True) - outdir_path = pathlib.Path("output/mms/") - - outdir_path.mkdir(parents = True, exist_ok = True) - - filepath = outdir_path.joinpath("e_vs_Delta_t").with_suffix(".png") + filepath = model.output_directory_path.joinpath( + "e_vs_Delta_t").with_suffix(".png") print("Writing plot to " + str(filepath)) diff --git a/fempy/model.py b/fempy/model.py index ac78f1d..374b335 100644 --- a/fempy/model.py +++ b/fempy/model.py @@ -1,13 +1,16 @@ """ An abstract class on which to base finite element models """ import firedrake as fe -import abc import pathlib import matplotlib.pyplot as plt -class Model(metaclass = abc.ABCMeta): - """ An abstract class on which to base finite element models. """ - def __init__(self): +class Model(object): + """ A class on which to base finite element models. """ + def __init__(self, quadrature_degree, spatial_order): + + self.quadrature_degree = quadrature_degree + + self.spatial_order = spatial_order self.init_mesh() @@ -15,7 +18,7 @@ def __init__(self): self.init_function_space() - self.init_solution() + self.init_solutions() self.init_integration_measure() @@ -31,30 +34,34 @@ def __init__(self): self.output_directory_path = pathlib.Path("output/") - @abc.abstractmethod + self.snes_iteration_counter = 0 + def init_mesh(self): """ Redefine this to set `self.mesh` to a `fe.Mesh`. """ + assert(False) - @abc.abstractmethod def init_element(self): """ Redefine this to set `self.element` to a `fe.FiniteElement` or `fe.MixedElement`. """ + assert(False) - @abc.abstractmethod def init_weak_form_residual(self): """ Redefine this to set `self.weak_form_residual` to a `fe.NonlinearVariationalForm`. """ + assert(False) def init_function_space(self): self.function_space = fe.FunctionSpace(self.mesh, self.element) - def init_solution(self): + def init_solutions(self): - self.solution = fe.Function(self.function_space) + self.solutions = [fe.Function(self.function_space)] + + self.solution = self.solutions[0] def init_dirichlet_boundary_conditions(self): """ Optionallay redefine this @@ -64,7 +71,7 @@ def init_dirichlet_boundary_conditions(self): def init_integration_measure(self): - self.integration_measure = fe.dx + self.integration_measure = fe.dx(degree = self.quadrature_degree) def init_problem(self): @@ -103,6 +110,8 @@ def assign_parameters(self, parameters): def solve(self): self.solver.solve() + + self.snes_iteration_counter += self.solver.snes.getIterationNumber() def unit_vectors(self): @@ -110,6 +119,10 @@ def unit_vectors(self): return tuple([fe.unit_vector(i, dim) for i in range(dim)]) + def write_solution(self, file): + + file.write(*self.solution.split(), time = self.time) + def plot(self): for i, f in enumerate(self.solution.split()): diff --git a/fempy/models/convection_diffusion.py b/fempy/models/convection_diffusion.py index c267eb3..7a7b291 100644 --- a/fempy/models/convection_diffusion.py +++ b/fempy/models/convection_diffusion.py @@ -5,11 +5,13 @@ class Model(fempy.model.Model): - def __init__(self): + def __init__(self, quadrature_degree, spatial_order): self.kinematic_viscosity = fe.Constant(1.) - super().__init__() + super().__init__( + quadrature_degree = quadrature_degree, + spatial_order = spatial_order) def init_element(self): diff --git a/fempy/models/enthalpy.py b/fempy/models/enthalpy.py index 0b0b2d5..0432c13 100644 --- a/fempy/models/enthalpy.py +++ b/fempy/models/enthalpy.py @@ -5,7 +5,7 @@ class Model(fempy.unsteady_model.Model): - def __init__(self): + def __init__(self, quadrature_degree, spatial_order, temporal_order): self.stefan_number = fe.Constant(1.) @@ -13,11 +13,15 @@ def __init__(self): self.smoothing = fe.Constant(1./32.) - super().__init__() + super().__init__( + quadrature_degree = quadrature_degree, + spatial_order = spatial_order, + temporal_order = temporal_order) def init_element(self): - self.element = fe.FiniteElement("P", self.mesh.ufl_cell(), 1) + self.element = fe.FiniteElement( + "P", self.mesh.ufl_cell(), self.spatial_order - 1) def porosity(self, T): @@ -30,20 +34,15 @@ def porosity(self, T): return 0.5*(1. + tanh((T - T_L)/s)) def init_time_discrete_terms(self): - """ Implicit Euler finite difference scheme """ - T = self.solution - - T_n = self.initial_values - - Delta_t = self.timestep_size - - T_t = (T - T_n)/Delta_t - phil = self.porosity + super().init_time_discrete_terms() - phil_t = (phil(T) - phil(T_n))/Delta_t + phil_t = fempy.time_discretization.bdf( + [self.porosity(T_n) for T_n in self.solutions], + order = self.temporal_order, + timestep_size = self.timestep_size) - self.time_discrete_terms = T_t, phil_t + self.time_discrete_terms = (self.time_discrete_terms, phil_t) def init_weak_form_residual(self): diff --git a/fempy/models/enthalpy_porosity.py b/fempy/models/enthalpy_porosity.py index 995b5fd..dbf7320 100644 --- a/fempy/models/enthalpy_porosity.py +++ b/fempy/models/enthalpy_porosity.py @@ -2,16 +2,12 @@ import firedrake as fe import fempy.unsteady_model import matplotlib.pyplot as plt -import fempy.autosmooth +import fempy.continuation class Model(fempy.unsteady_model.Model): - def __init__(self): - - self.liquid_dynamic_viscosity = fe.Constant(1.) - - self.solid_dynamic_viscosity = fe.Constant(1.e8) + def __init__(self, *args, **kwargs): self.grashof_number = fe.Constant(1.) @@ -23,28 +19,39 @@ def __init__(self): self.liquidus_temperature = fe.Constant(0.) + self.heat_capacity_solid_to_liquid_ratio = fe.Constant(1.) + + self.thermal_conductivity_solid_to_liquid_ratio = fe.Constant(1.) + + self.solid_velocity_relaxation_factor = fe.Constant(1.e-12) + self.smoothing = fe.Constant(1./256.) self.smoothing_sequence = None + self.save_smoothing_sequence = False + self.autosmooth_enable = True self.autosmooth_maxval = 4. - self.autosmooth_firstval = 1./4. - self.autosmooth_maxcount = 32 - super().__init__() + super().__init__(*args, **kwargs) self.backup_solution = fe.Function(self.solution) + # Initialize some attributes to be reported + self.liquid_area = None + def init_element(self): - + + rx = self.spatial_order + self.element = fe.MixedElement( - fe.FiniteElement("P", self.mesh.ufl_cell(), 1), - fe.VectorElement("P", self.mesh.ufl_cell(), 2), - fe.FiniteElement("P", self.mesh.ufl_cell(), 1)) + fe.FiniteElement("P", self.mesh.ufl_cell(), rx - 1), + fe.VectorElement("P", self.mesh.ufl_cell(), rx - 1), + fe.FiniteElement("P", self.mesh.ufl_cell(), rx - 1)) def porosity(self, T): """ Regularization from @cite{zimmerman2018monolithic} """ @@ -56,51 +63,68 @@ def porosity(self, T): return 0.5*(1. + tanh((T - T_L)/s)) + def phase_dependent_material_property(self, solid_to_liquid_ratio): + + a_s = solid_to_liquid_ratio + + def a(phil): + + return a_s + (1. - a_s)*phil + + return a + def buoyancy(self, T): """ Boussinesq buoyancy """ _, _, T = fe.split(self.solution) Gr = self.grashof_number - _, jhat = self.unit_vectors() - - ghat = fe.Constant(-jhat) + ghat = fe.Constant(-self.unit_vectors()[1]) return Gr*T*ghat - def dynamic_viscosity(self): - - _, _, T = fe.split(self.solution) + def solid_velocity_relaxation(self, T): - mu_s = self.solid_dynamic_viscosity + phil = self.porosity(T) - mu_l = self.liquid_dynamic_viscosity + phis = 1. - phil - phil = self.porosity(T) + tau = self.solid_velocity_relaxation_factor - return mu_s + (mu_l - mu_s)*phil + return 1./tau*phis def init_time_discrete_terms(self): - """ Implicit Euler finite difference scheme """ - p, u, T = fe.split(self.solution) - p_n, u_n, T_n = fe.split(self.initial_values) + super().init_time_discrete_terms() - Delta_t = self.timestep_size + temperature_solutions = [] - u_t = (u - u_n)/Delta_t + for solution in self.solutions: - T_t = (T - T_n)/Delta_t + temperature_solutions.append(fe.split(solution)[2]) phil = self.porosity - phil_t = (phil(T) - phil(T_n))/Delta_t + cp = self.phase_dependent_material_property( + self.heat_capacity_solid_to_liquid_ratio) + + cpT_t = fempy.time_discretization.bdf( + [cp(phil(T))*T for T in temperature_solutions], + order = self.temporal_order, + timestep_size = self.timestep_size) + + cpphil_t = fempy.time_discretization.bdf( + [cp(phil(T))*self.porosity(T) for T in temperature_solutions], + order = self.temporal_order, + timestep_size = self.timestep_size) - self.time_discrete_terms = u_t, T_t, phil_t + for w_i_t in (cpT_t, cpphil_t): + + self.time_discrete_terms.append(w_i_t) def mass(self): - p, u, _ = fe.split(self.solution) + _, u, _ = fe.split(self.solution) psi_p, _, _ = fe.TestFunctions(self.function_space) @@ -108,29 +132,25 @@ def mass(self): mass = psi_p*div(u) - gamma = self.pressure_penalty_factor - - stabilization = gamma*psi_p*p - - return mass + stabilization + return mass def momentum(self): p, u, T = fe.split(self.solution) - u_t, _, _ = self.time_discrete_terms + _, u_t, _, _, _ = self.time_discrete_terms b = self.buoyancy(T) - mu = self.dynamic_viscosity() + d = self.solid_velocity_relaxation(T) _, psi_u, _ = fe.TestFunctions(self.function_space) inner, dot, grad, div, sym = \ fe.inner, fe.dot, fe.grad, fe.div, fe.sym - return dot(psi_u, u_t + grad(u)*u + b) \ - - div(psi_u)*p + 2.*inner(sym(grad(psi_u)), mu*sym(grad(u))) + return dot(psi_u, u_t + grad(u)*u + b + d*u) \ + - div(psi_u)*p + 2.*inner(sym(grad(psi_u)), sym(grad(u))) def enthalpy(self): @@ -140,14 +160,32 @@ def enthalpy(self): _, u, T = fe.split(self.solution) - _, T_t, phil_t = self.time_discrete_terms + phil = self.porosity(T) + + cp = self.phase_dependent_material_property( + self.heat_capacity_solid_to_liquid_ratio)(phil) + + k = self.phase_dependent_material_property( + self.thermal_conductivity_solid_to_liquid_ratio)(phil) + + _, _, _, cpT_t, cpphil_t = self.time_discrete_terms _, _, psi_T = fe.TestFunctions(self.function_space) dot, grad = fe.dot, fe.grad - return psi_T*(T_t + 1./Ste*phil_t) \ - + dot(grad(psi_T), 1./Pr*grad(T) - T*u) + return psi_T*(cpT_t + dot(u, cp*grad(T)) + 1./Ste*cpphil_t) \ + + dot(grad(psi_T), k/Pr*grad(T)) + + def stabilization(self): + + p, _, _ = fe.split(self.solution) + + psi_p, _, _ = fe.TestFunctions(self.function_space) + + gamma = self.pressure_penalty_factor + + return gamma*psi_p*p def init_weak_form_residual(self): """ Weak form from @cite{zimmerman2018monolithic} """ @@ -157,24 +195,31 @@ def init_weak_form_residual(self): enthalpy = self.enthalpy() - self.weak_form_residual = mass + momentum + enthalpy - - def init_integration_measure(self): - - self.integration_measure = fe.dx(degree = 4) + stabilization = self.stabilization() + self.weak_form_residual = mass + momentum + enthalpy + stabilization + def solve(self): if self.autosmooth_enable: - fempy.autosmooth.solve(self, - firstval = self.autosmooth_firstval, - maxval = self.autosmooth_maxval, + smoothing_sequence = fempy.continuation.solve( + model = self, + solver = self.solver, + continuation_parameter = self.smoothing, + continuation_sequence = self.smoothing_sequence, + leftval = self.autosmooth_maxval, + rightval = self.smoothing.__float__(), + startleft = True, maxcount = self.autosmooth_maxcount) + + if self.save_smoothing_sequence: + + self.smoothing_sequence = smoothing_sequence elif self.smoothing_sequence == None: - self.solver.solve() + super().solve() else: @@ -185,12 +230,22 @@ def solve(self): self.smoothing.assign(s) - self.solver.solve() + super().solve() if not self.quiet: print("Solved with s = " + str(s)) + def report(self, write_header): + + p, u, T = self.solution.split() + + phil = self.porosity(T) + + self.liquid_area = fe.assemble(phil*self.integration_measure) + + super().report(write_header = write_header) + def plot(self): V = fe.FunctionSpace( @@ -203,9 +258,9 @@ def plot(self): timestr = str(self.time.__float__()) for f, label, filename in zip( - (self.mesh, p, u, T, phil), - ("\\Omega_h", "p", "\\mathbf{u}", "T", "\\phi_l"), - ("mesh", "p", "u", "T", "phil")): + (p, u, T, phil), + ("p", "\\mathbf{u}", "T", "\\phi_l"), + ("p", "u", "T", "phil")): fe.plot(f) @@ -229,80 +284,3 @@ def plot(self): plt.savefig(str(filepath)) plt.close() - - -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 - - def bdf2(u_np1, u_n, u_nm1): - - return (3.*u_np1 - 4.*u_n + u_nm1)/(2.*Delta_t) - - p_np1, u_np1, T_np1 = fe.split(self.solution) - - p_n, u_n, T_n = fe.split(self.initial_values[0]) - - p_nm1, u_nm1, T_nm1 = fe.split(self.initial_values[1]) - - u_t = bdf2(u_np1, u_n, u_nm1) - - T_t = bdf2(T_np1, T_n, T_nm1) - - phil = self.porosity - - phil_t = bdf2(phil(T_np1), phil(T_n), phil(T_nm1)) - - self.time_discrete_terms = u_t, T_t, phil_t - - -class ModelWithDarcyResistance(Model): - - def __init__(self): - - self.darcy_resistance_factor = fe.Constant(1.e6) - - self.small_number_to_avoid_division_by_zero = fe.Constant(1.e-8) - - super().__init__() - - delattr(self, "solid_dynamic_viscosity") - - def darcy_resistance(self, T): - """ Resistance to flow based on permeability of the porous media """ - D = self.darcy_resistance_factor - - epsilon = self.small_number_to_avoid_division_by_zero - - phil = self.porosity(T) - - return D*(1. - phil)**2/(phil**3 + epsilon) - - def dynamic_viscosity(self): - - return self.liquid_dynamic_viscosity - - def momentum(self): - - _, u, T = fe.split(self.solution) - - u_t, _, _ = self.time_discrete_terms - - b = self.buoyancy(T) - - d = self.darcy_resistance(T) - - _, psi_u, _ = fe.TestFunctions(self.function_space) - - dot = fe.dot - - return super().momentum() + dot(psi_u, d*u) - - \ No newline at end of file diff --git a/fempy/models/heat.py b/fempy/models/heat.py index af0ad44..1ccd41f 100644 --- a/fempy/models/heat.py +++ b/fempy/models/heat.py @@ -7,17 +7,8 @@ class Model(fempy.unsteady_model.Model): def init_element(self): - self.element = fe.FiniteElement("P", self.mesh.ufl_cell(), 1) - - def init_time_discrete_terms(self): - """ Implicit Euler finite difference scheme """ - u = self.solution - - un = self.initial_values - - Delta_t = self.timestep_size - - self.time_discrete_terms = (u - un)/Delta_t + self.element = fe.FiniteElement( + "P", self.mesh.ufl_cell(), self.spatial_order - 1) def init_weak_form_residual(self): diff --git a/fempy/models/laplace.py b/fempy/models/laplace.py index 872a9dc..209012b 100644 --- a/fempy/models/laplace.py +++ b/fempy/models/laplace.py @@ -7,7 +7,8 @@ class Model(fempy.model.Model): def init_element(self): - self.element = fe.FiniteElement("P", self.mesh.ufl_cell(), 1) + self.element = fe.FiniteElement( + "P", self.mesh.ufl_cell(), self.spatial_order - 1) def init_weak_form_residual(self): diff --git a/fempy/models/navier_stokes.py b/fempy/models/navier_stokes.py index 9ed539d..12ab897 100644 --- a/fempy/models/navier_stokes.py +++ b/fempy/models/navier_stokes.py @@ -5,15 +5,19 @@ class Model(fempy.model.Model): - def __init__(self): + def __init__(self, quadrature_degree, spatial_order): - super().__init__() + super().__init__( + quadrature_degree = quadrature_degree, + spatial_order = spatial_order) def init_element(self): self.element = fe.MixedElement( - fe.VectorElement('P', self.mesh.ufl_cell(), 2), - fe.FiniteElement('P', self.mesh.ufl_cell(), 1)) + fe.VectorElement( + 'P', self.mesh.ufl_cell(), self.spatial_order), + fe.FiniteElement( + 'P', self.mesh.ufl_cell(), self.spatial_order - 1)) def init_weak_form_residual(self): diff --git a/fempy/models/navier_stokes_boussinesq.py b/fempy/models/navier_stokes_boussinesq.py index 5cb1ada..ab39ee8 100644 --- a/fempy/models/navier_stokes_boussinesq.py +++ b/fempy/models/navier_stokes_boussinesq.py @@ -5,26 +5,27 @@ class Model(fempy.model.Model): - def __init__(self): - - self.dynamic_viscosity = fe.Constant(1.) + def __init__(self, quadrature_degree, spatial_order): self.grashof_number = fe.Constant(1.) self.prandtl_number = fe.Constant(1.) - super().__init__() + super().__init__( + quadrature_degree = quadrature_degree, + spatial_order = spatial_order) def init_element(self): self.element = fe.MixedElement( - fe.FiniteElement("P", self.mesh.ufl_cell(), 1), - fe.VectorElement("P", self.mesh.ufl_cell(), 2), - fe.FiniteElement("P", self.mesh.ufl_cell(), 1)) + fe.FiniteElement( + "P", self.mesh.ufl_cell(), self.spatial_order - 1), + fe.VectorElement( + "P", self.mesh.ufl_cell(), self.spatial_order), + fe.FiniteElement( + "P", self.mesh.ufl_cell(), self.spatial_order - 1)) def init_weak_form_residual(self): - - mu = self.dynamic_viscosity Gr = self.grashof_number @@ -46,7 +47,7 @@ def init_weak_form_residual(self): mass = psi_p*div(u) momentum = dot(psi_u, grad(u)*u + Gr*T*ghat) \ - - div(psi_u)*p + 2.*mu*inner(sym(grad(psi_u)), sym(grad(u))) + - 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)) diff --git a/fempy/models/unsteady_navier_stokes.py b/fempy/models/unsteady_navier_stokes.py index 8b1d084..2bde820 100644 --- a/fempy/models/unsteady_navier_stokes.py +++ b/fempy/models/unsteady_navier_stokes.py @@ -8,20 +8,10 @@ class Model(fempy.unsteady_model.Model): def init_element(self): self.element = fe.MixedElement( - fe.VectorElement('P', self.mesh.ufl_cell(), 2), - fe.FiniteElement('P', self.mesh.ufl_cell(), 1)) - - def init_time_discrete_terms(self): - """ Implicit Euler finite difference scheme """ - u, p = fe.split(self.solution) - - u_n, p_n = fe.split(self.initial_values) - - Delta_t = self.timestep_size - - u_t = (u - u_n)/Delta_t - - self.time_discrete_terms = u_t + fe.VectorElement( + 'P', self.mesh.ufl_cell(), self.spatial_order), + fe.FiniteElement( + 'P', self.mesh.ufl_cell(), self.spatial_order - 1)) def init_weak_form_residual(self): @@ -30,7 +20,7 @@ def init_weak_form_residual(self): u, p = fe.split(self.solution) - u_t = self.time_discrete_terms + u_t, _ = self.time_discrete_terms psi_u, psi_p = fe.TestFunctions(self.solution.function_space()) diff --git a/fempy/time_discretization.py b/fempy/time_discretization.py new file mode 100644 index 0000000..18524a3 --- /dev/null +++ b/fempy/time_discretization.py @@ -0,0 +1,48 @@ +""" Time discretization formulas """ + + +def bdf(ys, order, timestep_size): + """ Backward difference formulas """ + assert(len(ys) == (order + 1)) + + """ Table of BDF method coefficients """ + if order == 1: + + alphas = (1., -1.) + + elif order == 2: + + alphas = (3./2., -2., 1./2.) + + elif order == 3: + + alphas = (11./6., -3., 3./2., -1./3.) + + elif order == 4: + + alphas = (25./12., -4., 3., -4./3., 1./4.) + + elif order == 5: + + alphas = (137./60., -5., 5., -10./3., 5./4., -1./5.) + + elif order == 6: + + alphas = (147./60., -6., 15./2., -20./3., 15./4., -6./5., 1./6.) + + else: + + raise("BDF is not zero-stable with order > 6.") + + + y_t = alphas[-1]*ys[-1] + + for alpha, y in zip(alphas[:-1], ys[:-1]): + + y_t += alpha*y + + y_t /= timestep_size + + + return y_t + \ No newline at end of file diff --git a/fempy/unsteady_model.py b/fempy/unsteady_model.py index 8b63a42..7179613 100644 --- a/fempy/unsteady_model.py +++ b/fempy/unsteady_model.py @@ -1,17 +1,18 @@ """ An abstract class on which to base finite element models -with auxiliary data for unsteady (i.e. time-dependent) simulations. +with auxiliary data for unsteady simulations. """ import firedrake as fe import fempy.model -import abc +import fempy.time_discretization import matplotlib.pyplot as plt +import csv class Model(fempy.model.Model): """ An abstract class on which to base finite element models with auxiliary data for unsteady (i.e. time-dependent) simulations. """ - def __init__(self): + def __init__(self, quadrature_degree, spatial_order, temporal_order): self.time = fe.Constant(0.) @@ -19,50 +20,108 @@ def __init__(self): self.time_tolerance = 1.e-8 - super().__init__() + self.temporal_order = temporal_order + + super().__init__( + quadrature_degree = quadrature_degree, + spatial_order = spatial_order) self.solution_file = None - def init_initial_values(self): + def init_solutions(self): + + super().init_solutions() + + for i in range(self.temporal_order): - self.initial_values = fe.Function(self.function_space) + self.solutions.append(fe.Function(self.function_space)) + + self.init_time_discrete_terms() - def init_solution(self): + def update_initial_values(self): - super().init_solution() + initial_values = self.initial_values() - self.init_initial_values() + for solution in self.solutions: - self.init_time_discrete_terms() + solution.assign(initial_values) - def push_back_initial_values(self): + def initial_values(self): + """ Redefine this to return a fe.Function(self.function_space) + with the initial values. + """ + assert(False) + + def init_time_discrete_terms(self): - if not((type(self.initial_values) == type((0,))) - or (type(self.initial_values) == type([0,]))): + solutions = self.solutions - self.initial_values.assign(self.solution) - + time_discrete_terms = [ + fempy.time_discretization.bdf( + [fe.split(solutions[n])[i] for n in range(len(solutions))], + order = self.temporal_order, + timestep_size = self.timestep_size) + for i in range(len(fe.split(solutions[0])))] + + if len(time_discrete_terms) == 1: + + self.time_discrete_terms = time_discrete_terms[0] + else: - for i in range(len(self.initial_values) - 1): + self.time_discrete_terms = time_discrete_terms + + def push_back_solutions(self): + + for i in range(len(self.solutions[1:])): + + self.solutions[-(i + 1)].assign( + self.solutions[-(i + 2)]) - self.initial_values[-i - 1].assign( - self.initial_values[-i - 2]) - - self.initial_values[0].assign(self.solution) + def run(self, + endtime, + report = True, + write_solution = False, + plot = False, + update_initial_values = True): + + self.output_directory_path.mkdir( + parents = True, exist_ok = True) + + solution_filepath = self.output_directory_path.joinpath("solution").with_suffix(".pvd") + + solution_file = fe.File(str(solution_filepath)) + + if update_initial_values: + + self.update_initial_values() + + if report: + + self.report(write_header = True) + + if write_solution: - def run(self, endtime, plot = False): + self.write_solution(solution_file) if plot: - + self.plot() - + while self.time.__float__() < (endtime - self.time_tolerance): self.time.assign(self.time + self.timestep_size) self.solve() + if report: + + self.report(write_header = False) + + if write_solution: + + self.write_solution(solution_file) + if plot: self.plot() @@ -72,12 +131,33 @@ def run(self, endtime, plot = False): self.solution_file.write( *self.solution.split(), time = self.time.__float__()) - self.push_back_initial_values() + self.push_back_solutions() if not self.quiet: print("Solved at time t = " + str(self.time.__float__())) + def report(self, write_header = True): + + repvars = vars(self).copy() + + for key in repvars.keys(): + + if type(repvars[key]) is type(fe.Constant(0.)): + + repvars[key] = repvars[key].__float__() + + with self.output_directory_path.joinpath( + "report").with_suffix(".csv").open("a+") as csv_file: + + writer = csv.DictWriter(csv_file, fieldnames = repvars.keys()) + + if write_header: + + writer.writeheader() + + writer.writerow(repvars) + def plot(self): self.output_directory_path.mkdir( diff --git a/setup.py b/setup.py index d33f842..0fcd905 100644 --- a/setup.py +++ b/setup.py @@ -3,6 +3,6 @@ setuptools.setup( name = "fempy", - version = "0.2.0a0", + version = "0.3.0a0", packages = setuptools.find_packages(), ) diff --git a/tests/test__convection_diffusion.py b/tests/test__convection_diffusion.py index a90485c..46c8728 100755 --- a/tests/test__convection_diffusion.py +++ b/tests/test__convection_diffusion.py @@ -1,57 +1,58 @@ import firedrake as fe import fempy.models.convection_diffusion - -def test__verify_convergence_order_via_mms( - mesh_sizes = (16, 32), tolerance = 0.1, quadrature_degree = 2): - - class Model(fempy.models.convection_diffusion.Model): +class VerifiableModel(fempy.models.convection_diffusion.Model): - def __init__(self, meshsize = 4): - - self.meshsize = meshsize - - super().__init__() - - self.kinematic_viscosity.assign(0.1) - - def init_mesh(self): - - self.mesh = fe.UnitSquareMesh(self.meshsize, self.meshsize) + def __init__(self, quadrature_degree, spatial_order, meshsize): - def init_integration_measure(self): - - self.integration_measure = fe.dx(degree = 2) + self.meshsize = meshsize - def init_manufactured_solution(self): + super().__init__( + quadrature_degree = quadrature_degree, + spatial_order = spatial_order) - x = fe.SpatialCoordinate(self.mesh) - - sin, pi = fe.sin, fe.pi - - self.manufactured_solution = sin(2.*pi*x[0])*sin(pi*x[1]) - - ihat, jhat = self.unit_vectors() - - self.advection_velocity = sin(2.*pi*x[0])*sin(4.*pi*x[1])*ihat \ - + sin(pi*x[0])*sin(2.*pi*x[1])*jhat + self.kinematic_viscosity.assign(0.1) - def strong_form_residual(self, solution): - - x = fe.SpatialCoordinate(self.mesh) - - u = solution - - a = self.advection_velocity - - nu = self.kinematic_viscosity + def init_mesh(self): + + self.mesh = fe.UnitSquareMesh(self.meshsize, self.meshsize) + + def init_manufactured_solution(self): + + x = fe.SpatialCoordinate(self.mesh) + + sin, pi = fe.sin, fe.pi + + self.manufactured_solution = sin(2.*pi*x[0])*sin(pi*x[1]) + + ihat, jhat = self.unit_vectors() + + self.advection_velocity = sin(2.*pi*x[0])*sin(4.*pi*x[1])*ihat \ + + sin(pi*x[0])*sin(2.*pi*x[1])*jhat + + def strong_form_residual(self, solution): + + x = fe.SpatialCoordinate(self.mesh) + + u = solution + + a = self.advection_velocity + + nu = self.kinematic_viscosity + + dot, grad, div = fe.dot, fe.grad, fe.div + + return dot(a, grad(u)) - div(nu*grad(u)) - dot, grad, div = fe.dot, fe.grad, fe.div - return dot(a, grad(u)) - div(nu*grad(u)) - +def test__verify_convergence_order_via_mms( + mesh_sizes = (16, 32), tolerance = 0.1, quadrature_degree = 2): + fempy.mms.verify_spatial_order_of_accuracy( - Model = Model, + Model = VerifiableModel, + constructor_kwargs = { + "quadrature_degree": 2, + "spatial_order": 2}, expected_order = 2, mesh_sizes = mesh_sizes, tolerance = tolerance) diff --git a/tests/test__enthalpy.py b/tests/test__enthalpy.py index 5eb868e..0b946e9 100644 --- a/tests/test__enthalpy.py +++ b/tests/test__enthalpy.py @@ -5,22 +5,23 @@ class VerifiableModel(fempy.models.enthalpy.Model): - def __init__(self, meshsize): + def __init__(self, + quadrature_degree, + spatial_order, + temporal_order, + meshsize): self.meshsize = meshsize - super().__init__() - - self.update_initial_values() + super().__init__( + quadrature_degree = quadrature_degree, + spatial_order = spatial_order, + temporal_order = temporal_order) def init_mesh(self): self.mesh = fe.UnitIntervalMesh(self.meshsize) - def init_integration_measure(self): - - self.integration_measure = fe.dx(degree = 4) - def strong_form_residual(self, solution): T = solution @@ -44,23 +45,19 @@ def init_manufactured_solution(self): sin, pi, exp = fe.sin, fe.pi, fe.exp self.manufactured_solution = 0.5*sin(2.*pi*x)*(1. - 2*exp(-3.*t**2)) - - def update_initial_values(self): - - initial_values = fe.interpolate( - self.manufactured_solution, self.function_space) - - self.initial_values.assign(initial_values) - + + def test__verify_spatial_convergence_order_via_mms( mesh_sizes = (4, 8, 16, 32), timestep_size = 1./256., - tolerance = 0.1, - plot_errors = False, - plot_solution = False): + tolerance = 0.1): fempy.mms.verify_spatial_order_of_accuracy( Model = VerifiableModel, + constructor_kwargs = { + "quadrature_degree": 4, + "spatial_order": 2, + "temporal_order": 1}, parameters = { "stefan_number": 0.1, "smoothing": 1./32.}, @@ -69,74 +66,48 @@ def test__verify_spatial_convergence_order_via_mms( timestep_size = timestep_size, endtime = 1., tolerance = tolerance, - plot_errors = plot_errors, - plot_solution = plot_solution) + plot_errors = False, + plot_solution = False, + report = False) def test__verify_temporal_convergence_order_via_mms( meshsize = 256, timestep_sizes = (1./16., 1./32., 1./64., 1./128.), - tolerance = 0.1, - plot_errors = False, - plot_solution = False): + tolerance = 0.1): fempy.mms.verify_temporal_order_of_accuracy( Model = VerifiableModel, + constructor_kwargs = { + "quadrature_degree": 4, + "spatial_order": 2, + "temporal_order": 1}, + parameters = { + "stefan_number": 0.1, + "smoothing": 1./32.}, expected_order = 1, meshsize = meshsize, endtime = 1., timestep_sizes = timestep_sizes, tolerance = tolerance, - plot_errors = plot_errors, - plot_solution = plot_solution) - - -class SecondOrderVerifiableModel(VerifiableModel): + plot_errors = False, + plot_solution = False, + report = False) - def init_initial_values(self): - - self.initial_values = [fe.Function(self.function_space) - for i in range(2)] - - def update_initial_values(self): - - initial_values = fe.interpolate( - self.manufactured_solution, self.function_space) - - for iv in self.initial_values: - - iv.assign(initial_values) - - def init_time_discrete_terms(self): - """ Gear/BDF2 finite difference scheme - with constant time step size. """ - T_np1 = self.solution - - T_n = self.initial_values[0] - - T_nm1 = self.initial_values[1] - - Delta_t = self.timestep_size - - T_t = 1./Delta_t*(3./2.*T_np1 - 2.*T_n + 0.5*T_nm1) - - phil = self.porosity - - phil_t = 1./Delta_t*\ - (3./2.*phil(T_np1) - 2.*phil(T_n) + 0.5*phil(T_nm1)) - - self.time_discrete_terms = T_t, phil_t - -def test__verify_temporal_convergence_order_via_mms__bdf2( - meshsize = pow(2, 13), - timestep_sizes = (1./16., 1./32., 1./64., 1./128.), - tolerance = 0.1, +def test__verify_temporal_convergence__second_order__via_mms( + meshsize = 128, + timestep_sizes = (1./64., 1./128., 1./256.), + tolerance = 0.3, plot_errors = False, plot_solution = False): fempy.mms.verify_temporal_order_of_accuracy( - Model = SecondOrderVerifiableModel, + Model = VerifiableModel, + constructor_kwargs = { + "quadrature_degree": None, + "spatial_order": 3, + "temporal_order": 2}, parameters = { "stefan_number": 0.1, "smoothing": 1./32.}, diff --git a/tests/test__enthalpy_porosity.py b/tests/test__enthalpy_porosity.py deleted file mode 100644 index 3fb64d1..0000000 --- a/tests/test__enthalpy_porosity.py +++ /dev/null @@ -1,211 +0,0 @@ -import firedrake as fe -import fempy.mms -import fempy.models.enthalpy_porosity -import fempy.benchmarks.melting_octadecane - - -def test__melting_octadecane_benchmark__regression(): - - endtime, expected_liquid_area, tolerance = 30., 0.24, 0.01 - - nx = 32 - - Delta_t = 10. - - 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/" + "nx" + str(nx) + "_Deltat" + str(Delta_t) - + "_s" + str(s) + "_tf" + str(endtime) + "/") - - model.run(endtime = endtime, plot = False) - - p, u, T = model.solution.split() - - phil = model.porosity(T) - - liquid_area = fe.assemble(phil*fe.dx) - - print("Liquid area = " + str(liquid_area)) - - assert(abs(liquid_area - expected_liquid_area) < tolerance) - - phil_h = fe.interpolate(phil, model.function_space.sub(2)) - - max_phil = phil_h.vector().max() - - print("Maximum phil = " + str(max_phil)) - - assert(abs(max_phil - 1.) < tolerance) - - -def test__melting_octadecane_benchmark_with_darcy_resistance__regression(): - - endtime, expected_liquid_area, tolerance = 30., 0.24, 0.01 - - D = 1.e12 - - model = fempy.benchmarks.melting_octadecane.ModelWithDarcyResistance(meshsize = 32) - - model.timestep_size.assign(10.) - - model.smoothing.assign(1./256.) - - model.darcy_resistance_factor.assign(D) - - model.output_directory_path = model.output_directory_path.joinpath( - "melting_octadecane_with_darcy_resistance/D" + str(D) + "/") - - model.run(endtime = endtime, plot = False) - - p, u, T = model.solution.split() - - phil = model.porosity(T) - - liquid_area = fe.assemble(phil*fe.dx) - - print("Liquid area = " + str(liquid_area)) - - assert(abs(liquid_area - expected_liquid_area) < tolerance) - - phil_h = fe.interpolate(phil, model.function_space.sub(2)) - - max_phil = phil_h.vector().max() - - print("Maximum phil = " + str(max_phil)) - - assert(abs(max_phil - 1.) < tolerance) - - -class VerifiableModel(fempy.models.enthalpy_porosity.ModelWithBDF2): - - def __init__(self, meshsize): - - self.meshsize = meshsize - - super().__init__() - - def init_mesh(self): - - self.mesh = fe.UnitSquareMesh(self.meshsize, self.meshsize) - - def init_integration_measure(self): - - self.integration_measure = fe.dx(degree = 4) - - def strong_form_residual(self, solution): - - gamma = self.pressure_penalty_factor - - mu_s = self.solid_dynamic_viscosity - - mu_l = self.liquid_dynamic_viscosity - - Pr = self.prandtl_number - - Ste = self.stefan_number - - t = self.time - - grad, dot, div, sym, diff = fe.grad, fe.dot, fe.div, fe.sym, fe.diff - - p, u, T = solution - - b = self.buoyancy(T) - - phil = self.porosity(T) - - mu = mu_s + (mu_l - mu_s)*phil - - 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 + 1./Ste*phil, t) + div(T*u) - 1./Pr*div(grad(T)) - - return r_p, r_u, r_T - - def init_manufactured_solution(self): - - pi, sin, cos, exp = fe.pi, fe.sin, fe.cos, fe.exp - - x = fe.SpatialCoordinate(self.mesh) - - t = self.time - - t_f = fe.Constant(1.) - - ihat, jhat = self.unit_vectors() - - 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 = -sin(pi*x[0])*sin(2.*pi*x[1]) - - T = 0.5*sin(2.*pi*x[0])*sin(pi*x[1])*(1. - exp(-t**2)) - - self.manufactured_solution = p, u, T - - def update_initial_values(self): - - for u_m, V in zip( - self.manufactured_solution, self.function_space): - - self.initial_values.assign(fe.interpolate(u_m, V)) - - def solve(self): - - self.solver.parameters["snes_monitor"] = False - - super().solve() - - print("Solved at time t = " + str(self.time.__float__())) - - -def test__verify_spatial_convergence_order_via_mms( - parameters = { - "grashof_number": 2., - "prandtl_number": 5., - "stefan_number": 0.2, - "smoothing": 1./16.}, - mesh_sizes = (8, 16, 32), - timestep_size = 1./256., - tolerance = 0.4): - - fempy.mms.verify_spatial_order_of_accuracy( - Model = VerifiableModel, - parameters = parameters, - expected_order = 2, - mesh_sizes = mesh_sizes, - tolerance = tolerance, - timestep_size = timestep_size, - 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) - \ No newline at end of file diff --git a/tests/test__heat.py b/tests/test__heat.py index da080e9..2294523 100755 --- a/tests/test__heat.py +++ b/tests/test__heat.py @@ -2,15 +2,20 @@ import fempy.models.heat -class Model(fempy.models.heat.Model): +class VerifiableModel(fempy.models.heat.Model): - def __init__(self, meshsize): + def __init__(self, + quadrature_degree, + spatial_order, + temporal_order, + meshsize): self.meshsize = meshsize - super().__init__() - - self.update_initial_values() + super().__init__( + quadrature_degree = quadrature_degree, + spatial_order = spatial_order, + temporal_order = temporal_order) def init_mesh(self): @@ -35,126 +40,82 @@ def init_manufactured_solution(self): sin, pi, exp = fe.sin, fe.pi, fe.exp self.manufactured_solution = sin(2.*pi*x)*exp(-pow(t, 2)) - - def update_initial_values(self): - initial_values = fe.interpolate( - self.manufactured_solution, self.function_space) - - self.initial_values.assign(initial_values) - def init_solver(self, solver_parameters = {"ksp_type": "cg"}): self.solver = fe.NonlinearVariationalSolver( self.problem, solver_parameters = solver_parameters) -def test__verify_spatial_convergence_order_via_mms( +def test__verify_spatial_convergence__second_order__via_mms( mesh_sizes = (4, 8, 16, 32), timestep_size = 1./64., tolerance = 0.1): fempy.mms.verify_spatial_order_of_accuracy( - Model = Model, + Model = VerifiableModel, + constructor_kwargs = { + "quadrature_degree": None, "spatial_order": 2, "temporal_order": 1}, expected_order = 2, mesh_sizes = mesh_sizes, tolerance = tolerance, timestep_size = timestep_size, - endtime = 1.) + endtime = 1., + plot_solution = False, + report = False) -def test__verify_temporal_convergence_order_via_mms( +def test__verify_temporal_convergence__first_order__via_mms( meshsize = 256, timestep_sizes = (1./4., 1./8., 1./16., 1./32.), tolerance = 0.1): fempy.mms.verify_temporal_order_of_accuracy( - Model = Model, + Model = VerifiableModel, + constructor_kwargs = { + "quadrature_degree": None, "spatial_order": 2, "temporal_order": 1}, expected_order = 1, meshsize = meshsize, endtime = 1., timestep_sizes = timestep_sizes, - tolerance = tolerance) - - -class SecondOrderModel(Model): - - def init_initial_values(self): - - self.initial_values = [fe.Function(self.function_space) - for i in range(2)] - - def update_initial_values(self): - - initial_values = fe.interpolate( - self.manufactured_solution, self.function_space) - - for iv in self.initial_values: - - iv.assign(initial_values) - - def init_time_discrete_terms(self): - """ Gear/BDF2 finite difference scheme - with constant time step size. """ - unp1 = self.solution - - un = self.initial_values[0] - - unm1 = self.initial_values[1] - - Delta_t = self.timestep_size - - self.time_discrete_terms = 1./Delta_t*(3./2.*unp1 - 2.*un + 0.5*unm1) - - def init_manufactured_solution(self): - - x = fe.SpatialCoordinate(self.mesh)[0] - - t = self.time - - sin, pi, exp = fe.sin, fe.pi, fe.exp - - self.manufactured_solution = sin(2.*pi*x)*exp(-pow(t, 3)) + tolerance = tolerance, + plot_solution = False, + report = False) -def test__verify_temporal_convergence_order_via_mms__bdf2( - meshsize = 256, - timestep_sizes = (1./2., 1./4., 1./8., 1./16.), +def test__verify_temporal_convergence__second_order__via_mms( + meshsize = 128, + timestep_sizes = ( + 1./4., 1./8., 1./16., 1./32., 1./64., 1./128.), tolerance = 0.1): fempy.mms.verify_temporal_order_of_accuracy( - Model = SecondOrderModel, + Model = VerifiableModel, + constructor_kwargs = { + "quadrature_degree": None, "spatial_order": 3, "temporal_order": 2}, expected_order = 2, meshsize = meshsize, endtime = 1., timestep_sizes = timestep_sizes, tolerance = tolerance, - plot_solution = True) - - -class ModelWithWave(Model): - - def init_manufactured_solution(self): - - x = fe.SpatialCoordinate(self.mesh)[0] - - t = self.time - - sin, pi, exp = fe.sin, fe.pi, fe.exp + plot_solution = False, + report = False) - self.manufactured_solution = 0.5*sin(2.*pi*x - pi/4.*(2.*t + 1.)) - -def __fails__test_verify_spatial_convergence_order_via_mms_with_wave_solution( - mesh_sizes = (4, 8, 16, 32), - timestep_size = 1./64., +def test__verify_temporal_convergence__third_order__via_mms( + meshsize = 128, + timestep_sizes = (1./4., 1./8., 1./16., 1./32.), tolerance = 0.1): - fempy.mms.verify_spatial_order_of_accuracy( - Model = ModelWithWave, - expected_order = 2, - mesh_sizes = mesh_sizes, + fempy.mms.verify_temporal_order_of_accuracy( + Model = VerifiableModel, + constructor_kwargs = { + "quadrature_degree": None, "spatial_order": 3, "temporal_order": 3}, + expected_order = 3, + meshsize = meshsize, + endtime = 1., + timestep_sizes = timestep_sizes, tolerance = tolerance, - timestep_size = timestep_size, - endtime = 1.) + plot_solution = False, + report = False) \ No newline at end of file diff --git a/tests/test__laplace.py b/tests/test__laplace.py index 7ee8e58..7308aac 100755 --- a/tests/test__laplace.py +++ b/tests/test__laplace.py @@ -15,13 +15,13 @@ def strong_form_residual(self, solution): class OneDimMMSModel(VerifiableModel): - def __init__(self, meshsize): + def __init__(self, quadrature_degree, spatial_order, meshsize): self.meshsize = meshsize - super().__init__() - - self.integration_measure = fe.dx(degree = 2) + super().__init__( + quadrature_degree = quadrature_degree, + spatial_order = spatial_order) def init_mesh(self): @@ -37,9 +37,10 @@ def init_manufactured_solution(self): def test__verify_convergence_order_via_mms( - mesh_sizes = (16, 32), tolerance = 0.1): + mesh_sizes = (8, 16, 32), tolerance = 0.1): fempy.mms.verify_spatial_order_of_accuracy( + constructor_kwargs = {"quadrature_degree": 2, "spatial_order": 2}, Model = OneDimMMSModel, expected_order = 2, mesh_sizes = mesh_sizes, diff --git a/tests/test__navier_stokes.py b/tests/test__navier_stokes.py index 6f36ef1..7173530 100755 --- a/tests/test__navier_stokes.py +++ b/tests/test__navier_stokes.py @@ -2,25 +2,20 @@ import fempy.models.navier_stokes -def test__verify_convergence_order_via_MMS( - mesh_sizes = (16, 32), tolerance = 0.1): - - class Model(fempy.models.navier_stokes.Model): +class VerifiableModel(fempy.models.navier_stokes.Model): - def __init__(self, meshsize): + def __init__(self, quadrature_degree, spatial_order, meshsize): self.meshsize = meshsize - super().__init__() + super().__init__( + quadrature_degree = quadrature_degree, + spatial_order = spatial_order) def init_mesh(self): self.mesh = fe.UnitSquareMesh(self.meshsize, self.meshsize) - def init_integration_measure(self): - - self.integration_measure = fe.dx(degree = 4) - def init_manufactured_solution(self): sin, pi = fe.sin, fe.pi @@ -47,9 +42,14 @@ def strong_form_residual(self, solution): r_p = div(u) return r_u, r_p - + + +def test__verify_convergence_order_via_mms( + mesh_sizes = (16, 32), tolerance = 0.1): + fempy.mms.verify_spatial_order_of_accuracy( - Model = Model, + Model = VerifiableModel, + constructor_kwargs = {"quadrature_degree": 4, "spatial_order": 2}, expected_order = 2, mesh_sizes = mesh_sizes, tolerance = tolerance) diff --git a/tests/test__navier_stokes_boussinesq.py b/tests/test__navier_stokes_boussinesq.py index b66734f..c5d7e00 100755 --- a/tests/test__navier_stokes_boussinesq.py +++ b/tests/test__navier_stokes_boussinesq.py @@ -6,20 +6,18 @@ class VerifiableModel(fempy.models.navier_stokes_boussinesq.Model): - def __init__(self, meshsize): + def __init__(self, quadrature_degree, spatial_order, meshsize): self.meshsize = meshsize - super().__init__() + super().__init__( + quadrature_degree = quadrature_degree, + spatial_order = spatial_order) def init_mesh(self): self.mesh = fe.UnitSquareMesh(self.meshsize, self.meshsize) - def init_integration_measure(self): - - self.integration_measure = fe.dx(degree = 4) - def init_manufactured_solution(self): sin, pi = fe.sin, fe.pi @@ -42,8 +40,6 @@ def init_manufactured_solution(self): def strong_form_residual(self, solution): - mu = self.dynamic_viscosity - Gr = self.grashof_number Pr = self.prandtl_number @@ -56,7 +52,7 @@ def strong_form_residual(self, solution): r_p = div(u) - r_u = grad(u)*u + grad(p) - 2.*div(mu*sym(grad(u))) + Gr*T*ghat + r_u = grad(u)*u + grad(p) - 2.*div(sym(grad(u))) + Gr*T*ghat r_T = dot(u, grad(T)) - 1./Pr*div(grad(T)) @@ -75,8 +71,9 @@ def test__verify_convergence_order_via_mms( fempy.mms.verify_spatial_order_of_accuracy( Model = VerifiableModel, + constructor_kwargs = { + "quadrature_degree": 4, "spatial_order": 2}, parameters = { - "dynamic_viscosity": 0.1, "grashof_number": Ra/Pr, "prandtl_number": Pr}, expected_order = 2, @@ -156,7 +153,8 @@ def unsteadiness(model): def test__verify_against_heat_driven_cavity_benchmark(): - model = fempy.benchmarks.heat_driven_cavity.Model(meshsize = 40) + model = fempy.benchmarks.heat_driven_cavity.Model( + quadrature_degree = 4, spatial_order = 2, meshsize = 40) model.solver.solve() diff --git a/tests/test__unsteady_navier_stokes.py b/tests/test__unsteady_navier_stokes.py index f533d6e..52adae5 100644 --- a/tests/test__unsteady_navier_stokes.py +++ b/tests/test__unsteady_navier_stokes.py @@ -4,22 +4,20 @@ class Model(fempy.models.unsteady_navier_stokes.Model): - def __init__(self, meshsize): + def __init__(self, + quadrature_degree, spatial_order, temporal_order, meshsize): self.meshsize = meshsize - super().__init__() - - self.update_initial_values() + super().__init__( + quadrature_degree = quadrature_degree, + spatial_order = spatial_order, + temporal_order = temporal_order) def init_mesh(self): self.mesh = fe.UnitSquareMesh(self.meshsize, self.meshsize) - def init_integration_measure(self): - - self.integration_measure = fe.dx(degree = 4) - def init_manufactured_solution(self): exp, sin, pi = fe.exp, fe.sin, fe.pi @@ -50,22 +48,17 @@ def strong_form_residual(self, solution): r_p = div(u) return r_u, r_p - - def update_initial_values(self): - - for u_m, V in zip( - self.manufactured_solution, self.function_space): - - self.initial_values.assign(fe.interpolate(u_m, V)) def test__verify_spatial_convergence_order_via_mms( - mesh_sizes = (4, 8, 16, 32), - timestep_size = 1./64., - tolerance = 0.2): + mesh_sizes = (3, 6, 12, 24), + timestep_size = 1./32., + tolerance = 0.3): fempy.mms.verify_spatial_order_of_accuracy( Model = Model, + constructor_kwargs = { + "quadrature_degree": 4, "spatial_order": 2, "temporal_order": 1}, expected_order = 2, mesh_sizes = mesh_sizes, tolerance = tolerance, @@ -75,11 +68,13 @@ def test__verify_spatial_convergence_order_via_mms( def test__verify_temporal_convergence_order_via_mms( meshsize = 32, - timestep_sizes = (1., 1./2., 1./4., 1./8.), + timestep_sizes = (1./2., 1./4., 1./8.), tolerance = 0.1): fempy.mms.verify_temporal_order_of_accuracy( Model = Model, + constructor_kwargs = { + "quadrature_degree": 4, "spatial_order": 2, "temporal_order": 1}, expected_order = 1, meshsize = meshsize, endtime = 1., diff --git a/tests/test__validate__enthalpy_porosity.py b/tests/test__validate__enthalpy_porosity.py new file mode 100644 index 0000000..06e2a2e --- /dev/null +++ b/tests/test__validate__enthalpy_porosity.py @@ -0,0 +1,158 @@ +import firedrake as fe +import fempy.mms +import fempy.benchmarks.melt_octadecane +import fempy.benchmarks.freeze_water + + +def test__regression__validate__melt_octadecane(): + + endtime = 80. + + topwall_heatflux_switchtime = 40. + + q = -0.02 + + + s = 1./64. + + nx = 32 + + Delta_t = 10. + + tau = 1.e-12 + + + expected_liquid_area = 0.61 + + tolerance = 0.01 + + + model = fempy.benchmarks.melt_octadecane.Model( + quadrature_degree = 4, + spatial_order = 2, + temporal_order = 2, + meshsize = nx) + + model.timestep_size.assign(Delta_t) + + model.smoothing.assign(s) + + model.solid_velocity_relaxation_factor.assign(tau) + + model.topwall_heatflux_switchtime = topwall_heatflux_switchtime + \ + 2.*model.time_tolerance + + model.topwall_heatflux_postswitch = q + + model.save_smoothing_sequence = True + + model.output_directory_path = model.output_directory_path.joinpath( + "melt_octadecane/with_heatflux/switchtime" + + str(topwall_heatflux_switchtime) + + "_tf" + str(endtime) + "/" + + "q" + str(q) + "/" + + "second_order_" + + "nx" + str(nx) + "_Deltat" + str(Delta_t) + + "_s" + str(s) + "_tau" + str(tau) + "/" + "tf" + str(endtime) + "/") + + model.run(endtime = endtime, plot = True, report = True) + + + p, u, T = model.solution.split() + + phil = model.porosity(T) + + liquid_area = fe.assemble(phil*fe.dx) + + print("Liquid area = " + str(liquid_area)) + + assert(abs(liquid_area - expected_liquid_area) < tolerance) + + +def test__regression__validate__freeze_water(): + + mu_l__SI = 8.90e-4 # [Pa s] + + rho_l__SI = 999.84 # [kg / m^3] + + nu_l__SI = mu_l__SI/rho_l__SI # [m^2 / s] + + t_f__SI = 2340. # [s] + + L__SI = 0.038 # [m] + + Tau = pow(L__SI, 2)/nu_l__SI + + t_f = t_f__SI/Tau + + """ For Kowalewski's water freezing experiment, + at t_f__SI 2340 s, t_f = 1.44. + """ + plot = False + + write_solution = False + + report = False + + spatial_dimensions = 2 + + + s = 1./200. + + tau = 1.e-12 + + + rx = 2 + + nx = 32 + + + rt = 2 + + nt = 4 + + + q = 4 + + + expected_liquid_area = 0.69 + + tolerance = 0.01 + + + model = fempy.benchmarks.freeze_water.Model( + quadrature_degree = q, + spatial_order = rx, + temporal_order = rt, + spatial_dimensions = spatial_dimensions, + meshsize = nx) + + model.timestep_size.assign(t_f/float(nt)) + + model.solid_velocity_relaxation_factor.assign(tau) + + model.smoothing.assign(s) + + model.save_smoothing_sequence = False + + model.output_directory_path = model.output_directory_path.joinpath( + "freeze_water/" + + "s{0}_tau{1}/rx{2}_nx{3}_rt{4}_nt{5}/q{6}/tf{7}/dim{8}/".format( + s, tau, rx, nx, rt, nt, q, t_f, spatial_dimensions)) + + model.run( + endtime = t_f, + write_solution = write_solution, + plot = plot, + report = report) + + p, u, T = model.solution.split() + + phil = model.porosity(T) + + liquid_area = fe.assemble(phil*fe.dx) + + print("Liquid area = " + str(liquid_area)) + + assert(abs(liquid_area - expected_liquid_area) < tolerance) + \ No newline at end of file diff --git a/tests/test__verify__enthalpy_porosity.py b/tests/test__verify__enthalpy_porosity.py new file mode 100644 index 0000000..346f8e1 --- /dev/null +++ b/tests/test__verify__enthalpy_porosity.py @@ -0,0 +1,192 @@ +import firedrake as fe +import fempy.mms +import fempy.models.enthalpy_porosity + + +class VerifiableModel(fempy.models.enthalpy_porosity.Model): + + def __init__(self, *args, meshsize, **kwargs): + + self.meshsize = meshsize + + super().__init__(*args, **kwargs) + + def init_mesh(self): + + self.mesh = fe.UnitSquareMesh(self.meshsize, self.meshsize) + + def strong_form_residual(self, solution): + + Pr = self.prandtl_number + + Ste = self.stefan_number + + t = self.time + + grad, dot, div, sym, diff = fe.grad, fe.dot, fe.div, fe.sym, fe.diff + + p, u, T = solution + + b = self.buoyancy(T) + + d = self.solid_velocity_relaxation(T) + + phil = self.porosity(T) + + cp = self.phase_dependent_material_property( + self.heat_capacity_solid_to_liquid_ratio)(phil) + + k = self.phase_dependent_material_property( + self.thermal_conductivity_solid_to_liquid_ratio)(phil) + + r_p = div(u) + + r_u = diff(u, t) + grad(u)*u + grad(p) - 2.*div(sym(grad(u))) \ + + b + d*u + + r_T = diff(cp*T, t) + dot(u, grad(cp*T)) \ + - 1./Pr*div(k*grad(T)) + 1./Ste*diff(cp*phil, t) + + return r_p, r_u, r_T + + def init_manufactured_solution(self): + + pi, sin, cos, exp = fe.pi, fe.sin, fe.cos, fe.exp + + x = fe.SpatialCoordinate(self.mesh) + + t = self.time + + t_f = fe.Constant(1.) + + ihat, jhat = self.unit_vectors() + + 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 = -sin(pi*x[0])*sin(2.*pi*x[1]) + + T = 0.5*sin(2.*pi*x[0])*sin(pi*x[1])*(1. - exp(-t**2)) + + self.manufactured_solution = p, u, T + + +def test__verify__second_order_spatial_convergence__via_mms( + constructor_kwargs = { + "quadrature_degree": 4, + "spatial_order": 2, + "temporal_order": 2}, + parameters = { + "grashof_number": 2., + "prandtl_number": 5., + "stefan_number": 0.2, + "heat_capacity_solid_to_liquid_ratio": 0.500, + "thermal_conductivity_solid_to_liquid_ratio": 2.14/0.561, + "smoothing": 1./16.}, + mesh_sizes = (2, 4, 8), + timestep_size = 1./128., + tolerance = 0.02): + + fempy.mms.verify_spatial_order_of_accuracy( + Model = VerifiableModel, + constructor_kwargs = constructor_kwargs, + parameters = parameters, + expected_order = 2, + mesh_sizes = mesh_sizes, + tolerance = tolerance, + timestep_size = timestep_size, + endtime = 0.5, + plot_solution = False, + plot_errors = False, + report = False) + + +def test__verify__second_order_temporal_convergence__via_mms( + constructor_kwargs = { + "quadrature_degree": 4, + "spatial_order": 2, + "temporal_order": 2}, + parameters = { + "grashof_number": 2., + "prandtl_number": 5., + "stefan_number": 0.2, + "heat_capacity_solid_to_liquid_ratio": 0.500, + "thermal_conductivity_solid_to_liquid_ratio": 2.14/0.561, + "smoothing": 1./16.}, + meshsize = 24, + timestep_sizes = (1./8., 1./16., 1./32.), + tolerance = 0.4): + + fempy.mms.verify_temporal_order_of_accuracy( + Model = VerifiableModel, + constructor_kwargs = constructor_kwargs, + parameters = parameters, + expected_order = 2, + meshsize = meshsize, + tolerance = tolerance, + timestep_sizes = timestep_sizes, + endtime = 0.5, + plot_solution = False, + plot_errors = False, + report = False) + + +def test__verify__third_order_spatial_convergence__via_mms( + constructor_kwargs = { + "quadrature_degree": 2, + "spatial_order": 3, + "temporal_order": 3}, + parameters = { + "grashof_number": 2., + "prandtl_number": 5., + "stefan_number": 0.2, + "heat_capacity_solid_to_liquid_ratio": 0.500, + "thermal_conductivity_solid_to_liquid_ratio": 2.14/0.561, + "smoothing": 1./16.}, + mesh_sizes = (3, 6, 12, 24), + timestep_size = 1./64., + tolerance = 0.32): + + fempy.mms.verify_spatial_order_of_accuracy( + Model = VerifiableModel, + constructor_kwargs = constructor_kwargs, + parameters = parameters, + expected_order = 3, + mesh_sizes = mesh_sizes, + tolerance = tolerance, + timestep_size = timestep_size, + endtime = 0.32, + plot_solution = False, + plot_errors = False, + report = False) + + +def test__verify__third_order_temporal_convergence__via_mms( + constructor_kwargs = { + "quadrature_degree": 8, + "spatial_order": 3, + "temporal_order": 3}, + parameters = { + "grashof_number": 2., + "prandtl_number": 5., + "stefan_number": 0.2, + "heat_capacity_solid_to_liquid_ratio": 0.500, + "thermal_conductivity_solid_to_liquid_ratio": 2.14/0.561, + "smoothing": 1./16.}, + meshsize = 16, + timestep_sizes = (1./2., 1./4., 1./8., 1./16.), + tolerance = 0.15): + + fempy.mms.verify_temporal_order_of_accuracy( + Model = VerifiableModel, + constructor_kwargs = constructor_kwargs, + parameters = parameters, + expected_order = 3, + meshsize = meshsize, + tolerance = tolerance, + timestep_sizes = timestep_sizes, + endtime = 0.5, + plot_solution = False, + plot_errors = False, + report = False) + \ No newline at end of file