From 5148078c144a25cd60ae83ab91618abf145cc67c Mon Sep 17 00:00:00 2001 From: "Alexander G. Zimmerman" Date: Tue, 21 Jan 2020 10:59:26 +0100 Subject: [PATCH] Simplify the Simulation interface, improve MMS verification, and clean 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 | 51 +++-- sapphire/benchmarks/heat_driven_cavity.py | 32 +-- .../benchmarks/melt_octadecane_in_cavity.py | 73 ++----- sapphire/mms.py | 169 +++++++++------- sapphire/simulation.py | 41 ++-- .../convection_coupled_phasechange.py | 129 +++++++----- sapphire/simulations/convection_diffusion.py | 7 +- sapphire/simulations/heat.py | 13 +- sapphire/simulations/laplace.py | 11 +- .../simulations/navier_stokes_boussinesq.py | 33 ++- sapphire/simulations/phasechange.py | 23 ++- setup.py | 2 +- .../test__convection_coupled_phasechange.py | 127 +++--------- .../test__convection_coupled_phasechange.py | 190 ++++++++++++------ .../test__convection_diffusion.py | 16 +- tests/model_verification/test__heat.py | 70 +++---- tests/model_verification/test__laplace.py | 3 +- .../model_verification/test__navier_stokes.py | 3 +- .../test__navier_stokes_boussinesq.py | 63 +++++- tests/model_verification/test__phasechange.py | 80 +++++--- .../test__unsteady_navier_stokes.py | 47 +++-- 21 files changed, 652 insertions(+), 531 deletions(-) diff --git a/sapphire/benchmarks/freeze_water_in_cavity.py b/sapphire/benchmarks/freeze_water_in_cavity.py index 787b5c1..745874c 100644 --- a/sapphire/benchmarks/freeze_water_in_cavity.py +++ b/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,19 +151,30 @@ 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, @@ -170,20 +182,15 @@ def __init__(self, *args, meshsize, **kwargs): 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) \ No newline at end of file diff --git a/sapphire/benchmarks/heat_driven_cavity.py b/sapphire/benchmarks/heat_driven_cavity.py index 746cc66..56a61c8 100644 --- a/sapphire/benchmarks/heat_driven_cavity.py +++ b/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) \ No newline at end of file diff --git a/sapphire/benchmarks/melt_octadecane_in_cavity.py b/sapphire/benchmarks/melt_octadecane_in_cavity.py index a28b366..5d86e16 100644 --- a/sapphire/benchmarks/melt_octadecane_in_cavity.py +++ b/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 \ No newline at end of file diff --git a/sapphire/mms.py b/sapphire/mms.py index e2e8fab..3da8ccc 100644 --- a/sapphire/mms.py +++ b/sapphire/mms.py @@ -59,8 +59,8 @@ def mms_initial_values(sim, manufactured_solution): return initial_values -def mms_dirichlet_boundary_conditions(sim, manufactured_solution): - +def default_mms_dirichlet_boundary_conditions(sim, manufactured_solution): + """ By default, apply Dirichlet BC's to every component on every boundary. """ W = sim.function_space if type(sim.element) is fe.FiniteElement: @@ -76,13 +76,23 @@ def mms_dirichlet_boundary_conditions(sim, manufactured_solution): def make_mms_verification_sim_class( sim_module, - manufactured_solution): + manufactured_solution, + strong_residual = None, + mms_dirichlet_boundary_conditions = None): + if strong_residual is None: + + strong_residual = sim_module.strong_residual + def initial_values(sim): return mms_initial_values( sim = sim, manufactured_solution = manufactured_solution(sim)) + + if mms_dirichlet_boundary_conditions is None: + + mms_dirichlet_boundary_conditions = default_mms_dirichlet_boundary_conditions def dirichlet_boundary_conditions(sim): @@ -101,7 +111,7 @@ def __init__(self, *args, **kwargs): self.variational_form_residual -= mms_source( sim = self, - strong_residual = sim_module.strong_residual, + strong_residual = strong_residual, manufactured_solution = manufactured_solution)\ *fe.dx(degree = self.quadrature_degree) @@ -112,41 +122,26 @@ def write_outputs(self, *args, **kwargs): return MMSVerificationSimulation -def errornorm(w, wh, *args, **kwargs): - """ Extends fe.errornorm to handle mixed FEM functions """ - if len(wh.split()) == 1: - - return fe.errornorm(w, wh.split()[0], *args, **kwargs) - - else: - - e = 0. - - for w_i, wh_i in zip(w, wh.split()): - - e += fe.errornorm(w_i, wh_i, *args, **kwargs) - - return e - - def verify_spatial_order_of_accuracy( sim_module, manufactured_solution, meshes, - expected_order, + norms, + expected_orders, tolerance, - sim_constructor_kwargs = {}, - parameters = {}, - timestep_size = 1.e32, + sim_parameters = {}, endtime = 0., - starttime = 0., + strong_residual = None, + dirichlet_boundary_conditions = None, outfile = None): MMSVerificationSimulation = make_mms_verification_sim_class( sim_module = sim_module, - manufactured_solution = manufactured_solution) + manufactured_solution = manufactured_solution, + strong_residual = strong_residual, + mms_dirichlet_boundary_conditions = dirichlet_boundary_conditions) - table = sapphire.table.Table(("h", "L2_error", "spatial_order")) + table = sapphire.table.Table(("h", "cellcount", "dofcount", "errors", "spatial_orders")) print("") @@ -154,45 +149,55 @@ def verify_spatial_order_of_accuracy( h = mesh.cell_sizes((0.,)*mesh.geometric_dimension()) - sim = MMSVerificationSimulation(mesh = mesh, **sim_constructor_kwargs) - - sim = sim.assign_parameters(parameters) + sim = MMSVerificationSimulation(mesh = mesh, **sim_parameters) if sim.time is not None: - sim.time = sim.time.assign(starttime) - - sim.timestep_size = sim.timestep_size.assign(timestep_size) - sim.solutions, _ = sim.run(endtime = endtime) else: sim.solution = sim.solve() - table.append({ - "h": h, - "L2_error": - errornorm( - manufactured_solution(sim), - sim.solution, - norm_type = "L2")}) + errors = [] + + w = manufactured_solution(sim) + + wh = sim.solution + + if type(w) is not type((0,)): + + w = (w,) + + for w_i, wh_i, norm in zip(w, wh.split(), norms): + + errors.append(fe.errornorm(w_i, wh_i, norm_type = norm)) + + cellcount = mesh.topology.num_cells() + + dofcount = len(wh.vector().array()) + + table.append({"h": h, "cellcount": cellcount, "dofcount": dofcount, "errors": errors}) if len(table) > 1: - h, e = table.data["h"], table.data["L2_error"] + h, e = table.data["h"], table.data["errors"] log = math.log - r = h[-2]/h[-1] + orders = [] + + for i in range(len(sim.solution.split())): + + r = h[-2]/h[-1] - order = log(e[-2]/e[-1])/log(r) + orders.append(log(e[-2][i]/e[-1][i])/log(r)) - table.data["spatial_order"][-1] = order + table.data["spatial_orders"][-1] = orders print(str(table)) - print("Last observed spatial order of accuracy is {0}".format(order)) + print("Last observed spatial orders of accuracy are {}".format(orders)) if outfile: @@ -200,38 +205,43 @@ def verify_spatial_order_of_accuracy( outfile.write(str(table)) - assert(abs(order - expected_order) < tolerance) + for order, expected_order in zip(orders, expected_orders): + + if expected_order is not None: + + assert(abs(order - expected_order) < tolerance) def verify_temporal_order_of_accuracy( sim_module, manufactured_solution, - mesh, timestep_sizes, endtime, - expected_order, + norms, + expected_orders, tolerance, - sim_constructor_kwargs = {}, - parameters = {}, + sim_parameters = {}, starttime = 0., + strong_residual = None, + dirichlet_boundary_conditions = None, outfile = None): MMSVerificationSimulation = make_mms_verification_sim_class( sim_module = sim_module, - manufactured_solution = manufactured_solution) + manufactured_solution = manufactured_solution, + strong_residual = strong_residual, + mms_dirichlet_boundary_conditions = dirichlet_boundary_conditions) - table = sapphire.table.Table(("Delta_t", "L2_error", "temporal_order")) + table = sapphire.table.Table(("Delta_t", "errors", "temporal_orders")) print("") for timestep_size in timestep_sizes: - sim = MMSVerificationSimulation(mesh = mesh, **sim_constructor_kwargs) + sim = MMSVerificationSimulation(**sim_parameters) assert(len(sim.solutions) > 1) - sim = sim.assign_parameters(parameters) - sim.timestep_size = sim.timestep_size.assign(timestep_size) sim.time = sim.time.assign(starttime) @@ -242,28 +252,41 @@ def verify_temporal_order_of_accuracy( sim.solutions, _, = sim.run(endtime = endtime) - table.append({ - "Delta_t": timestep_size, - "L2_error": errornorm( - manufactured_solution(sim), - sim.solution, - norm_type = "L2")}) + errors = [] + + w = manufactured_solution(sim) + + wh = sim.solution + + if type(w) is not type((0,)): + + w = (w,) + + for w_i, wh_i, norm in zip(w, wh.split(), norms): + + errors.append(fe.errornorm(w_i, wh_i, norm_type = norm)) + table.append({"Delta_t": timestep_size, "errors": errors}) + if len(table) > 1: - Delta_t, e = table.data["Delta_t"], table.data["L2_error"] + Delta_t, e = table.data["Delta_t"], table.data["errors"] log = math.log - r = Delta_t[-2]/Delta_t[-1] - - order = log(e[-2]/e[-1])/log(r) - - table.data["temporal_order"][-1] = order + orders = [] + + for i in range(len(sim.solution.split())): + + r = Delta_t[-2]/Delta_t[-1] + + orders.append(log(e[-2][i]/e[-1][i])/log(r)) + + table.data["temporal_orders"][-1] = orders print(str(table)) - print("Last observed temporal order of accuracy is {0}".format(order)) + print("Last observed temporal orders of accuracy are {}".format(orders)) if outfile: @@ -271,5 +294,9 @@ def verify_temporal_order_of_accuracy( outfile.write(str(table)) - assert(abs(order - expected_order) < tolerance) + for order, expected_order in zip(orders, expected_orders): + + if expected_order is not None: + + assert(abs(order - expected_order) < tolerance) \ No newline at end of file diff --git a/sapphire/simulation.py b/sapphire/simulation.py index a6ca762..cbe693e 100644 --- a/sapphire/simulation.py +++ b/sapphire/simulation.py @@ -17,7 +17,15 @@ def __init__(self, initial_values, quadrature_degree = None, time_dependent = True, + timestep_size = 1., time_stencil_size = 2, + solver_parameters = { + "snes_type": "newtonls", + "snes_monitor": None, + "ksp_type": "preonly", + "pc_type": "lu", + "mat_type": "aij", + "pc_factor_mat_solver_type": "mumps"}, output_directory_path = "output/"): self.mesh = mesh @@ -43,7 +51,7 @@ def __init__(self, self.time = fe.Constant(0.) - self.timestep_size = fe.Constant(1.) + self.timestep_size = fe.Constant(timestep_size) else: @@ -51,6 +59,7 @@ def __init__(self, self.timestep_size = None + self.solver_parameters = solver_parameters self.output_directory_path = pathlib.Path(output_directory_path) @@ -73,17 +82,9 @@ def __init__(self, self.dirichlet_boundary_conditions = \ dirichlet_boundary_conditions(sim = self) - self.snes_iteration_count = 0 - def solve(self, - parameters = { - "snes_type": "newtonls", - "snes_monitor": None, - "ksp_type": "preonly", - "pc_type": "lu", - "mat_type": "aij", - "pc_factor_mat_solver_type": "mumps"}): + def solve(self): problem = fe.NonlinearVariationalProblem( F = self.variational_form_residual, @@ -93,7 +94,7 @@ def solve(self, solver = fe.NonlinearVariationalSolver( problem = problem, - solver_parameters = parameters) + solver_parameters = self.solver_parameters) solver.solve() @@ -149,7 +150,7 @@ def run(self, while self.time.__float__() < (endtime - time_tolerance): - self.time.assign(self.time + self.timestep_size) + self.time = self.time.assign(self.time + self.timestep_size) self.solution = solve() @@ -161,22 +162,6 @@ def run(self, return self.solutions, self.time - def assign_parameters(self, parameters): - - for key, value in parameters.items(): - - attribute = getattr(self, key) - - if type(attribute) is type(fe.Constant(0.)): - - attribute.assign(value) - - else: - - setattr(self, key, value) - - return self - def unit_vectors(self): return unit_vectors(self.mesh) diff --git a/sapphire/simulations/convection_coupled_phasechange.py b/sapphire/simulations/convection_coupled_phasechange.py index 1a5061e..36d99f4 100644 --- a/sapphire/simulations/convection_coupled_phasechange.py +++ b/sapphire/simulations/convection_coupled_phasechange.py @@ -5,12 +5,19 @@ def element(cell, degree): - - scalar = fe.FiniteElement("P", cell, degree) - vector = fe.VectorElement("P", cell, degree) + if type(degree) is type(1): + + degree = (degree,)*3 + + pressure_element = fe.FiniteElement("P", cell, degree[0]) + + velocity_element = fe.VectorElement("P", cell, degree[1]) + + temperature_element = fe.FiniteElement("P", cell, degree[2]) - return fe.MixedElement(scalar, vector, scalar) + return fe.MixedElement( + pressure_element, velocity_element, temperature_element) erf, sqrt = fe.erf, fe.sqrt @@ -21,7 +28,7 @@ def liquid_volume_fraction(sim, temperature): T_L = sim.liquidus_temperature - sigma = sim.smoothing + sigma = sim.liquidus_smoothing_factor return 0.5*(1. + erf((T - T_L)/(sigma*sqrt(2.)))) @@ -102,6 +109,17 @@ def strong_residual(sim, solution, buoyancy = linear_boussinesq_buoyancy): return r_p, r_u, r_T +def strong_residual_with_pressure_penalty(sim, solution, buoyancy = linear_boussinesq_buoyancy): + + r_p, r_u, r_T = strong_residual(sim = sim, solution = solution, buoyancy = buoyancy) + + p, _, _ = solution + + gamma = sim.pressure_penalty_factor + + return r_p + gamma*p, r_u, r_T + + def time_discrete_terms(sim): _, u_t, _ = sapphire.simulation.time_discrete_terms( @@ -191,7 +209,7 @@ def energy(sim, solution): + 1./Pr*dot(grad(psi_T), k*grad(T)) -def stabilization(sim, solution): +def pressure_penalty(sim, solution): p, _, _ = fe.split(solution) @@ -206,7 +224,7 @@ def variational_form_residual(sim, solution): return sum( [r(sim = sim, solution = solution) - for r in (mass, momentum, energy, stabilization)])\ + for r in (mass, momentum, energy, pressure_penalty)])\ *fe.dx(degree = sim.quadrature_degree) @@ -216,14 +234,10 @@ def plotvars(sim, solution = None): solution = sim.solution - V = fe.FunctionSpace( - solution.function_space().mesh(), - fe.FiniteElement("P", sim.mesh.ufl_cell(), 1)) - p, u, T = solution.split() phil = fe.interpolate(liquid_volume_fraction( - sim = sim, temperature = T), V) + sim = sim, temperature = T), T.function_space()) return (p, u, T, phil), \ ("p", "\\mathbf{u}", "T", "\\phi_l"), \ @@ -232,27 +246,56 @@ def plotvars(sim, solution = None): class Simulation(sapphire.simulation.Simulation): - def __init__(self, *args, mesh, element_degree = 1, **kwargs): + def __init__( + self, + *args, + mesh, + element_degree = (1, 2, 2), + grashof_number = 1., + prandtl_number = 1., + stefan_number = 1., + pressure_penalty_factor = 1.e-7, + liquidus_temperature = 0., + density_solid_to_liquid_ratio = 1., + heat_capacity_solid_to_liquid_ratio = 1., + thermal_conductivity_solid_to_liquid_ratio = 1., + solid_velocity_relaxation_factor = 1.e-12, + liquidus_smoothing_factor = 0.01, + solver_parameters = { + "snes_type": "newtonls", + "snes_max_it": 24, + "snes_monitor": None, + "snes_rtol": 0., + "ksp_type": "preonly", + "pc_type": "lu", + "mat_type": "aij", + "pc_factor_mat_solver_type": "mumps"}, + **kwargs): - self.grashof_number = fe.Constant(1.) + self.grashof_number = fe.Constant(grashof_number) - self.prandtl_number = fe.Constant(1.) + self.prandtl_number = fe.Constant(prandtl_number) - self.stefan_number = fe.Constant(1.) + self.stefan_number = fe.Constant(stefan_number) - self.pressure_penalty_factor = fe.Constant(1.e-7) + self.pressure_penalty_factor = fe.Constant(pressure_penalty_factor) - self.liquidus_temperature = fe.Constant(0.) + self.liquidus_temperature = fe.Constant(liquidus_temperature) - self.density_solid_to_liquid_ratio = fe.Constant(1.) + self.density_solid_to_liquid_ratio = fe.Constant( + density_solid_to_liquid_ratio) - self.heat_capacity_solid_to_liquid_ratio = fe.Constant(1.) + self.heat_capacity_solid_to_liquid_ratio = fe.Constant( + heat_capacity_solid_to_liquid_ratio) - self.thermal_conductivity_solid_to_liquid_ratio = fe.Constant(1.) + self.thermal_conductivity_solid_to_liquid_ratio = fe.Constant( + thermal_conductivity_solid_to_liquid_ratio) - self.solid_velocity_relaxation_factor = fe.Constant(1.e-12) + self.solid_velocity_relaxation_factor = fe.Constant( + solid_velocity_relaxation_factor) - self.smoothing = fe.Constant(1./256.) + self.liquidus_smoothing_factor = fe.Constant( + liquidus_smoothing_factor) self.smoothing_sequence = None @@ -268,32 +311,19 @@ def __init__(self, *args, mesh, element_degree = 1, **kwargs): mesh = mesh, element = element( cell = mesh.ufl_cell(), degree = element_degree), - **kwargs) - - def solve(self, *args, **kwargs): - - return super().solve(*args, - parameters = { - "snes_type": "newtonls", - "snes_max_it": 24, - "snes_monitor": None, - "snes_rtol": 0., - "ksp_type": "preonly", - "pc_type": "lu", - "mat_type": "aij", - "pc_factor_mat_solver_type": "mumps"}, + solver_parameters = solver_parameters, **kwargs) def solve_with_auto_smoothing(self): - s0 = self.smoothing.__float__() + s0 = self.liquidus_smoothing_factor.__float__() def solve_with_over_regularization(self, startval): return sapphire.continuation.solve_with_over_regularization( solve = self.solve, solution = self.solution, - regularization_parameter = self.smoothing, + regularization_parameter = self.liquidus_smoothing_factor, startval = startval) def solve_with_bounded_regularization_sequence(self): @@ -303,7 +333,7 @@ def solve_with_bounded_regularization_sequence(self): solve = self.solve, solution = self.solution, backup_solution = self.backup_solution, - regularization_parameter = self.smoothing, + regularization_parameter = self.liquidus_smoothing_factor, initial_regularization_sequence = self.smoothing_sequence) if self.smoothing_sequence is None: @@ -311,7 +341,7 @@ def solve_with_bounded_regularization_sequence(self): self.solution, smax = solve_with_over_regularization( self, startval = None) - s = self.smoothing.__float__() + s = self.liquidus_smoothing_factor.__float__() if s == smax: @@ -331,12 +361,13 @@ def solve_with_bounded_regularization_sequence(self): self.solution, smax = solve_with_over_regularization( self, startval = self.smoothing_sequence[-1]) - self.smoothing_sequence = (smax, self.smoothing.__float__()) + self.smoothing_sequence = ( + smax, self.liquidus_smoothing_factor.__float__()) self.solution, self.smoothing_sequence = \ solve_with_bounded_regularization_sequence(self) - assert(self.smoothing.__float__() == s0) + assert(self.liquidus_smoothing_factor.__float__() == s0) return self.solution @@ -348,11 +379,19 @@ def run(self, *args, **kwargs): def postprocess(self): - _, _, T = self.solution.split() + p, u, T = self.solution.split() + + dx = fe.dx(degree = self.quadrature_degree) + + div = fe.div + + self.mean_pressure = fe.assemble(p*dx) + + self.velocity_divergence = fe.assemble(div(u)*dx) phil = liquid_volume_fraction(sim = self, temperature = T) - self.liquid_area = fe.assemble(phil*fe.dx) + self.liquid_area = fe.assemble(phil*dx) return self diff --git a/sapphire/simulations/convection_diffusion.py b/sapphire/simulations/convection_diffusion.py index 5c71811..cdc508d 100644 --- a/sapphire/simulations/convection_diffusion.py +++ b/sapphire/simulations/convection_diffusion.py @@ -41,10 +41,13 @@ def element(cell, degree): class Simulation(sapphire.simulation.Simulation): def __init__(self, *args, - mesh, advection_velocity, element_degree = 1, + mesh, + advection_velocity, + diffusion_coefficient = 1., + element_degree = 1, **kwargs): - self.diffusion_coefficient = fe.Constant(1.) + self.diffusion_coefficient = fe.Constant(diffusion_coefficient) self.advection_velocity = advection_velocity(mesh) diff --git a/sapphire/simulations/heat.py b/sapphire/simulations/heat.py index 8fb4299..36d77ac 100644 --- a/sapphire/simulations/heat.py +++ b/sapphire/simulations/heat.py @@ -37,16 +37,17 @@ def strong_residual(sim, solution): class Simulation(sapphire.simulation.Simulation): - def __init__(self, *args, mesh, element_degree = 1, **kwargs): + def __init__(self, *args, + mesh, + element_degree = 1, + solver_parameters = {"ksp_type": "cg"}, + **kwargs): super().__init__(*args, mesh = mesh, element = element( cell = mesh.ufl_cell(), degree = element_degree), variational_form_residual = variational_form_residual, + solver_parameters = solver_parameters, **kwargs) - - def solve(self, *args, **kwargs): - - return super().solve(*args, parameters = {"ksp_type": "cg"}, **kwargs) - \ No newline at end of file + \ No newline at end of file diff --git a/sapphire/simulations/laplace.py b/sapphire/simulations/laplace.py index fe06569..c294b54 100644 --- a/sapphire/simulations/laplace.py +++ b/sapphire/simulations/laplace.py @@ -32,7 +32,11 @@ def strong_residual(sim, solution): class Simulation(sapphire.simulation.Simulation): - def __init__(self, *args, mesh, element_degree = 1, **kwargs): + def __init__(self, *args, + mesh, + element_degree = 1, + solver_parameters = {"ksp_type": "cg"}, + **kwargs): super().__init__(*args, mesh = mesh, @@ -40,9 +44,6 @@ def __init__(self, *args, mesh, element_degree = 1, **kwargs): cell = mesh.ufl_cell(), degree = element_degree), variational_form_residual = variational_form_residual, time_dependent = False, + solver_parameters = solver_parameters, **kwargs) - - def solve(self, *args, **kwargs): - - return super().solve(*args, parameters = {"ksp_type": "cg"}, **kwargs) \ No newline at end of file diff --git a/sapphire/simulations/navier_stokes_boussinesq.py b/sapphire/simulations/navier_stokes_boussinesq.py index 7c23495..cb82c8e 100644 --- a/sapphire/simulations/navier_stokes_boussinesq.py +++ b/sapphire/simulations/navier_stokes_boussinesq.py @@ -35,9 +35,13 @@ def variational_form_residual( energy = psi_T*dot(u, grad(T)) + dot(grad(psi_T), 1./Pr*grad(T)) + gamma = sim.pressure_penalty_constant + + pressure_penalty = gamma*psi_p*p + dx = fe.dx(degree = sim.quadrature_degree) - return (mass + momentum + energy)*dx + return (mass + momentum + energy + pressure_penalty)*dx def strong_residual(sim, solution, buoyancy = linear_boussinesq_buoyancy): @@ -59,20 +63,35 @@ def strong_residual(sim, solution, buoyancy = linear_boussinesq_buoyancy): def element(cell, degree): - scalar = fe.FiniteElement("P", cell, degree) + if type(degree) is type(1): + + degree = (degree,)*3 + + pressure_element = fe.FiniteElement("P", cell, degree[0]) + + velocity_element = fe.VectorElement("P", cell, degree[1]) - vector = fe.VectorElement("P", cell, degree + 1) + temperature_element = fe.FiniteElement("P", cell, degree[2]) - return fe.MixedElement(scalar, vector, scalar) + return fe.MixedElement( + pressure_element, velocity_element, temperature_element) class Simulation(sapphire.simulation.Simulation): - def __init__(self, *args, mesh, element_degree = 1, **kwargs): + def __init__(self, *args, + mesh, + element_degree = (1, 2, 2), + grashof_number = 1., + prandtl_number = 1., + pressure_penalty_constant = 0., + **kwargs): + + self.grashof_number = fe.Constant(grashof_number) - self.grashof_number = fe.Constant(1.) + self.prandtl_number = fe.Constant(prandtl_number) - self.prandtl_number = fe.Constant(1.) + self.pressure_penalty_constant = fe.Constant(pressure_penalty_constant) super().__init__(*args, mesh = mesh, diff --git a/sapphire/simulations/phasechange.py b/sapphire/simulations/phasechange.py index e15188f..ee79ac2 100644 --- a/sapphire/simulations/phasechange.py +++ b/sapphire/simulations/phasechange.py @@ -12,7 +12,7 @@ def liquid_volume_fraction(sim, temperature): T_m = sim.liquidus_temperature - sigma = sim.smoothing + sigma = sim.liquidus_smoothing_factor return 0.5*(1. + erf((T - T_m)/(sigma*sqrt(2)))) @@ -65,22 +65,27 @@ def element(cell, degree): class Simulation(sapphire.simulation.Simulation): - def __init__(self, *args, mesh, element_degree = 1, **kwargs): + def __init__(self, *args, + mesh, + element_degree = 1, + stefan_number = 1., + liquidus_temperature = 0., + liquidus_smoothing_factor = 0.01, + solver_parameters = {"ksp_type": "cg"}, + **kwargs): - self.stefan_number = fe.Constant(1.) + self.stefan_number = fe.Constant(stefan_number) - self.liquidus_temperature = fe.Constant(0.) + self.liquidus_temperature = fe.Constant(liquidus_temperature) - self.smoothing = fe.Constant(1./32.) + self.liquidus_smoothing_factor = fe.Constant( + liquidus_smoothing_factor) super().__init__(*args, mesh = mesh, element = element( cell = mesh.ufl_cell(), degree = element_degree), variational_form_residual = variational_form_residual, + solver_parameters = solver_parameters, **kwargs) - - def solve(self, *args, **kwargs): - - return super().solve(*args, parameters = {"ksp_type": "cg"}, **kwargs) \ No newline at end of file diff --git a/setup.py b/setup.py index 9925704..cbce0f9 100644 --- a/setup.py +++ b/setup.py @@ -3,6 +3,6 @@ setuptools.setup( name = 'sapphire', - version = '0.5.0a0', + version = '0.5.1a0', packages = setuptools.find_packages(), ) diff --git a/tests/model_validation/test__convection_coupled_phasechange.py b/tests/model_validation/test__convection_coupled_phasechange.py index 88fe3d6..65147c0 100644 --- a/tests/model_validation/test__convection_coupled_phasechange.py +++ b/tests/model_validation/test__convection_coupled_phasechange.py @@ -5,115 +5,44 @@ import sapphire.test -datadir = sapphire.test.datadir +tempdir = sapphire.test.datadir -def test__validate__melt_octadecane__regression(): - - tf = 80. - - th = 40. - - h = -0.02 - - - s = 1./200. - - tau = 1.e-12 - - - rx = 2 - - nx = 24 - - - rt = 2 - - nt = 4 - - - q = 4 - - - expected_liquid_area = 0.43 - - tolerance = 0.01 - +def test__validate__melt_octadecane__regression(tempdir): sim = sapphire.benchmarks.melt_octadecane_in_cavity.Simulation( - quadrature_degree = q, - element_degree = rx - 1, - time_stencil_size = rt + 1, - meshsize = nx, - output_directory_path = "output/melt_octadecane/" - + "th{0}_h{1}/".format(th, h) - + "s{0}_tau{1}/".format(s, tau) - + "rx{0}_nx{1}_rt{2}_nt{3}/".format(rx, nx, rt, nt) - + "q{0}/".format(q)) - + element_degree = (1, 1, 1), + meshsize = 24, + quadrature_degree = 4, + time_stencil_size = 3, + timestep_size = 20., + liquidus_smoothing_factor = 1./200., + solid_velocity_relaxation_factor = 1.e-12, + output_directory_path = tempdir) - sim.timestep_size.assign(tf/float(nt)) + sim.run(endtime = 80.) - sim.smoothing.assign(s) + print("Liquid area = {}".format(sim.liquid_area)) - sim.solid_velocity_relaxation_factor.assign(tau) - - sim.topwall_heatflux_switchtime = th - - sim.topwall_heatflux_postswitch = h - - - sim.solutions, _, = sim.run( - endtime = tf, - topwall_heatflux_poststart = h, - topwall_heatflux_starttime = th) - - - print("Liquid area = {0}".format(sim.liquid_area)) - - assert(abs(sim.liquid_area - expected_liquid_area) < tolerance) + assert(abs(sim.liquid_area - 0.41) < 0.01) -def freeze_water(endtime, s, tau, rx, nx, rt, nt, q, outdir = ""): - - sim = sapphire.benchmarks.freeze_water_in_cavity.Simulation( - quadrature_degree = q, - element_degree = rx - 1, - time_stencil_size = rt + 1, - meshsize = nx, - output_directory_path = str(outdir.join( - "freeze_water/" - + "s{0}_tau{1}/".format(s, tau) - + "rx{0}_nx{1}_rt{2}_nt{3}/".format(rx, nx, rt, nt) - + "q{0}/".format(q)))) - - sim.timestep_size = sim.timestep_size.assign(endtime/float(nt)) +def test__validate__freeze_water__regression(tempdir): - sim.solid_velocity_relaxation_factor = \ - sim.solid_velocity_relaxation_factor.assign(tau) + endtime = 1.44 - sim.smoothing = sim.smoothing.assign(s) - - - sim.solutions, _, = sim.run(endtime = endtime) - - - print("Liquid area = {0}".format(sim.liquid_area)) - - return sim - - -def test__validate__freeze_water__regression(datadir): - - sim = freeze_water( - outdir = datadir, - endtime = 1.44, - s = 1./200., - tau = 1.e-12, - rx = 2, - nx = 24, - rt = 2, - nt = 4, - q = 4) + sim = sapphire.benchmarks.freeze_water_in_cavity.Simulation( + element_degree = (1, 2, 2), + meshsize = 24, + quadrature_degree = 4, + time_stencil_size = 3, + timestep_size = endtime/4., + solid_velocity_relaxation_factor = 1.e-12, + liquidus_smoothing_factor = 1./200., + output_directory_path = tempdir) + + sim.run(endtime = endtime) + + print("Liquid area = {}".format(sim.liquid_area)) assert(abs(sim.liquid_area - 0.69) < 0.01) \ No newline at end of file diff --git a/tests/model_verification/test__convection_coupled_phasechange.py b/tests/model_verification/test__convection_coupled_phasechange.py index 90586d5..3613bb0 100644 --- a/tests/model_verification/test__convection_coupled_phasechange.py +++ b/tests/model_verification/test__convection_coupled_phasechange.py @@ -8,64 +8,80 @@ tempdir = sapphire.test.datadir + pi, sin, exp, dot = fe.pi, fe.sin, fe.exp, fe.dot -def manufactured_solution(mesh, time): - x = fe.SpatialCoordinate(mesh) +def space_verification_solution(sim): - ihat, jhat = sapphire.simulation.unit_vectors(mesh) + x, y = fe.SpatialCoordinate(sim.mesh) - t = time + ihat, jhat = sapphire.simulation.unit_vectors(sim.mesh) - u = exp(t/2)*sin(2.*pi*x[0])*sin(pi*x[1])*ihat + \ - exp(t/2)*sin(pi*x[0])*sin(2.*pi*x[1])*jhat + u = exp(0.5)*sin(2.*pi*x)*sin(pi*y)*ihat + \ + exp(0.5)*sin(pi*x)*sin(2.*pi*y)*jhat - p = -sin(pi*x[0])*sin(2.*pi*x[1]) + # The pressure derivative is zero at x = 0, 1 and y = 0, 1 + # so that Dirichlet BC's do not have to be applied + # and non-homogeneous Neumann BC's don't need to be handled. + p = -sin(pi*x - pi/2.)*sin(2.*pi*y - pi/2.) - T = 0.5*sin(2.*pi*x[0])*sin(pi*x[1])*(1. - exp(-0.5*t**2)) + T = 0.5*sin(2.*pi*x)*sin(pi*y)*(1. - exp(-0.5)) return p, u, T - - -def space_verification_solution(sim): - - return manufactured_solution(mesh = sim.mesh, time = 1.) - + def time_verification_solution(sim): - - return manufactured_solution(mesh = sim.mesh, time = sim.time) - -Gr = 3.6e5 - -Pr = 7.0 - -Ste = 0.13 - -rhos_over_rhol = 0.92 + t = sim.time + + x, y = fe.SpatialCoordinate(sim.mesh) + + ihat, jhat = sapphire.simulation.unit_vectors(sim.mesh) + + u = exp(t/2)*sin(2.*pi*x)*sin(pi*y)*ihat + \ + exp(t/2)*sin(pi*x)*sin(2.*pi*y)*jhat + + # The pressure derivative is zero at x = 0, 1 and y = 0, 1 + # so that Dirichlet BC's do not have to be applied + # and non-homogeneous Neumann BC's don't need to be handled. + p = -sin(pi*x - pi/2.)*sin(2.*pi*y - pi/2.) + + T = 0.5*sin(2.*pi*x)*sin(pi*y)*(1. - exp(-0.5*t**2)) + + return p, u, T + -cs_over_cl = 0.50 +parameters = { + "grashof_number": 3.6e5, + "prandtl_number": 7.0, + "stefan_number": 0.13, + "density_solid_to_liquid_ratio": 0.92, + "heat_capacity_solid_to_liquid_ratio": 0.50, + "thermal_conductivity_solid_to_liquid_ratio": 3.8, + "liquidus_smoothing_factor": 0.1, + "solid_velocity_relaxation_factor": 1.e-12, + "pressure_penalty_factor": 1.e-4, + "quadrature_degree": None, + "element_degree": None, + } + -kappas_over_kappal = 3.8 -sigma = 0.1 endtime = 1. -def test__verify__second_order_spatial_convergence__via_mms( - tempdir, - parameters = { - "grashof_number": Gr, - "prandtl_number": Pr, - "stefan_number": Ste, - "density_solid_to_liquid_ratio": rhos_over_rhol, - "heat_capacity_solid_to_liquid_ratio": cs_over_cl, - "thermal_conductivity_solid_to_liquid_ratio": kappas_over_kappal, - "smoothing": sigma}, - mesh_sizes = (4, 8, 16, 32, 64), - tolerance = 0.1): +def test__verify__taylor_hood_second_order_spatial_convergence__via_mms( + tempdir): + """ + Demonstrate second order accuracy for velocity and temperature fields. + Pressure error shows superconvergence until nx = 64 + where the magnitude hits a floor around 0.5. + To keep the test cheap, only up to nx = 32 is shown here. + """ + parameters["element_degree"] = (1, 2, 2) + + parameters["timestep_size"] = endtime testdir = "{}/{}/".format( __name__.replace(".", "/"), sys._getframe().f_code.co_name) @@ -78,29 +94,81 @@ def test__verify__second_order_spatial_convergence__via_mms( sapphire.mms.verify_spatial_order_of_accuracy( sim_module = sim_module, + sim_parameters = parameters, manufactured_solution = space_verification_solution, - meshes = [fe.UnitSquareMesh(n, n) for n in mesh_sizes], - parameters = parameters, - expected_order = 2, - tolerance = tolerance, - timestep_size = endtime, + #strong_residual = sim_module.strong_residual_with_pressure_penalty, #Adding the penalty term to the strong residual removes the floor and maintains superconvergence. + meshes = [fe.UnitSquareMesh(n, n) for n in (4, 8, 16, 32)], + norms = ("L2", "H1", "H1"), + expected_orders = (None, 2, 2), + tolerance = 0.1, endtime = endtime, outfile = outfile) + + +def test__verify__taylor_hood_third_order_spatial_convergence__via_mms( + tempdir): + + parameters["element_degree"] = (2, 3, 3) + + parameters["timestep_size"] = endtime + + testdir = "{}/{}/".format( + __name__.replace(".", "/"), sys._getframe().f_code.co_name) + + outdir_path = pathlib.Path(tempdir) / testdir + + outdir_path.mkdir(parents = True, exist_ok = True) + + with open(outdir_path / "convergence.csv", "w") as outfile: + sapphire.mms.verify_spatial_order_of_accuracy( + sim_module = sim_module, + sim_parameters = parameters, + manufactured_solution = space_verification_solution, + meshes = [fe.UnitSquareMesh(n, n) for n in (4, 8, 16, 32)], + norms = ("L2", "H1", "H1"), + expected_orders = (None, 3, 3), + tolerance = 0.1, + endtime = endtime, + outfile = outfile) + + +def test__verify__equal_order_first_order_spatial_convergence__via_mms( + tempdir): + + parameters["element_degree"] = (1, 1, 1) + + parameters["timestep_size"] = endtime + + testdir = "{}/{}/".format( + __name__.replace(".", "/"), sys._getframe().f_code.co_name) + + outdir_path = pathlib.Path(tempdir) / testdir + + outdir_path.mkdir(parents = True, exist_ok = True) + + with open(outdir_path / "convergence.csv", "w") as outfile: + sapphire.mms.verify_spatial_order_of_accuracy( + sim_module = sim_module, + sim_parameters = parameters, + manufactured_solution = space_verification_solution, + meshes = [fe.UnitSquareMesh(n, n) for n in (4, 8, 16, 32)], + norms = ("L2", "H1", "H1"), + expected_orders = (None, 1, 1), + tolerance = 0.1, + endtime = endtime, + outfile = outfile) + + def test__verify__second_order_temporal_convergence__via_mms( - tempdir, - parameters = { - "grashof_number": Gr, - "prandtl_number": Pr, - "stefan_number": Ste, - "density_solid_to_liquid_ratio": rhos_over_rhol, - "heat_capacity_solid_to_liquid_ratio": cs_over_cl, - "thermal_conductivity_solid_to_liquid_ratio": kappas_over_kappal, - "smoothing": sigma}, - meshsize = 48, - timestep_sizes = (1/2, 1/4, 1/8, 1/16), - tolerance = 0.2): + tempdir): + + parameters["element_degree"] = (1, 2, 2) + + meshsize = 32 + + parameters["mesh"] = fe.UnitSquareMesh(meshsize, meshsize) testdir = "{}/{}/".format( __name__.replace(".", "/"), sys._getframe().f_code.co_name) @@ -113,11 +181,11 @@ def test__verify__second_order_temporal_convergence__via_mms( sapphire.mms.verify_temporal_order_of_accuracy( sim_module = sim_module, + sim_parameters = parameters, manufactured_solution = time_verification_solution, - mesh = fe.UnitSquareMesh(meshsize, meshsize), - parameters = parameters, - expected_order = 2, - tolerance = tolerance, - timestep_sizes = timestep_sizes, + norms = ("L2", "L2", "L2"), + expected_orders = (None, 2, 2), + tolerance = 0.2, + timestep_sizes = (1/2, 1/4, 1/8, 1/16), endtime = endtime, outfile = outfile) diff --git a/tests/model_verification/test__convection_diffusion.py b/tests/model_verification/test__convection_diffusion.py index e8e2d8e..e595361 100644 --- a/tests/model_verification/test__convection_diffusion.py +++ b/tests/model_verification/test__convection_diffusion.py @@ -22,14 +22,16 @@ def advection_velocity(mesh): + sin(pi*x[0])*sin(2.*pi*x[1])*jhat -def test__verify_convergence_order_via_mms( - mesh_sizes = (8, 16, 32), tolerance = 0.1): +def test__verify_convergence_order_via_mms(): sapphire.mms.verify_spatial_order_of_accuracy( sim_module = sim_module, + sim_parameters = { + "diffusion_coefficient": 0.1, + "advection_velocity": advection_velocity, + }, manufactured_solution = manufactured_solution, - meshes = [fe.UnitSquareMesh(n, n) for n in mesh_sizes], - sim_constructor_kwargs = {"advection_velocity": advection_velocity}, - parameters = {"diffusion_coefficient": 0.1}, - expected_order = 2, - tolerance = tolerance) + meshes = [fe.UnitSquareMesh(n, n) for n in (8, 16, 32)], + norms = ("H1",), + expected_orders = (1,), + tolerance = 0.1) diff --git a/tests/model_verification/test__heat.py b/tests/model_verification/test__heat.py index 125b407..bf6337d 100644 --- a/tests/model_verification/test__heat.py +++ b/tests/model_verification/test__heat.py @@ -14,68 +14,62 @@ def manufactured_solution(sim): return sin(2.*pi*x)*exp(-pow(t, 2)) -def test__verify_spatial_convergence__second_order__via_mms( - mesh_sizes = (4, 8, 16, 32), - timestep_size = 1./64., - tolerance = 0.1): +def test__verify_spatial_convergence__first_order__via_mms(): sapphire.mms.verify_spatial_order_of_accuracy( sim_module = sim_module, + sim_parameters = {"timestep_size": 1./64.}, manufactured_solution = manufactured_solution, - meshes = [fe.UnitIntervalMesh(n) for n in mesh_sizes], - expected_order = 2, - tolerance = tolerance, - timestep_size = timestep_size, + meshes = [fe.UnitIntervalMesh(n) for n in (4, 8, 16, 32)], + norms = ("H1",), + expected_orders = (1,), + tolerance = 0.1, endtime = 1.) -def test__verify_temporal_convergence__first_order__via_mms( - meshsize = 256, - timestep_sizes = (1./4., 1./8., 1./16., 1./32.), - tolerance = 0.1): +def test__verify_temporal_convergence__first_order__via_mms(): sapphire.mms.verify_temporal_order_of_accuracy( sim_module = sim_module, + sim_parameters = {"mesh": fe.UnitIntervalMesh(256)}, manufactured_solution = manufactured_solution, - mesh = fe.UnitIntervalMesh(meshsize), - expected_order = 1, + norms = ("L2",), + expected_orders = (1,), endtime = 1., - timestep_sizes = timestep_sizes, - tolerance = tolerance) + timestep_sizes = (1./4., 1./8., 1./16., 1./32.), + tolerance = 0.1) -def test__verify_temporal_convergence__second_order__via_mms( - meshsize = 128, - timestep_sizes = (1./16., 1./32., 1./64., 1./128.), - tolerance = 0.1): +def test__verify_temporal_convergence__second_order__via_mms(): sapphire.mms.verify_temporal_order_of_accuracy( sim_module = sim_module, - manufactured_solution = manufactured_solution, - mesh = fe.UnitIntervalMesh(meshsize), - sim_constructor_kwargs = { + sim_parameters = { + "mesh": fe.UnitIntervalMesh(128), "element_degree": 2, - "time_stencil_size": 3}, - expected_order = 2, + "time_stencil_size": 3, + }, + manufactured_solution = manufactured_solution, + norms = ("L2",), + expected_orders = (2,), endtime = 1., - timestep_sizes = timestep_sizes, - tolerance = tolerance) + timestep_sizes = (1./16., 1./32., 1./64., 1./128.), + tolerance = 0.1) -def test__verify_temporal_convergence__third_order__via_mms( - meshsize = 128, - timestep_sizes = (1./4., 1./8., 1./16., 1./32.), - tolerance = 0.1): +def test__verify_temporal_convergence__third_order__via_mms(): sapphire.mms.verify_temporal_order_of_accuracy( sim_module = sim_module, - manufactured_solution = manufactured_solution, - mesh = fe.UnitIntervalMesh(meshsize), - sim_constructor_kwargs = { + sim_parameters = { + "mesh": fe.UnitIntervalMesh(128), "element_degree": 2, - "time_stencil_size": 4}, - expected_order = 3, + "time_stencil_size": 4, + }, + manufactured_solution = manufactured_solution, + norms = ("L2",), + expected_orders = (3,), endtime = 1., - timestep_sizes = timestep_sizes, - tolerance = tolerance) + timestep_sizes = (1./4., 1./8., 1./16., 1./32.), + tolerance = 0.1) \ No newline at end of file diff --git a/tests/model_verification/test__laplace.py b/tests/model_verification/test__laplace.py index 14fc343..9c10bc8 100644 --- a/tests/model_verification/test__laplace.py +++ b/tests/model_verification/test__laplace.py @@ -19,6 +19,7 @@ def test__verify_convergence_order_via_mms( sim_module = sim_module, manufactured_solution = manufactured_solution, meshes = [fe.UnitIntervalMesh(n) for n in mesh_sizes], - expected_order = 2, + norms = ("L2",), + expected_orders = (2,), tolerance = tolerance) \ No newline at end of file diff --git a/tests/model_verification/test__navier_stokes.py b/tests/model_verification/test__navier_stokes.py index 9feffd4..aad059d 100644 --- a/tests/model_verification/test__navier_stokes.py +++ b/tests/model_verification/test__navier_stokes.py @@ -26,6 +26,7 @@ def test__verify_convergence_order_via_mms( sim_module = sim_module, manufactured_solution = manufactured_solution, meshes = [fe.UnitSquareMesh(n, n) for n in mesh_sizes], - expected_order = 2, + norms = ("H1", "L2"), + expected_orders = (2, 2), tolerance = tolerance) \ No newline at end of file diff --git a/tests/model_verification/test__navier_stokes_boussinesq.py b/tests/model_verification/test__navier_stokes_boussinesq.py index 1ac605f..92e5da7 100644 --- a/tests/model_verification/test__navier_stokes_boussinesq.py +++ b/tests/model_verification/test__navier_stokes_boussinesq.py @@ -25,7 +25,7 @@ def manufactured_solution(sim): return p, u, T -def test__verify_convergence_order_via_mms( +def test__verify_taylor_hood_accuracy_via_mms( mesh_sizes = (4, 8, 16, 32), tolerance = 0.1): @@ -35,10 +35,65 @@ def test__verify_convergence_order_via_mms( sapphire.mms.verify_spatial_order_of_accuracy( sim_module = sim_module, + sim_parameters = { + "grashof_number": Ra/Pr, + "prandtl_number": Pr, + "element_degree": (1, 2, 2), + }, manufactured_solution = manufactured_solution, meshes = [fe.UnitSquareMesh(n, n) for n in mesh_sizes], - parameters = {"grashof_number": Ra/Pr, "prandtl_number": Pr}, - expected_order = 2, + norms = ("L2", "H1", "H1"), + expected_orders = (2, 2, 2), + tolerance = tolerance) + + +def test__show_equal_order_pressure_penalty_inaccuracy_via_mms( + mesh_sizes = (4, 8, 16, 32, 64), + tolerance = 0.1): + + Ra = 10. + + Pr = 0.7 + + gamma = 1.e-7 + + sapphire.mms.verify_spatial_order_of_accuracy( + sim_module = sim_module, + sim_parameters = { + "grashof_number": Ra/Pr, + "prandtl_number": Pr, + "pressure_penalty_constant": gamma, + "element_degree": (1, 1, 1), + }, + manufactured_solution = manufactured_solution, + meshes = [fe.UnitSquareMesh(n, n) for n in mesh_sizes], + norms = ("L2", "H1", "H1"), + expected_orders = (None, 1, 1), + tolerance = tolerance) + + +def test__verify_taylor_hood_with_pressure_penalty_accuracy_via_mms( + mesh_sizes = (4, 8, 16, 32), + tolerance = 0.1): + + Ra = 10. + + Pr = 0.7 + + gamma = 1.e-7 + + sapphire.mms.verify_spatial_order_of_accuracy( + sim_module = sim_module, + sim_parameters = { + "grashof_number": Ra/Pr, + "prandtl_number": Pr, + "pressure_penalty_constant": gamma, + "element_degree": (1, 2, 2), + }, + manufactured_solution = manufactured_solution, + meshes = [fe.UnitSquareMesh(n, n) for n in mesh_sizes], + norms = ("L2", "H1", "H1"), + expected_orders = (2, 2, 2), tolerance = tolerance) @@ -107,7 +162,7 @@ def verify_scalar_solution_component( def test__verify_against_heat_driven_cavity_benchmark(): sim = sapphire.benchmarks.heat_driven_cavity.Simulation( - element_degree = 1, meshsize = 40) + element_degree = (1, 2, 2), meshsize = 40) sim.solution = sim.solve() diff --git a/tests/model_verification/test__phasechange.py b/tests/model_verification/test__phasechange.py index 6b4b7ad..5189beb 100644 --- a/tests/model_verification/test__phasechange.py +++ b/tests/model_verification/test__phasechange.py @@ -14,53 +14,73 @@ def manufactured_solution(sim): return 0.5*sin(2.*pi*x)*(1. - 2*exp(-3.*t**2)) -def test__verify_spatial_convergence__second_order__via_mms( - mesh_sizes = (4, 8, 16, 32), - timestep_size = 1./256., - tolerance = 0.1): +def test__verify_spatial_convergence__first_order__via_mms(): sapphire.mms.verify_spatial_order_of_accuracy( sim_module = sim_module, + sim_parameters = { + "stefan_number": 0.1, + "liquidus_smoothing_factor": 1./32., + "timestep_size": 1./256., + }, manufactured_solution = manufactured_solution, - meshes = [fe.UnitIntervalMesh(size) for size in mesh_sizes], - parameters = {"stefan_number": 0.1, "smoothing": 1./32.}, - expected_order = 2, - timestep_size = timestep_size, + meshes = [fe.UnitIntervalMesh(size) for size in (4, 8, 16, 32)], + norms = ("H1",), + expected_orders = (1,), endtime = 1., - tolerance = tolerance) + tolerance = 0.1) + +def test__verify_spatial_convergence__second_order__via_mms(): + + sapphire.mms.verify_spatial_order_of_accuracy( + sim_module = sim_module, + sim_parameters = { + "stefan_number": 0.1, + "liquidus_smoothing_factor": 1./32., + "element_degree": 2, + "timestep_size": 1./256., + }, + manufactured_solution = manufactured_solution, + meshes = [fe.UnitIntervalMesh(size) for size in (4, 8, 16, 32)], + norms = ("H1",), + expected_orders = (2,), + endtime = 1., + tolerance = 0.1) -def test__verify_temporal_convergence__first_order__via_mms( - meshsize = 256, - timestep_sizes = (1./16., 1./32., 1./64., 1./128.), - tolerance = 0.1): + +def test__verify_temporal_convergence__first_order__via_mms(): sapphire.mms.verify_temporal_order_of_accuracy( sim_module = sim_module, + sim_parameters = { + "stefan_number": 0.1, + "liquidus_smoothing_factor": 1./32., + "mesh": fe.UnitIntervalMesh(256), + }, manufactured_solution = manufactured_solution, - mesh = fe.UnitIntervalMesh(meshsize), - parameters = {"stefan_number": 0.1, "smoothing": 1./32.}, - expected_order = 1, + norms = ("L2",), + expected_orders = (1,), endtime = 1., - timestep_sizes = timestep_sizes, - tolerance = tolerance) + 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.1): +def test__verify_temporal_convergence__second_order__via_mms(): sapphire.mms.verify_temporal_order_of_accuracy( sim_module = sim_module, - manufactured_solution = manufactured_solution, - mesh = fe.UnitIntervalMesh(meshsize), - sim_constructor_kwargs = { + sim_parameters = { + "stefan_number": 0.1, + "liquidus_smoothing_factor": 1./32., "element_degree": 2, - "time_stencil_size": 3}, - parameters = {"stefan_number": 0.1, "smoothing": 1./32.}, - expected_order = 2, + "time_stencil_size": 3, + "mesh": fe.UnitIntervalMesh(128), + }, + manufactured_solution = manufactured_solution, + norms = ("L2",), + expected_orders = (2,), endtime = 1., - timestep_sizes = timestep_sizes, - tolerance = tolerance) + timestep_sizes = (1./64., 1./128., 1./256.), + tolerance = 0.1) \ No newline at end of file diff --git a/tests/model_verification/test__unsteady_navier_stokes.py b/tests/model_verification/test__unsteady_navier_stokes.py index 977eda5..4b05237 100644 --- a/tests/model_verification/test__unsteady_navier_stokes.py +++ b/tests/model_verification/test__unsteady_navier_stokes.py @@ -22,40 +22,39 @@ def manufactured_solution(sim): return u, p -def test__verify_spatial_convergence__second_order__via_mms( - mesh_sizes = (3, 6, 12, 24), - timestep_size = 1./32., - tolerance = 0.3): +parameters = { + "quadrature_degree": 4, + "element_degree": 1, + "time_stencil_size": 2} + +def test__verify_spatial_convergence__second_order__via_mms(): + + parameters["timestep_size"] = 1./32. sapphire.mms.verify_spatial_order_of_accuracy( sim_module = sim_module, + sim_parameters = parameters, manufactured_solution = manufactured_solution, - meshes = [fe.UnitSquareMesh(n, n) for n in mesh_sizes], - sim_constructor_kwargs = { - "quadrature_degree": 4, - "element_degree": 1, - "time_stencil_size": 2}, - expected_order = 2, - tolerance = tolerance, - timestep_size = timestep_size, + meshes = [fe.UnitSquareMesh(n, n) for n in (3, 6, 12, 24)], + norms = ("H1", "L2"), + expected_orders = (2, 2), + tolerance = 0.3, endtime = 1.) -def test__verify_temporal_convergence__first_order__via_mms( - meshsize = 32, - timestep_sizes = (1./2., 1./4., 1./8.), - tolerance = 0.1): +def test__verify_temporal_convergence__first_order__via_mms(): + + meshsize = 32 + + parameters["mesh"] = fe.UnitSquareMesh(meshsize, meshsize) sapphire.mms.verify_temporal_order_of_accuracy( sim_module = sim_module, + sim_parameters = parameters, manufactured_solution = manufactured_solution, - mesh = fe.UnitSquareMesh(meshsize, meshsize), - sim_constructor_kwargs = { - "quadrature_degree": 4, - "element_degree": 1, - "time_stencil_size": 2}, - expected_order = 1, + norms = ("L2", "L2"), + expected_orders = (None, 1), endtime = 1., - timestep_sizes = timestep_sizes, - tolerance = tolerance) + timestep_sizes = (1./2., 1./4., 1./8.), + tolerance = 0.1) \ No newline at end of file