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