recover christopher's testcases

modules/slimParticles/CMakeLists.txt

@@ -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)

modules/slimParticles/dgParticleTracker2D.cpp

@@ -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];

modules/slimParticles/slimParticle.i

+%module slimParticle
+
+%{
+  #undef HAVE_DLOPEN
+  #include "dgParticleTracker2D.h"
+%}
+%import(module="dgpy.dgCommon") "dgCommon.i"
+
+%include "dgParticleTracker2D.h"

modules/slimParticles/tests/parTracker_basic/stommel_square.geo

+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;
+

modules/slimParticles/tests/parTracker_basic/test.py

+# 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)

modules/slimParticles/tests/parTracker_seedWithDof/stommel_square.geo

+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;
+

modules/slimParticles/tests/parTracker_seedWithDof/test.py

+# 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)

modules/slimParticles/tests/parTracker_wind/stommel_square.geo

+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;
+

modules/slimParticles/tests/parTracker_wind/test.py

+# 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('-----------------------')
+print('')
+print('Output directory:', LP_OutputDir)
+print('')
+#-------------------------------------------------------------------------------
+
+#------- SET UP HYDRO DOFCONTAINERS --------------------------------------------
+simSolution = dgDofContainer(groups, 3)
+simSolution.setAll(2.0)