IMEX: Isentropic Vortex on a Geometry-Split Mesh
This notebook demonstrates the geometry-split IMEX scheme on a periodic unit square \([-0.5, 0.5]^2\). The mesh is partitioned into a fine implicit region in the centre (HDG, solved implicitly) and two coarse explicit regions on the sides (DG, advanced explicitly). An isentropic vortex is injected as the initial condition and advected across the domain.
Scheme pair used: ARS-333 — 'sdirk33' (implicit HDG) + 'rk_ars33' (explicit DG), order 3.
[1]:
import numpy as np
import ngsolve as ngs
from dream.compressible_flow import CompressibleFlowSolver, flowfields, Initial, InterfaceBC, FarField
from dream.mesh import get_rectangular_mesh
from dream.time import SynchronizedIMEXTimeRoutine
ngs.SetNumThreads(12)
Mesh
Build a structured rectangular mesh with three named domains:
'explicit'— coarse regions on the left and right'implicit'— fine region in the centre (\(x \in [-0.125, 0.125]\))
Periodic boundary conditions are applied in \(x\); top and bottom boundaries are farfield.
[2]:
Nx_coarse = 10
Nx_fine = 20
Ny = 16
xl = np.linspace(-0.5, -0.125, Nx_coarse // 2 + 1)
xm = np.linspace(-0.125, 0.125, Nx_fine + 1)
xr = np.linspace( 0.125, 0.5, Nx_coarse // 2 + 1)
y = np.linspace(-0.5, 0.5, Ny + 1)
domains = (
('explicit', (xl, y)),
('implicit', (xm, y)),
('explicit', (xr, y)),
)
boundaries = (
('bottom', ((xl.min(), xr.max()), y.min())),
('right', (xr.max(), y)),
('top', ((xl.min(), xr.max()), y.max())),
('left', (xl.min(), y)),
('interface', (xm.min(), y)),
('interface', (xm.max(), y)),
)
mesh = get_rectangular_mesh(domains, boundaries, quads=True, periodic_x=True, periodic_y=False)
implicit_mesh = ngs.Mesh(mesh.ngmesh.GetSubMesh('implicit', 'implicit'))
explicit_mesh = ngs.Mesh(mesh.ngmesh.GetSubMesh('explicit', 'explicit'))
Implicit Solver (HDG, sdirk33)
The implicit sub-solver runs on implicit_mesh using a 3rd-order SDIRK scheme. A Newton solver with frozen Jacobian (reused each time step) is used.
[4]:
cfg_imp = CompressibleFlowSolver(implicit_mesh, **cfg.to_dict())
cfg_imp.fem = 'conservative_hdg'
cfg_imp.fem.scheme = 'sdirk33'
cfg_imp.fem.scheme.compile = True
cfg_imp.fem.order = 3
cfg_imp.fem.bonus_int_order = 6
cfg_imp.fem.viscous_treatment = None
cfg_imp.fem.solver = 'direct'
cfg_imp.fem.solver.method = 'newton'
cfg_imp.fem.solver.method.damping_factor = 1
cfg_imp.fem.solver.method.max_iterations = 20
cfg_imp.fem.solver.method.freeze_jacobian = 'step'
cfg_imp.fem.solver.method.convergence_criterion = 1e-10
Explicit Solver (DG, rk_ars33)
The explicit sub-solver runs on explicit_mesh using the matching ARS-333 ERK scheme.
[5]:
cfg_exp = CompressibleFlowSolver(explicit_mesh, **cfg.to_dict())
cfg_exp.fem = 'conservative_dg'
cfg_exp.fem.scheme = 'rk_ars33'
cfg_exp.fem.scheme.compile = True
cfg_exp.fem.order = 3
cfg_exp.fem.bonus_int_order = 6
cfg_exp.fem.viscous_treatment = None
Initial Condition: Isentropic Vortex
An isentropic vortex centred at the origin is superimposed on the uniform farfield state \(\mathbf{U}_\infty\) at Mach 0.3. The vortex Mach number is \(M_t = 0.01\) and the core radius is \(R = 0.1\).
[6]:
Uinf = cfg_imp.get_farfield_fields((1, 0))
M = cfg_imp.mach_number
gamma = cfg_imp.equation_of_state.heat_capacity_ratio
Mt = 0.01 # vortex Mach number
R = 0.1 # vortex core radius
r = ngs.sqrt(ngs.x**2 + ngs.y**2)
vt = Mt / M * cfg_imp.scaling.velocity
psi = vt * R * ngs.exp((R**2 - r**2) / (2 * R**2))
u_0 = Uinf.u + ngs.CF((psi.Diff(ngs.y), -psi.Diff(ngs.x)))
p_0 = Uinf.p * (1 - (gamma - 1) / 2 * Mt**2 * ngs.exp((R**2 - r**2) / R**2))**(gamma / (gamma - 1))
rho_0 = Uinf.rho * (1 - (gamma - 1) / 2 * Mt**2 * ngs.exp((R**2 - r**2) / R**2))**(1 / (gamma - 1))
initial_fields = flowfields(rho=rho_0, u=u_0, p=p_0)
Boundary and Initial Conditions
Interface boundary conditions link the two sub-solvers across \(\Gamma_i\). The fields attribute is populated after the gridfunctions are created.
[7]:
imp_interface = InterfaceBC(fields=None)
cfg_imp.dcs['implicit'] = Initial(initial_fields)
cfg_imp.bcs['interface'] = imp_interface
cfg_imp.bcs['top|bottom'] = FarField(Uinf)
exp_interface = InterfaceBC(fields=None)
cfg_exp.dcs['explicit'] = Initial(initial_fields)
cfg_exp.bcs['interface'] = exp_interface
cfg_exp.bcs['top|bottom'] = FarField(Uinf)
cfg_exp.bcs['left|right'] = 'periodic'
Initialise, Link, and Solve
SynchronizedIMEXTimeRoutine advances both sub-solvers stage-by-stage, exchanging interface data at each stage so that both sides sit at the same physical time.
[8]:
with ngs.TaskManager():
for cfg_sub in [cfg_imp, cfg_exp]:
cfg_sub.fem.initialize_finite_element_spaces()
cfg_sub.fem.initialize_trial_and_test_functions()
cfg_sub.fem.initialize_gridfunctions()
cfg_sub.fem.initialize_time_scheme_gridfunctions()
# Link interface states after gridfunctions exist
imp_interface.fields = cfg_exp.get_all_solution_fields()
exp_interface.fields = cfg_imp.get_all_solution_fields()
with ngs.TaskManager():
for cfg_sub in [cfg_imp, cfg_exp]:
cfg_sub.fem.set_boundary_conditions()
cfg_sub.fem.set_initial_conditions()
cfg_sub.fem.initialize_symbolic_forms()
[9]:
from ngsolve.webgui import Draw
drawing_imp = cfg_imp.get_solution_fields('p', default_fields=False)
drawing_exp = cfg_exp.get_solution_fields('p', default_fields=False)
p_ref = Uinf.p * (1 - (gamma - 1) / 2 * Mt**2 * ngs.exp(1))**(gamma / (gamma - 1))
drawing_imp['p*'] = (drawing_imp.p - Uinf.p) / (p_ref - Uinf.p)
drawing_exp['p*'] = (drawing_exp.p - Uinf.p) / (p_ref - Uinf.p)
p_star = mesh.MaterialCF({"implicit": drawing_imp['p*'], "explicit": drawing_exp['p*']})
Draw(p_star, mesh)
[9]:
WebGLScene
[10]:
times = [0.0, 0.5, 1.0]
for i in range(len(times) - 1):
cfg_imp.time.timer.interval = (times[i], times[i + 1])
cfg_exp.time.timer.interval = (times[i], times[i + 1])
with ngs.TaskManager():
time_routine = SynchronizedIMEXTimeRoutine(cfg_imp, cfg_exp)
time_routine.solve()
print(f"Completed time interval [{times[i]}, {times[i + 1]}]")
p_star = mesh.MaterialCF({"implicit": drawing_imp['p*'], "explicit": drawing_exp['p*']})
Draw(p_star, mesh)
Completed time interval [0.0, 0.5]
Completed time interval [0.5, 1.0]