LPT_Caller_DualReefMap.py

# Lagrangian particle-tracker
# Requires simulationInfo.txt file with hydrodynamics simulation parameters: 
#  -> line1: simulation dt, line2: total nb of iters, line3: nb of iters between each export, line4: simulation start time, line5: output directory, line6: mesh name
# Mesh must be in planar coords in Meshes/ folder. Stereo-coords: WIP
# Makefile will download reef map into Data/ folder
# Hydrodynamics should be in Output/ folder
# Bathymetry (in msh format, same mesh as hydrodynamics) should be in Bath/ folder
#
# !!! NB: Must update t_preComp in settleParams, EVEN IF settleType=0 !!!
#      --> Because _t_preComp is defined here, even for exposure time-connectivity method

from dgpy import *
import sys, time
import os
import os.path
from termcolor import colored
from ProjectMesh import *

print ''
print(colored('Lagrangian Advection-Diffusion particle tracking module for SLIM. Working offline.', "yellow"))
print ''

#--------- SETUP LIBRARIES  ----------------------------------------------------
if (Msg.GetCommRank()==0):
  print ''
  print(colored('Compiling Libraries ...', "red"))
  functionC.buildLibraryFromFile ("Libraries/GBR.cc", "Libraries/lib_gbr.so");
#  functionC.buildLibraryFromFile ("Libraries/REEF.cc", "Libraries/lib_reef.so");
  
glib="Libraries/lib_gbr.so"
#rlib="Libraries/lib_reef.so"
#-------------------------------------------------------------------------------

#---------- SET INPUT & OUTPUT FOLDERS -----------------------------------------
#simOutputDir = '/scratch/cthomas/Hydro_Realistic/gbr500K_sma_1svIn_bDf10_allTides_windOn_28112007_5wks/'
simOutputDir = '/home/cthomas/GBR_slim/Saved_Data/gbr460KTests/460K_tOnwCSFRbNCJ0p8SvSth_FWind1o50Every_bDn016Reefs10'

#LP_OutputDir = '/scratch/cthomas/temp'#simOutputDir +'/LPT_Output/SettProg_Mort007'
LP_OutputDir = 'temp'

##
# Add /CPUX to LP_OutputDir if -c option used
cpuExt=""
if ( len(sys.argv)>1 ):
  if sys.argv[1] == "-c":
    cpuExt = "/CPU"+sys.argv[2]
    time.sleep(float(sys.argv[2])*2) #Avoid simultaneous simulations crashing each other

LP_OutputDir = LP_OutputDir + cpuExt
if(not os.path.exists(LP_OutputDir)):
  try : os.mkdir(LP_OutputDir)
  except: 0;
##
#-------------------------------------------------------------------------------


#reefMapShallowFilenameTrunk = "SGBR_ReefMap_2121_Solid9200.dat"
#reefMapDeepFilenameTrunk = "SGBR_DeepReefMap_2121_Solid9200.dat"
reefMapShallowFilenameTrunk = "GBR_ReefMap_2966_Solid9200"
reefMapDeepFilenameTrunk = "GBR_DeepReefMap_2966_Solid9200"
reefMapShallowFilename = reefMapShallowFilenameTrunk+".dat"
reefMapDeepFilename = reefMapDeepFilenameTrunk+".dat"

# Radius: 6.37101e6 for old meshes (sma_zoomed), else 6.371e6 (used for reefsObj and reefsDistanceObj, the former only if useProj=0)
EarthRadius=6.371e6
# Set whether 1: to use PROJ projection library (set parameters in GBR.cc) or 0: to use old projection
usePROJ = 1
# If using PROJ, set projection systems:
lonlat_projection = "+proj=latlong +ellps=WGS84"
mesh_projection = "+proj=utm +ellps=WGS84 +zone=55 +south"
##

print(colored('Loading simulation parameters ...',"red"))
simInfoFile = open(simOutputDir+'/simulationInfo.txt','r')
simInfo = simInfoFile.readlines()
simDt = float(simInfo[0])
simTotIter = int(simInfo[1])
simExportGap = int(simInfo[2])
simStartTime = int(float(simInfo[3]))
#simOutputDir = simInfo[4].strip()
simMesh = simInfo[5].strip()
simInfoFile.close()

# Lauch Makefile to generate meshes and download forcings
#if(Msg.GetCommRank() == 0):
  #try : os.system("make "+reefMapFilename);
  #except: 0;
  #try : os.mkdir(LP_OutputDir);
  #except: 0;

if(not os.path.exists('./LPT_Functions')):
  try : os.mkdir('./LPT_Functions')
  except: 0;
  
simSolutions = [simExportGap*solutionNumber for solutionNumber in range(0,simTotIter/simExportGap+1)]
#-------------------------------------------------------------------------------

#--------- SET LPT TIME PARAMETERS ---------------------------------------------
print(colored('Loading particle-tracking parameters ...',"red"))
LP_dt = 90.                                            #(double) in seconds
LP_TotTime = 4*7*24.0                                   #(double) in hours
LP_Export = int( (12.0*3600.0)/LP_dt )                        #(int) export particle positions every ... steps
LP_Export_ConnectivityMatrix = int( (24.0*3600.0)/LP_dt )     #(int) export connectivity matrix every ... steps
LP_TotIter = int(LP_TotTime*3600.0/LP_dt)
LP_StartTime = simStartTime + (7.0*24.0*3600.0)              #in seconds
LP_EndTime = LP_StartTime + LP_TotTime*3600.0
simItersPerLPIter = LP_dt/simDt
#-------------------------------------------------------------------------------

#--------- SEEDING, MORTALITY, SETTLING AND SWIMMING PARAMETERS ----------------
# SELECT REEF TYPES TO SEED AND SETTLE OVER
#0: shallow reefs, 1: deep reefs, 2: both
seedOverWhichReefs = 0
settleOverWhichReefs = 0
# SET SEEDING CONCENTRATION (nb particles/kmsq)
concentration = 2.0
# SET MINIMUM NUMBER OF PARTICLES TO SEED ON EACH REEF AT EACH SEEDING EVENT
minNbParticlesPerReef = 200
  # Do we seed all in one go (0) or over a time interval (1)
seedOverTime = 0
  # If seedOverTime==1, specify the total time interval (seconds) over which to seed, and how many times to seed within this interval
seedOverInterval = 24.0*3600.0
seedNbTimes = 4
# WHICH REEFS SHALL WE SEED?
#  -> if reefsToSeed=0, then all reefs will be seeded, else, can set which reefs to seed
reefsToSeed = fullMatrixDouble(1,1)
reefsToSeed.set(0,0,0)
#reefsToSeed.set(1,1,216) etc.

# MORTALITY PARAMETERS
mortParams = fullMatrixDouble(4,4)
#[0]: 0: mortality off, 1: const. mortality, 2: Weibull mortality
mortParams.set(0,0, 1.0)
#[1]: (max if Weibull) mortality rate (PER DAY)
  #G. retiformis: 0.09, A. gemmifera: 0.067
  #    Reef fish: 0.18 (James et al. 2002)
mortParams.set(0,1, 0.067)
#[2]: Weibull only (else unused): lambda parameter <-- (Connolly & Baird 2010)lambda = 1/b(Pinder et al. 1978)
#A. millepora: 0.043, A. valida: 0.019, P. daedalea: 0.060
mortParams.set(0,2, 0.043)
#[3]: Weibull only (else unused): nu parameter <-- (Connolly & Baird 2010)nu = c(Pinder et al. 1978)
#A. millepora: 0.57, A. valida: 0.46, P. daedalea: 0.72
mortParams.set(0,3, 0.57)

# SETTLING PARAMETERS
#settleParams = [0.0,0.0,0.0,0.0,0.0]
settleParams = fullMatrixDouble(8,8)
#array settleParams: 
#[0]: 0: settling off, 1: const. settling, 2: triangle settling, 3: progressive competence acquisition/loss
settleParams.set(0,0, 3.0)
#[1]: t_preComp (seconds) == for type-3: t_timeLag
# G. retiformis: 0.0, A. gemmifera: 3.47, A. millepora: 3.239, A. valida: 0.0, P. daedalea: 2.937
settleParams.set(0,1, 3.47*24.0*3600.0)
#[2]: t_postComp (for ALL settling types)
settleParams.set(0,2, 100.0*24.0*3600.0)
#[3]: t_settMax (t of triangle peak, -1 if [0] != 2)
settleParams.set(0,3, -1.0)
#[4]: (max if [0]=2) settling rate (PER dt) - **! must convert to .dt(-1) !**
settleParams.set(0,4, 1.0)
#[5]: (for type-3 only) rate of acquisition of competence (d-1)
# G. retiformis: 0.58, A. gemmifera: 0.39, A. millepora: 0.18, A. valida: 0.22, P. daedalea: 0.39
settleParams.set(0,5, 0.39)
#[6]: (for type-3 only) rate of loss of competence (d-1)
# G. retiformis: 0.096, A. gemmifera: 0.145, A. millepora: 0.05, A. valida: 0.031, P. daedalea: 0.099
settleParams.set(0,6, 0.145)
#[7]: settle over shallow, deep or both reef types
settleParams.set(0,7, settleOverWhichReefs)

# SWIMMING PARAMETERS
swimParams = fullMatrixDouble(4,4)
#[0]: 0: swimming off, 1: swimming on
swimParams.set(0,0, 0.0)
#[1]: swimming speed (m/s)
swimParams.set(0,1, 0.1)
#[2]: swim when within .. kms of a reef
swimParams.set(0,2, 2.0)
#[3]: swimming starts after .. hours
swimParams.set(0,3, 0.0)
#-------------------------------------------------------------------------------

#--------- LOAD MODEL & BATHYMETRY ---------------------------------------------
print(colored('Loading model ...',"red"))
model = GModel()
if usePROJ == 0:
  print " (Using old projection)"
  model.load("./Meshes/"+simMesh+"_tan.msh")
else:
  print " (Using PROJ for projections)"
  model.load("./Meshes/"+simMesh+"_utm.msh")
groups = dgGroupCollection(model, 2, 1)   #model, dimension, order
XYZ = function.getCoordinates()
XYZ_DC = dgDofContainer(groups, 3)
XYZ_DC.interpolate(XYZ)
XYZ_DC.exportMsh("XYZ_DC",0,0)

# Print size of smallest element:
maxEdge = 100000.0
maxMaxEdge = 0.0
minSurface = 1000000.0
maxSurface = 0.0
nbGroups = groups.getNbElementGroups()
for group in range( 0, nbGroups):
  groupOfElements = groups.getElementGroup(group)
  nbElements = groupOfElements.getNbElements()
  for element in range( 0, nbElements):
    el = groupOfElements.getElement(element)
    if el.getVolume() < minSurface:
      minSurface = el.getVolume()
    if el.getVolume() > maxSurface:
      maxSurface = el.getVolume()
    if el.maxEdge() < maxEdge:
      maxEdge = el.maxEdge()
    if el.maxEdge() > maxMaxEdge:
      maxMaxEdge = el.maxEdge()
print 'Smallest element\'s maxEdge is', maxEdge, 'm, Largest element\'s maxEdge is', maxMaxEdge, 'm'
print 'Smallest element\'s surface is', minSurface, 'msq, Largest element\'s surface is', maxSurface, 'msq'

bathDC = dgDofContainer(groups, 1);
if(os.path.exists("./Bath/"+simMesh+"_bath_smooth/"+simMesh+"_bath_smooth.idx")):
  print(colored('Importing bathymetry ...',"red"))
  bathDC.importIdx("./Bath/"+simMesh+"_bath_smooth/"+simMesh+"_bath_smooth.idx")
else:
  print(colored('**** ERROR: Bathymetry file does not exist! ****',"red"))
#-------------------------------------------------------------------------------

#------- SET UP DIFFUSIVITY ----------------------------------------------------
print(colored('Defining diffusivity ...',"red"))

CCode = """
#include "fullMatrix.h"
#include "function.h"
#include "dgGroupOfElements.h"
#include "math.h"
#include "MElement.h"

extern "C" {
  void diff(dataCacheMap *m, fullMatrix<double> &diff, fullMatrix<double> &alpha) {
    double delta, K;
    int ElementId;
    const dgGroupOfElements *group = m->getGroupOfElements();
    ElementId = m->getElementId();

    MElement *e = group->getElement(ElementId);
    delta = e->maxEdge();
    K = alpha(0,0) * 0.00020551 * pow(delta,1.15);
    for (int i=0; i<diff.size1(); i++) {
      diff.set (i, 0, K);
    }
  }
}

"""
tmpCLib = "./tempCLib.dylib"
if (Msg.GetCommRank() == 0 ) :
  functionC.buildLibrary (CCode, tmpCLib)
  
alpha = functionConstant(1.0)
diffFC = functionC(tmpCLib,"diff",1,[alpha])
diffDC = dgDofContainer(groups,1)
diffDC.L2Projection(diffFC)
diffDC.exportMsh("./LPT_Functions/DIFFUSIVITY_RAW_"+simMesh,0,0)

# Smooth diffusivity:
sys = linearSystemPETScDouble()
dofMan = dgDofManager.newCG (groups, 1, sys)
diffProjector = L2ProjectionContinuous(dofMan)
diffProjector.apply(diffDC,diffFC)
diffDC.exportMsh("./LPT_Functions/DIFFUSIVITY_smooooth_"+simMesh,0,0)
#-------------------------------------------------------------------------------

#------- LOAD REEF MAPS --------------------------------------------------------
reefsShallowDC = dgDofContainer(groups,1)
reefsDeepDC = dgDofContainer(groups,1)
reefsShallowDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapShallowFilenameTrunk+"_"+simMesh
reefsDeepDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapDeepFilenameTrunk+"_"+simMesh

#1. Load reef maps from data if they don't already exist
UTMProjector = None
if usePROJ == 1:
  UTMProjector = slimProjector(mesh_projection, lonlat_projection)
if(not os.path.exists(reefsShallowDCfilename+"_COMP_0.msh")):
  print(colored('Creating shallow reef map msh file ...',"red"))
  # 1. Create a reefs object (initialises object and creates reefs function)
  if usePROJ == 0:
    reefsObj = reefsRaw_oldProj("./Data/"+reefMapShallowFilename, EarthRadius)
  else:
    reefsObj = reefsRaw("./Data/"+reefMapShallowFilename, UTMProjector)
  # 2. Interpolate the function onto DC
  reefsShallowDC.interpolate(reefsObj)
  reefsShallowDC.exportMsh(reefsShallowDCfilename,0,0)
  # 3. Delete reefs object data
  reefsObj.clearReefsData()
else:
  print(colored('Loading shallow reef map from file ...',"red"))
  reefsShallowDC.importMsh(reefsShallowDCfilename)

if(not os.path.exists(reefsDeepDCfilename+"_COMP_0.msh")):
  print(colored('Creating deep reef map msh file ...',"red"))
  # 1. Create a reefs object (initialises object and creates reefs function)
  if usePROJ == 0:
    reefsObj = reefsRaw_oldProj("./Data/"+reefMapDeepFilename, EarthRadius)
  else:
    reefsObj = reefsRaw("./Data/"+reefMapDeepFilename, UTMProjector)
  # 2. Interpolate the function onto DC
  reefsDeepDC.interpolate(reefsObj)
  reefsDeepDC.exportMsh(reefsDeepDCfilename,0,0)
  # 3. Delete reefs object data
  reefsObj.clearReefsData()
else:
  print(colored('Loading deep reef map from file ...',"red"))
  reefsDeepDC.importMsh(reefsDeepDCfilename)
  
#2. Load reefs data
minReefNbFM = fullMatrixDouble()
maxReefNbFM = fullMatrixDouble()
reefsShallowDC.minmax(minReefNbFM,maxReefNbFM)
nbOfShallowReefs = int(maxReefNbFM(0,0))
print " --> HIGHEST SHALLOW REEF ID: ",nbOfShallowReefs, "<-- "
reefsDeepDC.minmax(minReefNbFM,maxReefNbFM)
nbOfDeepReefs = int(maxReefNbFM(0,0))
nbOfReefs = nbOfShallowReefs + nbOfDeepReefs
print "Maximum possible number of reefs: ",nbOfReefs, "(shallow:", nbOfShallowReefs, ", deep:", nbOfDeepReefs, ")"
#-------------------------------------------------------------------------------

#------- WRITE SIMULATION INFO TO FILE -----------------------------------------
LP_SimInfo = [str(nbOfDeepReefs),'\n'+str(LP_dt),'\n',str((LP_StartTime-simStartTime)/3600),'\n', str(LP_TotTime),'\n',str(LP_Export),'\n', str(simItersPerLPIter),'\n',str(mortParams(0,0)),'\n' ,str(mortParams(0,1)/(24.0*3600.0)),'\n',str(settleParams(0,0)),'\n','Pre-comptetency period:'+str(settleParams(0,1)/(24.0*3600.0))+' (days)\n', 'Post-competent at:'+str(settleParams(0,2)/(24.0*3600.0))+' days)\n', '(max) settling rate:'+str(settleParams(0,4))+' per dt\nConcentration:'+str(concentration)+'/kmsq\n', 'Seeding over: '+str(seedOverWhichReefs)+'\n', 'Settling over: '+str(settleOverWhichReefs)+'\n', 'Highest shallow reef Id: '+str(nbOfShallowReefs)+'\n']
LP_SimInfoFile = open(LP_OutputDir+'/LP_runInfo.txt','w')
LP_SimInfoFile.writelines(LP_SimInfo)
LP_SimInfoFile.close()
#-------------------------------------------------------------------------------

#------- LOAD MAPS OF DISTANCE AND DIRECTION TO NEAREST REEF -------------------
distToReefsDC = None
directionToReefsDC = None
reefsDualDC = None
nbOfShallowReefsFConst = None
if swimParams(0,0) > 0.1:
  print(colored('Calculating distances to nearest shallow reefs ...', "red"))
  distToReefsDC = dgDofContainer(groups, 1)
  directionToReefsDC = dgDofContainer(groups, 3)

  if(os.path.exists("./LPT_Functions/DISTANCE_TO_REEFS_"+simMesh+"_COMP_0.msh") and os.path.exists("./LPT_Functions/DIRECTION_TO_REEFS_"+simMesh+"_COMP_0.msh")):
    print'./LPT_Functions/DISTANCE_TO_REEFS_'+simMesh+'_COMP_0.msh and ./LPT_Functions/DIRECTION_TO_REEFS_'+simMesh+'_COMP_0.msh files already exist: importing them.'
    distToReefsDC.importMsh("./LPT_Functions/DISTANCE_TO_REEFS_"+simMesh)
    directionToReefsDC.importMsh("./LPT_Functions/DIRECTION_TO_REEFS_"+simMesh)
  else:
    #TODO: Fix distance_to_reef generation in stereo
    print'./LPT_Functions/DISTANCE_TO_REEFS_'+simMesh+'_COMP_0.msh and/or ./LPT_Functions/DIRECTION_TO_REEFS_'+simMesh+'_COMP_0.msh files do not exist. Generating dof container with distance and directions to nearest reefs ...'
    reefDistanceObj = reefDistance(groups, reefsShallowDC, XYZ_DC, EarthRadius)
    reefDistanceObj.getNearestReefMap(distToReefsDC, directionToReefsDC)
    distToReefsDC.exportMsh("./LPT_Functions/DISTANCE_TO_REEFS_"+simMesh,0,0)
    directionToReefsDC.exportMsh("./LPT_Functions/DIRECTION_TO_REEFS_"+simMesh,0,0)
    #Also save DIR-TO-REEFS as function, for visualisation
    directionToReefsDC.exportFunctionMsh(directionToReefsDC.getFunction(), "./LPT_Functions/DIRECTION_TO_REEFS_vector_"+simMesh,0,0)

nbOfShallowReefsFConst = functionConstant(nbOfShallowReefs)
reefsDualDC = dgDofContainer(groups, 1)
reefsDualFC = functionC(glib, "overwrite2", 1, [bathDC.getFunction(), reefsShallowDC.getFunction(), reefsDeepDC.getFunction(), nbOfShallowReefsFConst])
reefsDualDC.interpolate(reefsDualFC)
reefsDualDC.exportMsh("LPT_Functions/REEFS_RAW_LPT_DUAL_"+simMesh,0,0)
#-------------------------------------------------------------------------------

#--------- SEED PARTICLES ------------------------------------------------------
print(colored('Seeding particles ...',"red"))

secondsBetweenSeeds = 0
itersBetweenSeeds = 0
if seedOverTime == 1:
  secondsBetweenSeeds = seedOverInterval/float(seedNbTimes)
  itersBetweenSeeds = secondsBetweenSeeds/LP_dt
  concentration = concentration/float(seedNbTimes)
  print'  We are seeding particles PROGRESSIVELY OVER TIME:'
  print'  -> Seeding',seedNbTimes,'times, at intervals of', secondsBetweenSeeds/3600.0, 'hours ('+str(itersBetweenSeeds)+' iters)'
  print'  -> Total concentration:', concentration*seedNbTimes, ', concentration per seeding event:', concentration

# ***SEED OVER ALL REEFS***
initCM = fullMatrixDouble( nbOfReefs, nbOfReefs )
particles = particleArray()
particles.addParticlesWithReefMap(groups, reefsDualDC, initCM, reefsToSeed, concentration, minNbParticlesPerReef, seedOverWhichReefs, nbOfShallowReefs)
#particles.addParticlesWithDofContainer(groups, reefsDualDC, concentration, 200)
particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
dummyPrinterICM = connectivityMatrixPrint()
if settleParams(0,0) == 0:
  dummyPrinterICM.CMPrinter(initCM, LP_OutputDir+'/connectivityMatrix_initial' )
else:
  dummyPrinterICM.CMPrinter(initCM, LP_OutputDir+'/settlementMatrix_initial' )
# ******

# ***SEED AT A POINT*** (NB: Look at method used to print results if changing seeding method)
#UTMProjectorLonLatToUTM = slimProjector(lonlat_projection, mesh_projection)
#reefPointsLatLon = [ [-14.7,145.4,0], [-14.7,145.4,1] ]
#reefPointsCart = [ [0]*4 for i in range(len(reefPointsLatLon)) ]
#pj_latlong = pj_init_plus(lonlat_projection)
#pj_merc = pj_init_plus(mesh_projection)
#for i in range(0,len(reefPointsLatLon)):
  #reefPointsCart[i][3] = reefPointsLatLon[i][2]
  #src_points = fullMatrixDouble(2,2)
  #dest_points = fullMatrixDouble(2,2)
  #src_points.set(0,0,reefPointsLatLon[i][1]/180.0*pi)
  #src_points.set(1,1,reefPointsLatLon[i][0]/180.0*pi)
  #UTMProjectorLonLatToUTM.projectPoint(src_points, dest_points)
  #reefPointsCart[i][0] = dest_points(0,0)
  #reefPointsCart[i][1] = dest_points(1,1)

#particles = particleArray()
#for location in range(0,len(reefPointsCart)):
  #particles.addParticlesAtPoint(groups, 100, reefPointsCart[location][0], reefPointsCart[location][1], 0, reefPointsCart[location][3])
#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 reefs in reef map:', nbOfReefs, '(highest reef Id on reefs DC)'
print 'Number of particles:', particles.printNbParticles()
print '-----------------------'
print ''
print '-----------------------'
print 'BIOLOGICAL PARAMETERS:'
print 'Mortality flag:', mortParams(0,0), '  (0: OFF, 1: CONST, 2: WEIBULL)'
if (int(mortParams(0,0)) == 1):
  print 'Mortality rate (const):', mortParams(0,1), '(per day)'
elif (int(mortParams(0,0)) == 2):
  print 'Weibull lambda:', mortParams(0,2), ', Weibull nu:', mortParams(0,3)
print '-----------------------'
print 'Settling flag: ', settleParams(0,0), '  (0: OFF, 1: CONST, 2: TRIANGLE, 3: PROGRESSIVE)'
if (settleParams(0,0) > 0):
  print 'Settlement competence begins: ', (settleParams(0,1)/(24.0*3600.0)), '(days) & ends: ', settleParams(0,2)/(24.0*3600.0), '(days)'
  print 'Settlement rate (max):', settleParams(0,4), '(per dt)'
if (settleParams(0,0) == 3):
  print 'Rate of acquisition of competence:', settleParams(0,5), '(per day) | Rate of loss of competence:', settleParams(0,6), '(per day)'
print 'Seeding over:', seedOverWhichReefs, '  (0: SHALLOW REEFS ONLY, 1: DEEP REEFS ONLY, 2: BOTH)'
print 'Settling over:', settleOverWhichReefs, '  (0: SHALLOW REEFS ONLY, 1: DEEP REEFS ONLY, 2: BOTH)'
print '-----------------------'
print 'Swimming flag:', swimParams(0,0), '  (0: OFF, 1: ON)'
if (swimParams(0,0) > 0):
  print 'Swimming speed: ', swimParams(0,1), 'm/s'
  print 'Swimming max distance: ',swimParams(0,2), 'km'
  print 'Swimming starts after: ', swimParams(0,3), 'hrs'
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)
simSolutionUnder = dgDofContainer(groups, 3)
simSolutionOver = dgDofContainer(groups, 3)
#Keep track of the solutions loaded into the three DCs:
loadedSolutions = [-1] * 3 #[0]:simSolution, [1]:simSolutionOver, [2]:simSolutionUnder

connectivityMatrix = fullMatrixDouble( nbOfReefs, nbOfReefs )
dummyPrinter = connectivityMatrixPrint()
#-------------------------------------------------------------------------------

#------ 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)
particleTracker.addMortality(mortParams)
particleTracker.addSettling(settleParams, reefsDualDC, nbOfShallowReefs)
#particleTracker.addSwimming(swimParams, reefsDualDC, distToReefsDC, directionToReefsDC)
#-------------------------------------------------------------------------------

#-------- PARTICLE LOOP --------------------------------------------------------
print 'Beginning loop ...'
nextIterToSeed=itersBetweenSeeds
seedNbTimes = seedNbTimes-1 #seedNbTimes now counts down how many more times we need to seed
startcpu = time.clock()

for n in range(1,LP_TotIter+1):
  #Check whether to seed more particles
  if seedOverTime == 1 and seedNbTimes > 0 and n >= nextIterToSeed:
    print '*** Seeding more particles ... ***'
    particles.addParticlesWithReefMap(groups, reefsDualDC, initCM, reefsToSeed, concentration, seedOverWhichReefs, nbOfShallowReefs)
    nextIterToSeed = nextIterToSeed + itersBetweenSeeds
    seedNbTimes = seedNbTimes - 1
    print '***', seedNbTimes,'seeding events left ***'
  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) | Loaded solutions:',loadedSolutions[0],loadedSolutions[1],loadedSolutions[2],'| CPU time:', time.clock()-startcpu
  
  if iterNumberWanted in simSolutions:
    
    if iterNumberWanted in loadedSolutions:
      if iterNumberWanted == loadedSolutions[0]:
        print '->  Exact solution for simSolution (iteration %06d'%( iterNumberWanted )+') already loaded in simSolution. Using this.'
      elif iterNumberWanted == loadedSolutions[1]:
        print '->  Exact solution for simSolution (iteration %06d'%( iterNumberWanted )+') already loaded in simSolutionOver. Using this.'
        simSolution.copy(simSolutionOver)
        loadedSolutions[0] = iterNumberWanted
      elif iterNumberWanted == loadedSolutions[2]:
        print '->  Exact solution for simSolution (iteration %06d'%( iterNumberWanted )+') already loaded in simSolutionUnder. Are we going back in time??? Using this anyway...'
        simSolution.copy(simSolutionUnder)
        loadedSolutions[0] = iterNumberWanted
      else:
        print 'Problem in loadedSolutions!'
    else:
      print '->  simSolution not already loaded. Loading it from file. m is:', m, '. File is (gbr-%06d'%( iterNumberWanted )+').'
      simSolution.importIdx(simOutputDir+'/mr-gbr/mr-gbr-%06d'%( iterNumberWanted )+'.idx')
      loadedSolutions[0] = iterNumberWanted
  
  else:
    
    l=0
    while (not iterNumberWanted-l in simSolutions):
      l=l+1
    
    if iterNumberWanted-l in loadedSolutions:
      if iterNumberWanted-l == loadedSolutions[0]:
        print '->  Exact solution for simSolutionUnder (iteration %06d'%( iterNumberWanted-l )+') already loaded in simSolution. Using this.'
        simSolutionUnder.copy(simSolution)
        loadedSolutions[2] = iterNumberWanted-l
        t_solUnder = simStartTime + float(iterNumberWanted-l)*simDt
      elif iterNumberWanted-l == loadedSolutions[1]:
        print '->  Exact solution for simSolutionUnder (iteration %06d'%( iterNumberWanted-l )+') already loaded in simSolutionOver. Using this.'
        simSolutionUnder.copy(simSolutionOver)
        loadedSolutions[2] = iterNumberWanted-l
        t_solUnder = simStartTime + float(iterNumberWanted-l)*simDt
      elif iterNumberWanted-l == loadedSolutions[2]:
        print '->  Exact solution for simSolutionUnder (iteration %06d'%( iterNumberWanted-l )+') already loaded in simSolutionUnder. Using this.'
      else:
        print 'Problem in loadedSolutions!'
    else:
      print '->  simSolutionUnder not already loaded. Loading it from file. l is:', l, '. File is (gbr-%06d'%( iterNumberWanted-l )+').'
      t_solUnder = simStartTime + float(iterNumberWanted-l)*simDt
      simSolutionUnder.importIdx(simOutputDir+'/mr-gbr/mr-gbr-%06d'%( iterNumberWanted-l )+'.idx')
      loadedSolutions[2] = iterNumberWanted-l
    
    m=0
    while (not iterNumberWanted+m in simSolutions):
      m=m+1
    
    if iterNumberWanted+m in loadedSolutions:
      if iterNumberWanted+m == loadedSolutions[0]:
        print '->  Exact solution for simSolutionOver (iteration %06d'%( iterNumberWanted+m )+') already loaded in simSolution. Using this.'
        simSolutionOver.copy(simSolution)
        loadedSolutions[1] = iterNumberWanted+m
        t_solOver = simStartTime + float(iterNumberWanted+m)*simDt
      elif iterNumberWanted+m == loadedSolutions[1]:
        print '->  Exact solution for simSolutionOver (iteration %06d'%( iterNumberWanted+m )+') already loaded in simSolutionOver. Using this.'
      elif iterNumberWanted+m == loadedSolutions[2]:
        print '->  Exact solution for simSolutionOver (iteration %06d'%( iterNumberWanted+m )+') already loaded in simSolutionUnder. Are we going back in time??? Using this anyway...'
        simSolutionOver.copy(simSolutionUnder)
        loadedSolutions[1] = iterNumberWanted+m
        t_solOver = simStartTime + float(iterNumberWanted+m)*simDt
      else:
        print 'Problem in loadedSolutions!'
    else:
      print '->  simSolutionOver not already loaded. Loading it from file. m is:', m, '. File is (gbr-%06d'%( iterNumberWanted+m )+').'
      t_solOver = simStartTime + float(iterNumberWanted+m)*simDt
      simSolutionOver.importIdx(simOutputDir+'/mr-gbr/mr-gbr-%06d'%( iterNumberWanted+m )+'.idx')
      loadedSolutions[1] = iterNumberWanted+m
    
    scaleFactor = 1.0 - (t-t_solUnder)/(t_solOver-t_solUnder)
    print ' Info: simSolutionUnder scale factor is', scaleFactor
    simSolution.copy(simSolutionUnder)                   #Put simSolUnder in simSolution
    simSolution.scale(scaleFactor)                       #Scale simSolUnder (dof * relative dist. in t of sol. to simSolUnder)
    simSolution.axpy(simSolutionOver,1.0-scaleFactor)    #Scale simSolOver and add to scaled simSolUnder
    loadedSolutions[0] = -1
    
    #simSolution.exportMsh('./_LP-%06d'%(n), t, n)
  
  particleTracker.particleMove(simSolution, connectivityMatrix, n)
  # Print output:
  if ( n%(LP_Export)  == 0 ):
    particles.printPositions(LP_OutputDir+'/particlesAlive_%06d'%( n ), 'pos', 0)
    if (int(settleParams(0,0)) != 0):
      particles.printPositions(LP_OutputDir+'/particlesSettled_%06d'%( n ), 'pos', 1)
    if (int(mortParams(0,0)) != 0):
      particles.printPositions(LP_OutputDir+'/particlesDead_%06d'%( n ), 'pos', 2)
  # *** More useful only for SEEDING AT SINGLE POINT ***
#  if ( n%LP_Export  == 0 ):
#    particles.printPositions(LP_OutputDir+'/particleOutput_%06d'%( n ), 'dat', 0)
#    particles.printPositionsFollowSingle(LP_OutputDir+'/particleOutputFollow', 'pos', n, LP_TotIter, LP_dt, 0)
#    particles.printPositionsFollowSingle(LP_OutputDir+'/particleOutputFollow', 'pos', n, LP_TotIter, LP_dt, 400)
#    particles.printPositionsFollowSingle(LP_OutputDir+'/particleOutputFollow', 'pos', n, LP_TotIter, LP_dt, 800)
  if ( n%LP_Export_ConnectivityMatrix  == 0 ):
    if settleParams(0,0) == 0:
      dummyPrinter.CMPrinter(connectivityMatrix, LP_OutputDir+'/connectivityMatrix%06d'%( n ) )
    else:
      dummyPrinter.CMPrinter(connectivityMatrix, LP_OutputDir+'/settlementMatrix%06d'%( n ) )
#-------------------------------------------------------------------------------

print(colored('Particle tracking finished. Closing module.',"yellow"))
Msg.Exit(0)