.. include:: ../../global.rst

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 :math:`\Omega = [0,5]\times [0,5]` and :math:`g=1`. The constant
bathymetry is set to :math:`h_b =0`. And the initial condition is

.. math::

   \begin{aligned}
   &\xi(x, y ,t)=\begin{cases}2+0.5\, e^{-15((x-2.5)^2+(y-2.5)^2)}, &\text{if } (x-2.5)^2+(y-2.5)^2<0.25,\\
   1,& \text{otherwise},
   \end{cases}\\
   &U(x,y,t)= 0, \; V(x,y,t) = 0.
   \end{aligned}

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
(:math:`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:

.. code:: python

    	config.default_dt = 0.0025 # time step size in seconds
        config.t_max = 3
        
Adaptivity and approximation order:

.. code:: python

	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:

.. code:: python
	
	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:

.. code:: python

        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:

.. code:: python

        # 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:

.. code:: python

        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:

.. code:: python

        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:

.. code:: python

        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 :math:`\xi` . In the top figure, the color coding shows  :math:`\xi`  and in the bottom one, the local approximation order.

.. raw:: html

	<video class="hide-in-ie8" width="720" height="394" autoplay muted controls loop>
		<source src="../../_static/_videos/dambreak.ogv" type="video/ogg" /> 
	Your browser does not support the video tag.
	</video>

.. line-block::
	Fig. 4: Free surface elevation simulation (3 seconds) for the radial dam break. Colour-coding: :math:`\xi` (left) and local approximation order (right).



.. [FA2022]  S. 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 <https://doi.org/10.1007/s13137-022-00208-3>`_.