Commit 82c5b420 authored by David Vincent's avatar David Vincent Committed by Jonathan Lambrechts
Browse files

posprocessing file for particle tracking

	--> will be set in slimPost.py once merged from fourier branch
parent 3ae1b336
"""Second-generation Louvain-la-Neuve Ice-ocean model - User interface
Contact: jonathan.lambrechts@uclouvain.be
webpage: www.climate.be/slim
SLIM solves the governing equations on adaptative unstructured meshes using the discontinuous Galerkin finite element method.
slimPost contains some post processing tools usable with results from slim (for simulations, see slim.py file), such as fourier series ...
To generate the post processing results, run the script my_postpro_script.py:
- On a single core: rundgpy my_prepro_script.py
"""
import dgpy
import numpy as np
import matplotlib.pyplot as plt
import time
import os
from dgpy.scripts import slim_private
class Particle_tracker:
""" predict the path of particles from the velocity field predicted by SLIM"""
def __init__(self, mesh_file_name, path_input, bathy, index_start, n_iter, n_period, dt, diffusivity = None) :
"""keyword arguments:
* mesh_file_name
path to the mesh file (.msh format)
* path_input
path to the output files provided by slim (sw2d folder ex: path_input/sw2d/eta/eta.idx)
* bathy
bathymetric file [in meters, positive] (.msh, .idx or .nc format)
* diffusivity
diffusivity of the model (.msh, .idx or .nc format) (default: 0)
* index_start
index of the first iteration
* n_iter
number of iteration
* n_period
number of times the loop of hydrodynamic solution will be run
* dt
time step
"""
self._mesh = mesh_file_name
self._path_input = path_input
self._index_start = index_start
self._n_iter = n_iter
self._nP = n_period
self._dt = dt
self._tot_iter = self._nP * self._n_iter
self._groups = dgpy.dgGroupCollection(self._mesh)
self._bath = dgpy.dgDofContainer(self._groups, 1)
slim_private._load(self._bath, bathy)
self._diffDof = dgpy.dgDofContainer(self._groups,1)
if diffusivity is None:
self._alpha = dgpy.functionConstant(0)
self._diffDof.interpolate(self._alpha)
else:
slim_private._load(self._diffDof, diffusivity)
self._particles = dgpy.particleArray()
self._parts = []
def addParticleAtPoint(self, x, y, number_of_particles, time_step, locationId, status=0, particleId=-1):
"""Seed particle(s) at location (x,y) and time step time_step
keyword arguments:
* x
x coordinates where the particle is seeded
* y
y coordinates where the particle is seeded
* number_of_particles
number of particles seeded at this location
*time_step
index of the time step at which the particle is seeded
*status
status of the particle (dead or alive, default: 0)
*LocationId
Identification of the source
*particleId
Identification of the particle (default: each particle is identified depending on the source)
"""
self._parts.append([x, y, number_of_particles, time_step, status, locationId, particleId])
def loop(self, path_output):
""" Return the fourier serie at a node of the mesh
keyword arguments:
* path_output
name of the .pos files which will contain the position of the particle(s)
"""
self._path_output = path_output
if(not os.path.exists(path_output)):
try : os.mkdir(path_output)
except: pass
self._hydro_sol_dof = dgpy.dgDofContainer(self._groups,3)
self._dof_eta = dgpy.dgDofContainer(self._groups,1)
self._dof_u = dgpy.dgDofContainer(self._groups,1)
self._dof_v = dgpy.dgDofContainer(self._groups,1)
self._dof_eta.importIdx(self._path_input+'/sw2d/eta/eta-%06d.idx'%self._index_start,{0:0})
self._dof_eta.scatter()
self._dof_u.importIdx(self._path_input+'/sw2d/uv/uv-%06d.idx'%self._index_start,{0:0})
self._dof_u.scatter()
self._dof_v.importIdx(self._path_input+'/sw2d/uv/uv-%06d.idx'%self._index_start,{0:1})
self._dof_v.scatter()
self._hydrosolF = slim_private._merger_3d(self._dof_eta.getFunction(),self._dof_u.getFunction(),self._dof_v.getFunction()).functor
self._hydro_sol_dof.interpolate(self._hydrosolF)
self._particleTracker = dgpy.dgParticleTracker2D(self._groups, self._particles, self._bath, self._diffDof, self._tot_iter, self._dt, self._path_output)
startcpu = time.clock()
connectivityMatrix = dgpy.fullMatrixDouble(1,1)
for part in self._parts:
if part[3]==0:
self._particles.addParticlesAtPoint(self._groups, part[2], part[0], part[1], part[3], part[4], part[5], part[6])
self._particles.printPositions(self._path_output+'/particlesAlive_%06d'%( 0 ), 'pos', 0)
for i in range(self._index_start, self._index_start+self._tot_iter+1):
for part in self._parts:
if i>self._index_start and part[3]==(i-self._index_start):
self._particles.addParticlesAtPoint(self._groups, part[2], part[0], part[1], part[3], part[4], part[5], part[6])
print ('LPT iteration:', i, 'of',self._index_start+self._tot_iter+1,'CPU time:', time.clock()-startcpu)
_iter_read = self._index_start +(i-self._index_start)%self._nP
self._dof_eta.importIdx(self._path_input+'/sw2d/eta/eta-%06d.idx'%_iter_read,{0:0})
self._dof_eta.scatter()
self._dof_u.importIdx(self._path_input+'/sw2d/uv/uv-%06d.idx'%_iter_read,{0:0})
self._dof_u.scatter()
self._dof_v.importIdx(self._path_input+'/sw2d/uv/uv-%06d.idx'%_iter_read,{0:1})
self._dof_v.scatter()
self._hydrosolF = slim_private._merger_3d(self._dof_eta.getFunction(),self._dof_u.getFunction(),self._dof_v.getFunction()).functor
self._hydro_sol_dof.interpolate(self._hydrosolF)
self._particleTracker.particleMove(self._hydro_sol_dof, connectivityMatrix, i)
self._particles.printPositions(self._path_output+'/particlesAlive_%06d'%( i ), 'pos', 0)
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment