Tracer source

In this example, a pollutant is released and convected in a river. The flow in the river is subcritical and boundary conditions are fluvial flowrate for inlet and fluvial height for outlet. Pollutant is defined as a tracer and its release into the river is modelled with fk.TracerSource external effect i.e. a quantity of tracer is realeased at given times in designated cells.

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

os.system('rm -r outputs')

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

Time loop:

FINAL_TIME = 50.

simutime = fk.SimuTime(final_time=FINAL_TIME, time_iteration_max=10000,
                       second_order=False)

Mesh:

dir_path = os.path.abspath(os.path.dirname(sys.argv[0]))
triangular_mesh = fk.TriangularMesh.from_msh_file(dir_path + '/inputs/river_polluant.msh')

if not args.nographics:
  fk_plt.plot_mesh(triangular_mesh, plot_labels=False)

Layers:

NL = 4
layer = fk.Layer(NL, triangular_mesh, topography=0.)

Primitives:

primitives = fk.Primitive(triangular_mesh, layer, free_surface=2.)

Tracer:

Important

Let us define two tracers, labeled polluant_0 and polluant_1. Labels are very important to link boundary conditions as well as source terms to the correct tracer in the code.

tracers = [fk.Tracer(triangular_mesh, layer, primitives, Tinit=0., label='polluant_0'),
           fk.Tracer(triangular_mesh, layer, primitives, Tinit=10., label='polluant_1')]

Boundary conditions:

Boundary conditions are the same as in reiver example excepts we define inlet flowrate in a input file.

times, flowrates = np.loadtxt('inputs/river_polluant_flowrate.txt', unpack=True)
flowrate = fk.TimeDependentFlowRate(times=times, flowrates=flowrates)

fluvial_flowrates = [fk.FluvialFlowRate(
    ref=1, time_dependent_flowrate=flowrate, x_flux_direction=1.0, y_flux_direction=0.0,
    tracers_inlet={'polluant_0':0., 'polluant_1':10.})]

fluvial_heights = [fk.FluvialHeight(
    ref=2, height=2., tracers_inlet={'polluant_0':0., 'polluant_1':10.})]

slides = [fk.Slide(ref=3)]

Note

Note that tracer labels have been used to provide tracer value at inlet boundary condition. In our case, we provide the same boundary condition inlet values as for tracer initialization.

External effects (tracer source terms):

ExternalEffects objects are called in problem.solve() at each time step. In this example external effects are TracerSource. There can be as much external effects as needed as long as they are stored in a dictionary.

TracerSource can have different unit based on the way you want it to be applied.

• If source_unit='[T]/s': tracer source is continuous over time, and is applied at each time step. In this case tracer source value is interpolated linearly from sources and times parameters.

• If source_unit='[T]': tracer source is sparse over time, and is applied only on specified times with the corresponding value of sources.

polluant_source_0 = fk.TracerSource(
    times=[1., FINAL_TIME],
    sources=[20, 20],
    source_unit='[T]/s',
    tracer_label='polluant_0',
    cells=[607])

polluant_source_1 = fk.TracerSource(
    times=[1, 10, 20, 30],
    sources=[20, 20, 20, 20],
    source_unit='[T]',
    tracer_label='polluant_1',
    cells=[617])

external_effects = {"continuous": polluant_source_0, "sparse": polluant_source_1}

Note

• Once again, tracer labels are used to apply source terms on each one of the tracers.

• Note that tracer sources can be set in kg instead of [T]=kg/m3 in the case of species. The corresponding concentration source term is computed based on cell volume and mass source.

Writter:

NB_VTK = 10
vtk_writer = fk.VTKWriter(triangular_mesh, scheduler=fk.schedules(count=NB_VTK),
                          scale_h=5.)

Problem definition:

problem = fk.Problem(simutime, triangular_mesh, layer, primitives,
                     slides=slides,
                     fluvial_heights=fluvial_heights,
                     fluvial_flowrates=fluvial_flowrates,
                     tracers=tracers,
                     external_effects=external_effects,
                     vtk_writer=vtk_writer)

===================================================================
|                          INITIALIZATION                         |
===================================================================
Problem size: Nodes=866, Layers=4, Triangles=1594,
Iter =        0 ; Dt = 0.0000s ; Time =     0.00s ; ETA =     0.00s

Problem solving:

problem.solve()

===================================================================
|                            TIME LOOP                            |
===================================================================
Iter =       68 ; Dt = 0.0150s ; Time =     1.02s ; ETA =    24.86s
Iter =      136 ; Dt = 0.0150s ; Time =     2.04s ; ETA =    21.10s
Iter =      204 ; Dt = 0.0150s ; Time =     3.07s ; ETA =    20.71s
Iter =      272 ; Dt = 0.0150s ; Time =     4.09s ; ETA =    19.99s
Iter =      340 ; Dt = 0.0150s ; Time =     5.11s ; ETA =    19.43s
Iter =      408 ; Dt = 0.0148s ; Time =     6.12s ; ETA =    19.18s
Iter =      478 ; Dt = 0.0145s ; Time =     7.15s ; ETA =    19.31s
Iter =      549 ; Dt = 0.0142s ; Time =     8.17s ; ETA =    19.29s
Iter =      621 ; Dt = 0.0141s ; Time =     9.18s ; ETA =    19.23s
Iter =      694 ; Dt = 0.0141s ; Time =    10.21s ; ETA =    18.51s
Iter =      766 ; Dt = 0.0141s ; Time =    11.23s ; ETA =    18.07s
Iter =      838 ; Dt = 0.0141s ; Time =    12.24s ; ETA =    17.32s
Iter =      911 ; Dt = 0.0141s ; Time =    13.28s ; ETA =    16.88s
Iter =      983 ; Dt = 0.0141s ; Time =    14.29s ; ETA =    16.42s
Iter =     1055 ; Dt = 0.0141s ; Time =    15.31s ; ETA =    16.11s
Iter =     1127 ; Dt = 0.0141s ; Time =    16.33s ; ETA =    15.57s
Iter =     1200 ; Dt = 0.0141s ; Time =    17.36s ; ETA =    15.22s
Iter =     1273 ; Dt = 0.0140s ; Time =    18.38s ; ETA =    15.01s
Iter =     1345 ; Dt = 0.0140s ; Time =    19.39s ; ETA =    14.45s
Iter =     1418 ; Dt = 0.0140s ; Time =    20.41s ; ETA =    14.60s
Iter =     1492 ; Dt = 0.0139s ; Time =    21.44s ; ETA =    14.05s
Iter =     1565 ; Dt = 0.0139s ; Time =    22.46s ; ETA =    13.41s
Iter =     1638 ; Dt = 0.0139s ; Time =    23.48s ; ETA =    13.26s
Iter =     1711 ; Dt = 0.0139s ; Time =    24.49s ; ETA =    12.55s
Iter =     1784 ; Dt = 0.0139s ; Time =    25.51s ; ETA =    12.01s
Iter =     1858 ; Dt = 0.0139s ; Time =    26.54s ; ETA =    11.77s
Iter =     1931 ; Dt = 0.0139s ; Time =    27.56s ; ETA =    11.00s
Iter =     2004 ; Dt = 0.0140s ; Time =    28.58s ; ETA =    10.65s
Iter =     2077 ; Dt = 0.0140s ; Time =    29.60s ; ETA =    10.02s
Iter =     2149 ; Dt = 0.0141s ; Time =    30.61s ; ETA =     9.31s
Iter =     2222 ; Dt = 0.0142s ; Time =    31.65s ; ETA =     8.95s
Iter =     2293 ; Dt = 0.0143s ; Time =    32.66s ; ETA =     8.20s
Iter =     2364 ; Dt = 0.0144s ; Time =    33.68s ; ETA =     7.37s
Iter =     2435 ; Dt = 0.0145s ; Time =    34.71s ; ETA =     6.93s
Iter =     2505 ; Dt = 0.0146s ; Time =    35.72s ; ETA =     6.85s
Iter =     2575 ; Dt = 0.0146s ; Time =    36.75s ; ETA =     6.16s
Iter =     2644 ; Dt = 0.0146s ; Time =    37.76s ; ETA =     5.78s
Iter =     2714 ; Dt = 0.0147s ; Time =    38.78s ; ETA =     5.41s
Iter =     2784 ; Dt = 0.0147s ; Time =    39.81s ; ETA =     4.91s
Iter =     2853 ; Dt = 0.0147s ; Time =    40.82s ; ETA =     4.36s
Iter =     2923 ; Dt = 0.0147s ; Time =    41.85s ; ETA =     3.89s
Iter =     2992 ; Dt = 0.0147s ; Time =    42.87s ; ETA =     3.16s
Iter =     3061 ; Dt = 0.0147s ; Time =    43.88s ; ETA =     2.87s
Iter =     3130 ; Dt = 0.0148s ; Time =    44.90s ; ETA =     2.40s
Iter =     3199 ; Dt = 0.0148s ; Time =    45.92s ; ETA =     1.93s
Iter =     3268 ; Dt = 0.0148s ; Time =    46.94s ; ETA =     1.41s
Iter =     3337 ; Dt = 0.0148s ; Time =    47.96s ; ETA =     0.95s
Iter =     3406 ; Dt = 0.0148s ; Time =    48.99s ; ETA =     0.48s
Iter =     3475 ; Dt = 0.0029s ; Time =    50.00s ; ETA =     0.00s
===================================================================
|                               END                               |
===================================================================
Problem.solve() completed in 24.927244663238525s (wall time)

Basic plots:

if not args.nographics:

  x = np.asarray(triangular_mesh.triangulation.x)
  y = np.asarray(triangular_mesh.triangulation.y)
  trivtx = np.asarray(triangular_mesh.triangulation.trivtx)
  triang = mtri.Triangulation(x, y, trivtx)

  LAYER_IDX = [int(NL/3), int(2*NL/3)]

  plt.rcParams["figure.figsize"] = [8, 10]
  fig = plt.figure()

  ax1 = fig.add_subplot(211)
  T0 = problem.tracers[0].T[:, LAYER_IDX[1]]
  v = np.linspace(min(T0), max(T0), 20, endpoint=True)
  ax1.triplot(x, y, trivtx, color='k', lw=0.5)
  im1 = ax1.tricontourf(x, y, trivtx, T0, v, cmap=plt.cm.jet)
  fig.colorbar(im1, ax=ax1, ticks=v)
  ax1.set_title("Tracer 0 on layer {}".format(LAYER_IDX[1]))

  ax1b = fig.add_subplot(212)
  T1 = problem.tracers[1].T[:, LAYER_IDX[1]]
  v = np.linspace(min(T1), max(T1), 20, endpoint=True)
  ax1b.triplot(x, y, trivtx, color='k', lw=0.5)
  im1b = ax1b.tricontourf(x, y, trivtx, T1, v, cmap=plt.cm.jet)
  fig.colorbar(im1b, ax=ax1b, ticks=v)
  ax1b.set_title("Tracer 1 on layer {}".format(LAYER_IDX[1]))
  plt.show()