p-adaptive radial dam break¶
In this example, we want to present a p-adaptive simulation of a radial dam beak on a randomly perturbed uniform mesh [FA2022].
In p-adaptive computations, the local approximation order of each element is dynamically adapted throughout the simulation to combine a good solution quality with an acceptable runtime.
In this case, we want to vary between constant and linear approximations. To decide which approximation order is appropriate for each element, an adaptivity indicator is evaluated in each time step. For demonstration purposes, we use an indicator comparing the jump of the free-surface elevation with a user-defined threshold.
Additionally, in this tutorial, we want to show how to generate a CUDA code that is run on a GPU.
Example setup and parameter choice¶
We choose the domain \(\Omega = [0,5]\times [0,5]\) and \(g=1\). The constant bathymetry is set to \(h_b =0\). And the initial condition is
We use a randomly perturbed uniform mesh with 32 768 triangles. Since all external boundaries use land boundary conditions, the wave is reflected at the walls.
The simulation is run for 15 seconds with a constant time step of 0.0025 seconds for the p-adaptive constant-linear DG discretization (\(p=0-1\)).
Now we look at the input parameters for this setup.
The parameters are set in the file src/config/dambreak.py.
Below we describe the most important parameters to be set for running this example:
Time step size and simulation time:
config.default_dt = 0.0025 # time step size in seconds
config.t_max = 3
Adaptivity and approximation order:
config.adapt_local_order = True
config.min_order = 0
config.max_order = 1
config.adaptivity_nlock = 10 # avoid frequent in- and decreases of local order
config.adaptivity_criterion = f"eta_jumps" # no_adaption: no adaption; eta: eta values being outside of certain range;
# element_grad: element gradients exceeding a threshold (Kubatko et al. 2009);
# neigh_grad: gradients to neighbour exceeding threshold;
# diff_to_lower_order: L^2 norm of the difference between the approximation of the water depth of order o and o - 1 (Esskilsson 2010)
# eta_jumps: sum of jumps in eta (Remacle2003);
# diff_to_lower_order_and_eta_jumps: combination of diff_to_lower_order and eta_jumps
config.adaptivity_threshold_eta_jump = 0.0003 # indicator threshold
config.use_float_for_indicator_evaluation = False # impacts indicator sensitivity and reproducibility of results in some cases
Grid and bathymetry: Randomly perturbed uniform mesh of size 5x5:
config.read_grid_from_file = False
config.read_bath_from_file = False
config.bathymetry = lambda pos: sp.S(0.0)
# size
config.uniform_domain_size_x = sp.S(5)
config.uniform_domain_size_y = sp.S(5)
# number of fragments and blocks
config.grid_gen_frag_scale_x = 1
config.grid_gen_frag_scale_y = 1
config.grid_gen_num_frag_per_block_x = 1
config.grid_gen_num_frag_per_block_y = 1
config.grid_gen_num_block_x = 1
config.grid_gen_num_block_y = 1
config.n_blocks = config.grid_gen_num_block_x * config.grid_gen_num_block_y # derived value
config.n_frags_per_block = config.grid_gen_num_frag_per_block_x * config.grid_gen_num_frag_per_block_y # derived value
config.minLevel = 7
config.maxLevel = config.minLevel # 2*2**maxLevel is the number of triangular elements per dimension
config.assume_uniform_grid = False # randomly perturbed mesh
config.store_normals = not config.assume_uniform_grid # pre-compute and store normals
config.store_edge_lengths = not config.assume_uniform_grid # pre-compute and store edge lenghts
Initial and boundary condition:
def compose_cond(orientation):
[[x0, y0], [x1, y1], [x2, y2]] = [node_position(Neighbors.elem_to_node(orientation, i)) for i in range(3)]
cent_x = (x0 + x1 + x2) / 3
cent_y = (y0 + y1 + y2) / 3
return Compare((cent_x - 2.5) * (cent_x - 2.5) + (cent_y - 2.5) * (cent_y - 2.5), 0.25, "<")
config.exact_eta = {o: [(compose_cond(o), lambda pos, t: 2 + 0.5 * sp.exp(15 * (-(pos[0] - 2.5) * (pos[0] - 2.5) - (pos[1] - 2.5) * (pos[1] - 2.5)))),
(Not(compose_cond(o)), lambda pos, t: sp.S(1))]
for o in TriOrientation}
config.exact_u = lambda pos, t: sp.S(0)
config.exact_v = lambda pos, t: sp.S(0)
# use land boundary condition on all 4 sides
config.available_boundary_conditions = [BoundaryCondition.Land]
config.grid_override_bc = [(0, 0, -100), (1, 0, -100), (2, 0, -100), (3, 0, -100)] # form: [(neigh_idx, old_bc_id, new_bc_id)]
Right-hand side forcings:
# coriolis option parameter 0: spatially constant, 1: spatially variable
config.ncor = 0
# constant coriolis coefficient
config.cori = sp.S(0.00001)
# linearized friction coefficient
config.tau = sp.S(0.0001)
# gravitational acceleration
config.default_g = sp.S(1)
Minimum-depth restriction to prevent elements from getting dry:
config.restrict_min_depth = True
config.min_depth_threshold = sp.Rational(0.01)
config.min_depth_zero_velocities = True # set velocity components to zero in case of min depth restriction
Parallelization and GPU usage:
config.omp_enabled = False # switch off OpenMP
config.use_cuda = True
config.use_cuda_condition = False # generate CUDA-code for everything
config.use_likwid = True # likwid pinning
vtk and recording station output:
config.print_interval = 0.1 # in seconds
config.print_vtk = True
Running the code and visualizing results¶
Setting the correct major and minor CUDA capacity in the .platform file is important. The code is run with the command python3 sweMain.py config.dambreak. It automatically generates layer 4 code and, subsequently the optimized CUDA code.
Here, we are using the Warp By Scalar filter of ParaView to create a 3D plot of the free surface elevation \(\xi\) . In the top figure, the color coding shows \(\xi\) and in the bottom one, the local approximation order.
Faghih-Naini and V. Aizinge. p-adaptive discontinuous Galerkin method for the shallow water equations with a parameter-free error indicator. In: International Journal on Geomathematics (GEM) Vol. 13 (2022) Issue 1 No. 18. doi.org/10.1007/s13137-022-00208-3.