Thacker 2d

In this example, 2D thacker simulation is carried out. It consists of a planar surface rotating in a paraboloid without friction. The free surface motion is periodic in time and remains planar. One could think of water rotating in a paraboloïd glass.

Note

Analytical solution is given in: https://hal.archives-ouvertes.fr/hal-00628246/document at section 4.2.2, paragraph “Planar surface in a paraboloid”.

import os, sys
import argparse
import matplotlib.pyplot as plt
import numpy as np
import freshkiss3d as fk
import freshkiss3d.extra.plots as fk_plt

#sphinx_gallery_thumbnail_number = 3

parser = argparse.ArgumentParser()
parser.add_argument('--nographics', action='store_true')
args = parser.parse_args()

Parameters:

A = 1.
H0 = 0.1
ETA = 0.5
L = 4.

#Numerical scheme:
IPRES = True
SECOND_ORDER = True
NL = 2

Time loop:

simutime = fk.SimuTime(final_time=4.5, time_iteration_max=10000,
                       second_order=SECOND_ORDER)

when = [0., 2.25, 4.5]
create_figure_scheduler = fk.schedules(times=when)

Mesh:

dir_path = os.path.abspath(os.path.dirname(sys.argv[0]))
TG, vertex_labels, boundary_labels = fk.read_msh(dir_path + '/inputs/thacker2d.msh')
x = np.asarray(TG.x)
y = np.asarray(TG.y)
trivtx = np.asarray(TG.trivtx)

x *= 0.4
y *= 0.4

triangular_mesh = fk.TriangularMesh(TG, vertex_labels, boundary_labels)

if not args.nographics:
  fk_plt.plot_mesh(triangular_mesh)

Analytic solution:

Analytic solution is defined in a Thacker2D class with h0 being the water depth at the central point of the domain for a zero elevation and a is the distance from this central point to the zero elevation of the shoreline.

thacker2d_analytic = fk.Thacker2D(triangular_mesh, a=A, h0=H0, eta=ETA, L=L,
                                  scheduler=fk.schedules(times=when), NL=NL)

The topography is defined in freshkiss3d.extra.analyticsol.Thacker2D so that:

\[z(r) = -h_0 \left( 1 - \frac{r^2}{a^2} \right)\]

\[r = \sqrt{(x-L/2)^2 + (y-L/2)^2}\]

thacker2d_analytic(0.)

Analytic solution is computed with the following formula:

\[\begin{split}\begin{cases} h(x,y,t) = \frac{\eta h_0}{a^2} \left( 2\left( x- \frac{L}{2} \right) \cos (\omega t) + 2\left( x- \frac{L}{2} \right) \sin (\omega t) - \eta \right) - z(x,y) \\ u(x,y,t) = - \eta \omega \sin( \omega t ) \\ v(x,y,t) = \eta \omega \cos( \omega t ) \end{cases}\end{split}\]

with:

\[\omega = \sqrt{2g h_0}/a \quad \text{and} \quad \eta = 0.5\]

Warning

Analytic solution is available in Thacker2d for single layer St Venant only

Layers:

Note

Topography defined in thacker2d is shared with layer

layer = fk.Layer(NL, triangular_mesh, topography=thacker2d_analytic.topography)

Primitives:

Note

H, U and V are initialized with analytic initial solution

primitives = fk.Primitive(triangular_mesh, layer,
                          height=thacker2d_analytic.H,
                          Uinit=thacker2d_analytic.U,
                          Vinit=thacker2d_analytic.V)

Topography and initial height computed in thacker2d_analytic are the following:

def Topo(x, y):
    r2 = (x - 0.5*L)**2 + (y-0.5*L)**2
    return - H0 * (1 - (r2/A**2))

def H_0(x, y):
    return ETA * H0 / A**2 * (2* (x-0.5*L) - ETA) - Topo(x, y)

if not args.nographics:
  fk_plt.plot_init_1d_slice(triangular_mesh, surf_xy=H_0, topo_xy=Topo, surface_type='height')

Boundary conditions:

fluvial_heights = [fk.FluvialHeight(ref=r, height=0.0) for r in [1, 2, 3, 4]]

Problem definition:

problem = fk.Problem(simutime, triangular_mesh, layer, primitives,
                     fluvial_heights=fluvial_heights,
                     analytic_sol=thacker2d_analytic,
                     numerical_parameters={'ipres':IPRES,'space_second_order':SECOND_ORDER},
                     custom_funct={'plot':fk_plt.plot_freesurface_3d_analytic_2},
                     custom_funct_scheduler=create_figure_scheduler)

===================================================================
|                          INITIALIZATION                         |
===================================================================
Problem size: Nodes=2601, Layers=2, Triangles=5000,
Iter =        0 ; Dt = 0.0000s ; Time =     0.00s ; ETA =     0.00s

Problem solving:

When a figure plot is scheduled thacker2d.compute is called to compute analytic solution.

problem.solve()

if not args.nographics:
  plt.tight_layout()
  plt.show()

• 
• 

===================================================================
|                            TIME LOOP                            |
===================================================================
===================================================================
|                      L2 errors at time 0 s:
|                           H : 0.000272
|                          QX : 0.000107
|                          QY : 0.008785
===================================================================
Iter =       11 ; Dt = 0.0086s ; Time =     0.10s ; ETA =     3.71s
Iter =       22 ; Dt = 0.0085s ; Time =     0.19s ; ETA =     3.66s
Iter =       33 ; Dt = 0.0084s ; Time =     0.28s ; ETA =     6.32s
Iter =       44 ; Dt = 0.0084s ; Time =     0.37s ; ETA =     3.69s
Iter =       55 ; Dt = 0.0085s ; Time =     0.47s ; ETA =     3.59s
Iter =       65 ; Dt = 0.0086s ; Time =     0.55s ; ETA =     3.58s
Iter =       76 ; Dt = 0.0085s ; Time =     0.65s ; ETA =     3.54s
Iter =       87 ; Dt = 0.0084s ; Time =     0.74s ; ETA =     3.58s
Iter =       98 ; Dt = 0.0084s ; Time =     0.83s ; ETA =     3.52s
Iter =      109 ; Dt = 0.0081s ; Time =     0.92s ; ETA =     3.58s
Iter =      121 ; Dt = 0.0075s ; Time =     1.02s ; ETA =     3.87s
Iter =      133 ; Dt = 0.0071s ; Time =     1.10s ; ETA =     3.96s
Iter =      146 ; Dt = 0.0073s ; Time =     1.20s ; ETA =     5.01s
Iter =      159 ; Dt = 0.0073s ; Time =     1.29s ; ETA =     3.76s
Iter =      171 ; Dt = 0.0073s ; Time =     1.38s ; ETA =     3.77s
Iter =      185 ; Dt = 0.0068s ; Time =     1.47s ; ETA =     3.88s
Iter =      199 ; Dt = 0.0061s ; Time =     1.57s ; ETA =     4.30s
Iter =      213 ; Dt = 0.0063s ; Time =     1.65s ; ETA =     4.07s
Iter =      228 ; Dt = 0.0069s ; Time =     1.75s ; ETA =     3.65s
Iter =      242 ; Dt = 0.0065s ; Time =     1.84s ; ETA =     3.79s
Iter =      255 ; Dt = 0.0073s ; Time =     1.93s ; ETA =     3.20s
Iter =      268 ; Dt = 0.0069s ; Time =     2.02s ; ETA =     3.31s
Iter =      281 ; Dt = 0.0073s ; Time =     2.12s ; ETA =     3.06s
Iter =      293 ; Dt = 0.0079s ; Time =     2.21s ; ETA =     2.74s
===================================================================
|                      L2 errors at time 2 s:
|                           H : 0.001378
|                          QX : 0.000432
|                          QY : 0.009114
===================================================================
Iter =      304 ; Dt = 0.0083s ; Time =     2.30s ; ETA =     2.49s
Iter =      315 ; Dt = 0.0085s ; Time =     2.39s ; ETA =     2.36s
Iter =      326 ; Dt = 0.0085s ; Time =     2.49s ; ETA =     2.26s
Iter =      337 ; Dt = 0.0085s ; Time =     2.57s ; ETA =     2.24s
Iter =      349 ; Dt = 0.0077s ; Time =     2.67s ; ETA =     2.32s
Iter =      361 ; Dt = 0.0064s ; Time =     2.76s ; ETA =     2.68s
Iter =      374 ; Dt = 0.0072s ; Time =     2.85s ; ETA =     2.31s
Iter =      387 ; Dt = 0.0073s ; Time =     2.94s ; ETA =     2.12s
Iter =      400 ; Dt = 0.0067s ; Time =     3.03s ; ETA =     2.22s
Iter =      413 ; Dt = 0.0072s ; Time =     3.13s ; ETA =     1.94s
Iter =      425 ; Dt = 0.0068s ; Time =     3.22s ; ETA =     1.88s
Iter =      439 ; Dt = 0.0066s ; Time =     3.31s ; ETA =     1.82s
Iter =      453 ; Dt = 0.0067s ; Time =     3.40s ; ETA =     1.65s
Iter =      465 ; Dt = 0.0081s ; Time =     3.49s ; ETA =     1.26s
Iter =      476 ; Dt = 0.0086s ; Time =     3.59s ; ETA =     1.09s
Iter =      487 ; Dt = 0.0050s ; Time =     3.67s ; ETA =     1.71s
Iter =      501 ; Dt = 0.0082s ; Time =     3.77s ; ETA =     0.92s
Iter =      514 ; Dt = 0.0078s ; Time =     3.86s ; ETA =     0.84s
Iter =      528 ; Dt = 0.0078s ; Time =     3.95s ; ETA =     0.72s
Iter =      543 ; Dt = 0.0070s ; Time =     4.04s ; ETA =     0.68s
Iter =      556 ; Dt = 0.0077s ; Time =     4.14s ; ETA =     0.49s
Iter =      568 ; Dt = 0.0069s ; Time =     4.23s ; ETA =     0.40s
Iter =      581 ; Dt = 0.0082s ; Time =     4.32s ; ETA =     0.22s
Iter =      592 ; Dt = 0.0074s ; Time =     4.41s ; ETA =     0.12s
Iter =      604 ; Dt = 0.0014s ; Time =     4.50s ; ETA =     0.00s
===================================================================
|                               END                               |
===================================================================
Problem.solve() completed in 5.828488111495972s (wall time)