Commit aeb1ad59 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

recover christopher's testcases

parent 12552a59
......@@ -3,5 +3,5 @@ set(SRC
)
dg_add_module(slimParticle "${SRC}" "")
#dg_add_swig_module(wave wave.i dgWave)
#dg_add_test_directory(tests axel.modave@gmail.com)
dg_add_swig_module(slimParticle slimParticle.i slimParticle)
dg_add_test_directory(tests jonathan.lambrechts@uclouvain.be)
......@@ -564,7 +564,7 @@ void dgParticleTracker2D::particleMove(dgDofContainer* solution, fullMatrix<doub
int index=0;
int indexMax = static_cast<int>(_dt/4.0);
while (index<indexMax) //Safety valve. WARNING: indexMax is arbitrary, may need changing for small-element meshes. Ok for 200m.
while (index<indexMax*10) //Safety valve. WARNING: indexMax is arbitrary, may need changing for small-element meshes. Ok for 200m.
{
index++;
//Get uvw (NB: need to have new uvw at start of each loop cycle)
......@@ -578,6 +578,12 @@ void dgParticleTracker2D::particleMove(dgDofContainer* solution, fullMatrix<doub
for (int i=0; i<3; i++) { uvw[i] += duvw[i]; }
double remaining = 0.0;
int edge_id=-1;
if (1-(uvw[0] + uvw[1]) < 2e-8) {
uvw[0] -= 2e-8;
uvw[1] -= 2e-8;
}
if (uvw[0] < 1e-8) uvw[0] = 2e-8;
if (uvw[1] < 1e-8) uvw[1] = 2e-8;
if((uvw[1] < 0.) && ((uvw[1] / duvw[1]) > remaining)) {
edge_id = 0;
remaining = uvw[1] / duvw[1];
......
%module slimParticle
%{
#undef HAVE_DLOPEN
#include "dgParticleTracker2D.h"
%}
%import(module="dgpy.dgCommon") "dgCommon.i"
%include "dgParticleTracker2D.h"
l = 5e5;
Point(1) = {-l, -l, 0};
Point(2) = {l, -l, 0};
Point(3) = {l, l, 0};
Point(4) = {-l, l, 0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(5) = {3, 4, 1, 2};
Field[1] = MathEval;
Field[1].F = "1e5";
Background Field = 1;
//Field[1] = MathEval;
//Field[1].F = "0.5e4+1*(abs(x)+abs(y))";
//Background Field = 1;
Plane Surface(6) = {5};
Physical Surface("Surface")={6};
Physical Line("Wall") = {1,2,3,4};
Mesh.CharacteristicLengthExtendFromBoundary=1;
Mesh.CharacteristicLengthFromPoints=1;
Mesh.SecondOrderIncomplete=0;
# Lagrangian particle-tracker
# Test for dgParticleTracker2D
from dgpy import *
import time
import os.path
import subprocess
from dgpy.scripts import Common
from dgpy.scripts.termcolor import colored
import numpy as np
print ('')
print(colored('Lagrangian Advection-Diffusion particle tracking module for SLIM. Test #1: basic', "yellow"))
print ('')
#---------- SET OUTPUT FOLDER --------------------------------------------------
LP_OutputDir = 'output'
if(not os.path.exists(LP_OutputDir)):
try : os.mkdir(LP_OutputDir)
except: pass
#-------------------------------------------------------------------------------
# -------- SET HYDRO SIMULATION PARAMETERS -------------------------------------
print(colored('Loading simulation parameters ...',"red"))
simInfoFile = open('simulationInfo.txt','r') # File with info on hydro simulation
simInfo = simInfoFile.readlines()
simDt = float(simInfo[0]) # Simulation dt
simTotIter = int(simInfo[1]) # Tot. number of iters
simExportGap = int(simInfo[2]) # Gap (in iters) between exports
simStartTime = int(float(simInfo[3])) # Simulation start time
simMesh = simInfo[5].strip() # Mesh (not used here)
simInfoFile.close()
#-------------------------------------------------------------------------------
#--------- SET LPT TIME PARAMETERS ---------------------------------------------
print(colored('Loading particle-tracking parameters ...',"red"))
LP_dt = 90. #(double) in seconds
LP_TotTime = 1.0 #(double) in hours
LP_Export = int( (1.0*(3600.0-90.0))/LP_dt ) #(int) export particle positions every ... steps
LP_StartTime = simStartTime
LP_TotIter = int(LP_TotTime*3600.0/LP_dt)
LP_EndTime = LP_StartTime + LP_TotTime*3600.0
simItersPerLPIter = LP_dt/simDt
#-------------------------------------------------------------------------------
#---------- GENERATE & LOAD MESH & BATHYMETRY ----------------------------------
print(colored('Generating and loading mesh and bathymetry ...',"red"))
# The mesh "stommel_square.msh" is generated using the geometry file "stommel_square.geo"
# The generated mesh is 2d (second argument) and of order 1 (third argument)
Common.genMesh('stommel_square', 2, 1)
# Polynomial order of the discretization, and number of dimensions
order = 1
# Groups are used to differentiate different elements of the mesh
# Surface elements and edges belong to different groups
# Group can be split into sub-groups (i.e. split groups of elements that belong to searate physical entities)
groups = dgGroupCollection("stommel_square.msh", order)
# Bathymetry
bathFC = functionConstant(20.0)
bathDC = dgDofContainer(groups, 1)
bathDC.interpolate(bathFC)
#-------------------------------------------------------------------------------
#------- SET UP DIFFUSIVITY ----------------------------------------------------
print(colored('Defining diffusivity ...',"red"))
def diff(cmap, val, alpha):
maxEdges = dgConservationLawFunction.maxEdge(cmap)
val[:] = alpha * 0.00020551 * np.power(maxEdges[:, 0], 1.15)
alpha = functionConstant(1.0)
diffFC = functionNumpy(1, diff, [alpha])
diffDC = dgDofContainer(groups,1)
# Smooth diffusivity:
sys = linearSystemPETScDouble()
dofMan = dgDofManager.newCG (groups, 1, sys)
diffProjector = L2ProjectionContinuous(dofMan)
diffProjector.apply(diffDC,diffFC)
#-------------------------------------------------------------------------------
#--------- SEED PARTICLES ------------------------------------------------------
print(colored('Seeding particles ...',"red"))
#Set max width of domain to seed one group of particles close to walls
LengthOfDomain = 5e5
seedPoints = [[0, 0], [-100000,-100000], [-100000,100000], [100000,-100000], [100000,100000], [LengthOfDomain-1, LengthOfDomain-1] ]
particles = particleArray()
for location in range(0,len(seedPoints)):
particles.addParticlesAtPoint(groups, 4000, seedPoints[location][0], seedPoints[location][1], 0, 0, location)
particles.printPositions('seededParticles','pos',0)
particles.printPositions('seededParticles','dat',0)
#-------------------------------------------------------------------------------
#-------- PRINT SIMULATION INFO ------------------------------------------------
print ('')
print ('-----------------------')
print ('HYDRODYNAMICS PARAMETERS:')
print ('Dt:', simDt,'s')
print ('Simulation length:', simTotIter, '(iter) |', simTotIter*simDt/3600.0, '(h)' )
print ('Simulation starts:', simStartTime/3600.0, '(hr since origin) | Ends:', (simStartTime+simTotIter*simDt)/3600.0, '(hr since origin)' )
print ('Data exported every', simExportGap, '(iter) |', simExportGap*simDt/60.0, '(mins)')
print ('')
print ('-----------------------')
print ('PARTICLE-TRACKING PARAMETERS:')
print ('Dt:', LP_dt, 's')
print ('Simulation length:', LP_TotIter, '(iter) |', LP_TotTime, '(h)')
print ('Simulation starts:', LP_StartTime/3600.0, '(hr since origin) | Ends:', LP_EndTime/3600.0, '(hr since origin)' )
print ('Number of particles:', particles.printNbParticles())
print ('-----------------------')
print ('')
LPsim_TimeOffset = LP_StartTime - simStartTime
print ('-----------------------')
print ('Particle-tracker / Simulation time offset is:', LPsim_TimeOffset/(24.0*3600.0),'(days)')
print ('-----------------------')
print ('')
print ('Output directory:', LP_OutputDir)
print ('')
#-------------------------------------------------------------------------------
#------- SET UP HYDRO DOFCONTAINERS --------------------------------------------
simSolution = dgDofContainer(groups, 3)
simSolution.setAll(2.0)
#-------------------------------------------------------------------------------
#------ SET UP PARTICLE TRACKER OBJECT -----------------------------------------
#Create particle tracker object
print ('Initialising particle tracker ...')
particleTracker = dgParticleTracker2D(groups, particles, bathDC, diffDC, LP_TotIter, LP_dt, LP_OutputDir)
#-------------------------------------------------------------------------------
#------ SET UP DUMMY CONNECTIVITY MATRIX ---------------------------------------
# Required argument of particleMove: NxN matrix
connectivityMatrix = fullMatrixDouble(1,1)
#-------------------------------------------------------------------------------
#-------- PARTICLE LOOP --------------------------------------------------------
print ('Beginning loop ...')
startcpu = time.clock()
for n in range(1,LP_TotIter+1):
t = simStartTime + LPsim_TimeOffset + float(n)*LP_dt
iterNumberWanted = int(LPsim_TimeOffset/simDt + n*simItersPerLPIter)
print ('')
print ('LPT iteration:', n, 'of',LP_TotIter,'| simIteration wanted:',iterNumberWanted, 'for t = %.2f'%( (t-simStartTime)/3600.0 )+ ' (hrs) | CPU time:', time.clock()-startcpu)
#Always give it same solution
particleTracker.particleMove(simSolution, connectivityMatrix, n)
# Print output:
if ( n%(LP_Export) == 0 ):
# Export .pos with positions
particles.printPositions(LP_OutputDir+'/particlesAlive_%06d'%( n ), 'pos', 0)
# Export .dat with positions
particles.printPositions(LP_OutputDir+'/particleOutput_%06d'%( n ), 'dat', 0)
#-------------------------------------------------------------------------------
print(colored('Particle tracking finished. Closing module.',"yellow"))
Msg.Exit(0)
l = 5e5;
Point(1) = {-l, -l, 0};
Point(2) = {l, -l, 0};
Point(3) = {l, l, 0};
Point(4) = {-l, l, 0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(5) = {3, 4, 1, 2};
Field[1] = MathEval;
Field[1].F = "1e5";
Background Field = 1;
//Field[1] = MathEval;
//Field[1].F = "0.5e4+1*(abs(x)+abs(y))";
//Background Field = 1;
Plane Surface(6) = {5};
Physical Surface("Surface")={6};
Physical Line("Wall") = {1,2,3,4};
Mesh.CharacteristicLengthExtendFromBoundary=1;
Mesh.CharacteristicLengthFromPoints=1;
Mesh.SecondOrderIncomplete=0;
# Lagrangian particle-tracker
# Test for dgParticleTracker2D
# This test = basicTest PLUS seeds at location given by a dof container rather than at specified points
from dgpy import *
import time
import os.path
from dgpy.scripts import Common
from termcolor import colored
import numpy as np
print ('')
print(colored('Lagrangian Advection-Diffusion particle tracking module for SLIM. Test #2: basic + seed with dof container', "yellow"))
print ('')
#---------- SET OUTPUT FOLDER --------------------------------------------------
LP_OutputDir = 'output'
if(not os.path.exists(LP_OutputDir)):
try : os.mkdir(LP_OutputDir)
except: 0
#-------------------------------------------------------------------------------
# -------- SET HYDRO SIMULATION PARAMETERS -------------------------------------
print(colored('Loading simulation parameters ...',"red"))
simInfoFile = open('simulationInfo.txt','r') # File with info on hydro simulation
simInfo = simInfoFile.readlines()
simDt = float(simInfo[0]) # Simulation dt
simTotIter = int(simInfo[1]) # Tot. number of iters
simExportGap = int(simInfo[2]) # Gap (in iters) between exports
simStartTime = int(float(simInfo[3])) # Simulation start time
simMesh = simInfo[5].strip() # Mesh (not used here)
simInfoFile.close()
#-------------------------------------------------------------------------------
#--------- SET LPT TIME PARAMETERS ---------------------------------------------
print(colored('Loading particle-tracking parameters ...',"red"))
LP_dt = 90. #(double) in seconds
LP_TotTime = 1.0 #(double) in hours
LP_Export = int( (1.0*(3600.0-90.0))/LP_dt ) #(int) export particle positions every ... steps
LP_StartTime = simStartTime
LP_TotIter = int(LP_TotTime*3600.0/LP_dt)
LP_EndTime = LP_StartTime + LP_TotTime*3600.0
simItersPerLPIter = LP_dt/simDt
#-------------------------------------------------------------------------------
#---------- GENERATE & LOAD MESH & BATHYMETRY ----------------------------------
print(colored('Generating and loading mesh and bathymetry ...',"red"))
# The mesh "stommel_square.msh" is generated using the geometry file "stommel_square.geo"
# The generated mesh is 2d (second argument) and of order 1 (third argument)
Common.genMesh('stommel_square', 2, 1)
# Polynomial order of the discretization, and number of dimensions
# Groups are used to differentiate different elements of the mesh
# Surface elements and edges belong to different groups
# Group can be split into sub-groups (i.e. split groups of elements that belong to searate physical entities)
groups = dgGroupCollection("stommel_square.msh")
# Bathymetry
bathFC = functionConstant(20.0)
bathDC = dgDofContainer(groups, 1)
bathDC.interpolate(bathFC)
#-------------------------------------------------------------------------------
#------- SET UP DIFFUSIVITY ----------------------------------------------------
print(colored('Defining diffusivity ...',"red"))
def diff(cmap, val, alpha):
maxEdges = dgConservationLawFunction.maxEdge(cmap)[:, 0]
val[:] = alpha * 0.00020551 * np.power(maxEdges, 1.15)
alpha = functionConstant(1.0)
diffFC = functionNumpy(1, diff, [alpha])
diffDC = dgDofContainer(groups,1)
# Smooth diffusivity:
sys = linearSystemPETScDouble()
dofMan = dgDofManager.newCG (groups, 1, sys)
diffProjector = L2ProjectionContinuous(dofMan)
diffProjector.apply(diffDC,diffFC)
#-------------------------------------------------------------------------------
#--------- SEED PARTICLES ------------------------------------------------------
print(colored('Seeding particles ...',"red"))
#Define function where we want to seed particles
XYZ = function.getCoordinates()
def whereToSeed(cmap, val, xyz):
val[:] = np.where(xyz[:,0]>0, 1, 0)
whereToSeed_FN = functionNumpy(1, whereToSeed, [XYZ])
whereToSeed_DC = dgDofContainer(groups, 1)
whereToSeed_DC.interpolate(whereToSeed_FN)
whereToSeed_DC.exportMsh("whereToSeed")
# Seed using function
particles = particleArray()
particles.addParticlesWithDofContainer(groups, whereToSeed_DC, 0.05, 1000, 0)
particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
particles.printPositions(LP_OutputDir+'/seededParticles','dat',0)
#-------------------------------------------------------------------------------
#-------- PRINT SIMULATION INFO ------------------------------------------------
print ('')
print ('-----------------------')
print ('HYDRODYNAMICS PARAMETERS:')
print ('Dt:', simDt,'s')
print ('Simulation length:', simTotIter, '(iter) |', simTotIter*simDt/3600.0, '(h)' )
print ('Simulation starts:', simStartTime/3600.0, '(hr since origin) | Ends:', (simStartTime+simTotIter*simDt)/3600.0, '(hr since origin)' )
print ('Data exported every', simExportGap, '(iter) |', simExportGap*simDt/60.0, '(mins)')
print ('')
print ('-----------------------')
print ('PARTICLE-TRACKING PARAMETERS:')
print ('Dt:', LP_dt, 's')
print ('Simulation length:', LP_TotIter, '(iter) |', LP_TotTime, '(h)')
print ('Simulation starts:', LP_StartTime/3600.0, '(hr since origin) | Ends:', LP_EndTime/3600.0, '(hr since origin)' )
print ('Number of particles:', particles.printNbParticles())
print ('-----------------------')
print ('')
LPsim_TimeOffset = LP_StartTime - simStartTime
print ('-----------------------')
print ('Particle-tracker / Simulation time offset is:', LPsim_TimeOffset/(24.0*3600.0),'(days)')
print ('-----------------------')
print ('')
print ('Output directory:', LP_OutputDir)
print ('')
#-------------------------------------------------------------------------------
#------- SET UP HYDRO DOFCONTAINERS --------------------------------------------
simSolution = dgDofContainer(groups, 3)
simSolution.setAll(2.0)
#-------------------------------------------------------------------------------
#------ SET UP PARTICLE TRACKER OBJECT -----------------------------------------
#Create particle tracker object
print ('Initialising particle tracker ...')
particleTracker = dgParticleTracker2D(groups, particles, bathDC, diffDC, LP_TotIter, LP_dt, LP_OutputDir)
#-------------------------------------------------------------------------------
#------ SET UP DUMMY CONNECTIVITY MATRIX ---------------------------------------
# Required argument of particleMove: NxN matrix
connectivityMatrix = fullMatrixDouble(1,1)
#-------------------------------------------------------------------------------
#-------- PARTICLE LOOP --------------------------------------------------------
print ('Beginning loop ...')
startcpu = time.clock()
for n in range(1,LP_TotIter+1):
t = simStartTime + LPsim_TimeOffset + float(n)*LP_dt
iterNumberWanted = int(LPsim_TimeOffset/simDt + n*simItersPerLPIter)
print ('')
print ('LPT iteration:', n, 'of',LP_TotIter,'| simIteration wanted:',iterNumberWanted, 'for t = %.2f'%( (t-simStartTime)/3600.0 )+ ' (hrs) | CPU time:', time.clock()-startcpu)
#Always give it same solution
particleTracker.particleMove(simSolution, connectivityMatrix, n)
# Print output:
if ( n%(LP_Export) == 0 ):
# Export .pos with positions
particles.printPositions(LP_OutputDir+'/particlesAlive_%06d'%( n ), 'pos', 0)
# Export .dat with positions
particles.printPositions(LP_OutputDir+'/particleOutput_%06d'%( n ), 'dat', 0)
#-------------------------------------------------------------------------------
print(colored('Particle tracking finished. Closing module.',"yellow"))
Msg.Exit(0)
l = 5e5;
Point(1) = {-l, -l, 0};
Point(2) = {l, -l, 0};
Point(3) = {l, l, 0};
Point(4) = {-l, l, 0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(5) = {3, 4, 1, 2};
Field[1] = MathEval;
Field[1].F = "1e5";
Background Field = 1;
//Field[1] = MathEval;
//Field[1].F = "0.5e4+1*(abs(x)+abs(y))";
//Background Field = 1;
Plane Surface(6) = {5};
Physical Surface("Surface")={6};
Physical Line("Wall") = {1,2,3,4};
Mesh.CharacteristicLengthExtendFromBoundary=1;
Mesh.CharacteristicLengthFromPoints=1;
Mesh.SecondOrderIncomplete=0;
# Lagrangian particle-tracker
# Test for dgParticleTracker2D
from dgpy import *
import time
import os.path
import subprocess
from dgpy.scripts import Common
from termcolor import colored
import numpy as np
print ('')
print((colored('Lagrangian Advection-Diffusion particle tracking module for SLIM. Test #3: direct wind forcing of particles', "yellow")))
print ('')
#---------- SET OUTPUT FOLDER --------------------------------------------------
LP_OutputDir = 'output'
if(not os.path.exists(LP_OutputDir)):
try : os.mkdir(LP_OutputDir)
except: 0
#-------------------------------------------------------------------------------
# -------- SET HYDRO SIMULATION PARAMETERS -------------------------------------
print((colored('Loading simulation parameters ...',"red")))
simInfoFile = open('simulationInfo.txt','r') # File with info on hydro simulation
simInfo = simInfoFile.readlines()
simDt = float(simInfo[0]) # Simulation dt
simTotIter = int(simInfo[1]) # Tot. number of iters
simExportGap = int(simInfo[2]) # Gap (in iters) between exports
simStartTime = int(float(simInfo[3])) # Simulation start time
simMesh = simInfo[5].strip() # Mesh (not used here)
simInfoFile.close()
#-------------------------------------------------------------------------------
#--------- SET LPT TIME PARAMETERS ---------------------------------------------
print((colored('Loading particle-tracking parameters ...',"red")))
LP_dt = 90. #(double) in seconds
LP_TotTime = 1.0 #(double) in hours
LP_Export = int( (1.0*(3600.0-90.0))/LP_dt ) #(int) export particle positions every ... steps
LP_StartTime = simStartTime
LP_TotIter = int(LP_TotTime*3600.0/LP_dt)
LP_EndTime = LP_StartTime + LP_TotTime*3600.0
simItersPerLPIter = LP_dt/simDt
#-------------------------------------------------------------------------------
#---------- GENERATE & LOAD MESH & BATHYMETRY ----------------------------------
print((colored('Generating and loading mesh and bathymetry ...',"red")))
# The mesh "stommel_square.msh" is generated using the geometry file "stommel_square.geo"
# The generated mesh is 2d (second argument) and of order 1 (third argument)
Common.genMesh('stommel_square', 2, 1)
# Polynomial order of the discretization, and number of dimensions
# Groups are used to differentiate different elements of the mesh
# Surface elements and edges belong to different groups
# Group can be split into sub-groups (i.e. split groups of elements that belong to searate physical entities)
groups = dgGroupCollection("stommel_square.msh")
# Bathymetry
bathFC = functionConstant(20.0)
bathDC = dgDofContainer(groups, 1)
bathDC.interpolate(bathFC)
#-------------------------------------------------------------------------------
#------- SET UP DIFFUSIVITY ----------------------------------------------------
print((colored('Defining diffusivity ...',"red")))
def diff(cmap, val, alpha):
maxEdges = dgConservationLawFunction.maxEdge(cmap)[:, 0]
val[:] = alpha * 0.00020551 * np.power(maxEdges, 1.15)
alpha = functionConstant(1.0)
diffFC = functionNumpy(1, diff, [alpha])
diffDC = dgDofContainer(groups,1)
# Smooth diffusivity:
sys = linearSystemPETScDouble()
dofMan = dgDofManager.newCG (groups, 1, sys)
diffProjector = L2ProjectionContinuous(dofMan)
diffProjector.apply(diffDC,diffFC)
#-------------------------------------------------------------------------------
#--------- SEED PARTICLES ------------------------------------------------------
print((colored('Seeding particles ...',"red")))
#Set max width of domain to seed one group of particles close to walls
LengthOfDomain = 5e5
seedPoints = [[0, 0], [-100000,-100000], [-100000,100000], [100000,-100000], [100000,100000], [LengthOfDomain-1, LengthOfDomain-1] ]
particles = particleArray()
for location in range(0,len(seedPoints)):
particles.addParticlesAtPoint(groups, 4000, seedPoints[location][0], seedPoints[location][1], 0, 0, location)
particles.printPositions('seededParticles','pos',0)
particles.printPositions('seededParticles','dat',0)
#-------------------------------------------------------------------------------
#-------- PRINT SIMULATION INFO ------------------------------------------------
print('')
print('-----------------------')
print('HYDRODYNAMICS PARAMETERS:')
print('Dt:', simDt,'s')
print('Simulation length:', simTotIter, '(iter) |', simTotIter*simDt/3600.0, '(h)')
print('Simulation starts:', simStartTime/3600.0, '(hr since origin) | Ends:', (simStartTime+simTotIter*simDt)/3600.0, '(hr since origin)')
print('Data exported every', simExportGap, '(iter) |', simExportGap*simDt/60.0, '(mins)')
print('')
print('-----------------------')
print('PARTICLE-TRACKING PARAMETERS:')
print('Dt:', LP_dt, 's')
print('Simulation length:', LP_TotIter, '(iter) |', LP_TotTime, '(h)')
print('Simulation starts:', LP_StartTime/3600.0, '(hr since origin) | Ends:', LP_EndTime/3600.0, '(hr since origin)')
print('Number of particles:', particles.printNbParticles())
print('-----------------------')
print('')
LPsim_TimeOffset = LP_StartTime - simStartTime
print('-----------------------')
print('Particle-tracker / Simulation time offset is:', LPsim_TimeOffset/(24.0*3600.0),'(days)')
print('-----------------------')