LPT_Caller_DualReefMap.py 28.8 KB
Newer Older
thomas's avatar
thomas committed
1
2
3
# 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
thomas's avatar
thomas committed
4
# Mesh must be in planar coords in Meshes/ folder. Stereo-coords: WIP
thomas's avatar
thomas committed
5
6
7
8
# 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
#
thomas's avatar
thomas committed
9
# !!! NB: Must update t_preComp in settleParams, EVEN IF settleType=0 !!!
thomas's avatar
thomas committed
10
11
12
13
14
15
16
17
18
19
20
21
22
#      --> 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 ''

thomas's avatar
thomas committed
23
#--------- SETUP LIBRARIES  ----------------------------------------------------
thomas's avatar
thomas committed
24
25
if (Msg.GetCommRank()==0):
  print ''
thomas's avatar
thomas committed
26
  print(colored('Compiling Libraries ...', "red"))
thomas's avatar
thomas committed
27
28
29
30
31
  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"
thomas's avatar
thomas committed
32
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
33

thomas's avatar
thomas committed
34
35
36
37
38
39
40
41
42
#---------- 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
thomas's avatar
thomas committed
43
44
45
cpuExt=""
if ( len(sys.argv)>1 ):
  if sys.argv[1] == "-c":
46
47
    cpuExt = "/CPU"+sys.argv[2]
    time.sleep(float(sys.argv[2])*2) #Avoid simultaneous simulations crashing each other
thomas's avatar
thomas committed
48

thomas's avatar
thomas committed
49
50
51
52
53
54
LP_OutputDir = LP_OutputDir + cpuExt
if(not os.path.exists(LP_OutputDir)):
  try : os.mkdir(LP_OutputDir)
  except: 0;
##
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
55
56


57
58
#reefMapShallowFilenameTrunk = "SGBR_ReefMap_2121_Solid9200.dat"
#reefMapDeepFilenameTrunk = "SGBR_DeepReefMap_2121_Solid9200.dat"
thomas's avatar
thomas committed
59
60
reefMapShallowFilenameTrunk = "GBR_ReefMap_2966_Solid9200"
reefMapDeepFilenameTrunk = "GBR_DeepReefMap_2966_Solid9200"
thomas's avatar
thomas committed
61
62
reefMapShallowFilename = reefMapShallowFilenameTrunk+".dat"
reefMapDeepFilename = reefMapDeepFilenameTrunk+".dat"
thomas's avatar
thomas committed
63

thomas's avatar
thomas committed
64
65
# 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
thomas's avatar
thomas committed
66
# Set whether 1: to use PROJ projection library (set parameters in GBR.cc) or 0: to use old projection
67
68
69
usePROJ = 1
# If using PROJ, set projection systems:
lonlat_projection = "+proj=latlong +ellps=WGS84"
thomas's avatar
thomas committed
70
mesh_projection = "+proj=utm +ellps=WGS84 +zone=55 +south"
thomas's avatar
thomas committed
71
72
73
74
75
76
77
78
79
80
##

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(simInfo[3])
#simOutputDir = simInfo[4].strip()
thomas's avatar
thomas committed
81
simMesh = simInfo[5].strip()
thomas's avatar
thomas committed
82
83
84
85
86
87
88
89
90
91
92
93
94
95
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)]
thomas's avatar
thomas committed
96
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
97

thomas's avatar
thomas committed
98
#--------- SET LPT TIME PARAMETERS ---------------------------------------------
thomas's avatar
thomas committed
99
print(colored('Loading particle-tracking parameters ...',"red"))
thomas's avatar
thomas committed
100
101
LP_dt = 90.                                            #(double) in seconds
LP_TotTime = 4*7*24.0                                   #(double) in hours
thomas's avatar
thomas committed
102
LP_Export = int( (12.0*3600.0)/LP_dt )                        #(int) export particle positions every ... steps
thomas's avatar
thomas committed
103
104
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)
thomas's avatar
thomas committed
105
LP_StartTime = simStartTime + (7.0*24.0*3600.0)              #in seconds
thomas's avatar
thomas committed
106
107
LP_EndTime = LP_StartTime + LP_TotTime*3600.0
simItersPerLPIter = LP_dt/simDt
thomas's avatar
thomas committed
108
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
109

thomas's avatar
thomas committed
110
111
112
113
114
115
116
117
118
119
120
121
122
#--------- SEEDING, MORTALITY, SETTLING AND SWIMMING PARAMETERS ----------------
# SELECT REEF TYPES TO SEED AND SETTLE OVER
#0: shallow reefs, 1: deep reefs, 2: both
seedOverWhichReefs = 2
settleOverWhichReefs = 2
# SET SEEDING CONCENTRATION (nb particles/kmsq)
concentration = 2.0
  # 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?
thomas's avatar
thomas committed
123
124
125
126
127
128
129
130
#  -> 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
thomas's avatar
thomas committed
131
mortParams.set(0,0, 1.0)
thomas's avatar
thomas committed
132
#[1]: (max if Weibull) mortality rate (PER DAY)
133
  #G. retiformis: 0.09, A. gemmifera: 0.067
thomas's avatar
thomas committed
134
  #    Reef fish: 0.18 (James et al. 2002)
thomas's avatar
thomas committed
135
mortParams.set(0,1, 0.067)
thomas's avatar
thomas committed
136
#[2]: Weibull only (else unused): lambda parameter <-- (Connolly & Baird 2010)lambda = 1/b(Pinder et al. 1978)
137
#A. millepora: 0.043, A. valida: 0.019, P. daedalea: 0.060
thomas's avatar
thomas committed
138
139
mortParams.set(0,2, 0.043)
#[3]: Weibull only (else unused): nu parameter <-- (Connolly & Baird 2010)nu = c(Pinder et al. 1978)
140
#A. millepora: 0.57, A. valida: 0.46, P. daedalea: 0.72
thomas's avatar
thomas committed
141
142
143
144
145
146
147
148
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)
thomas's avatar
thomas committed
149
#[1]: t_preComp (seconds) == for type-3: t_timeLag
150
# G. retiformis: 0.0, A. gemmifera: 3.47, A. millepora: 3.239, A. valida: 0.0, P. daedalea: 2.937
thomas's avatar
thomas committed
151
152
153
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)
thomas's avatar
thomas committed
154
155
156
157
158
#[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)
159
# G. retiformis: 0.58, A. gemmifera: 0.39, A. millepora: 0.18, A. valida: 0.22, P. daedalea: 0.39
thomas's avatar
thomas committed
160
settleParams.set(0,5, 0.39)
thomas's avatar
thomas committed
161
#[6]: (for type-3 only) rate of loss of competence (d-1)
162
# G. retiformis: 0.096, A. gemmifera: 0.145, A. millepora: 0.05, A. valida: 0.031, P. daedalea: 0.099
thomas's avatar
thomas committed
163
settleParams.set(0,6, 0.145)
thomas's avatar
thomas committed
164
165
166
167
168
169
170
171
172
173
174
175
176
#[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)
thomas's avatar
thomas committed
177
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
178

thomas's avatar
thomas committed
179
#--------- LOAD MODEL & BATHYMETRY ---------------------------------------------
thomas's avatar
thomas committed
180
181
print(colored('Loading model ...',"red"))
model = GModel()
182
183
if usePROJ == 0:
  print " (Using old projection)"
thomas's avatar
thomas committed
184
185
  model.load("./Meshes/"+simMesh+"_tan.msh")
else:
186
187
  print " (Using PROJ for projections)"
  model.load("./Meshes/"+simMesh+"_utm.msh")
thomas's avatar
thomas committed
188
189
190
191
192
193
194
195
196
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)

bathDC = dgDofContainer(groups, 1);

if(os.path.exists("./Bath/"+simMesh+"_bath_smooth/"+simMesh+"_bath_smooth.idx")):
thomas's avatar
thomas committed
197
  print(colored('Importing bathymetry ...',"red"))
thomas's avatar
thomas committed
198
199
200
  bathDC.readMsh("./Bath/"+simMesh+"_bath_smooth/"+simMesh+"_bath_smooth.idx")
else:
  print(colored('**** ERROR: Bathymetry file does not exist! ****',"red"))
thomas's avatar
thomas committed
201
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
202

thomas's avatar
thomas committed
203
#------- SET UP DIFFUSIVITY ----------------------------------------------------
thomas's avatar
thomas committed
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
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);
    }
  }
}

"""
230
tmpCLib = "./tempCLib.dylib"
thomas's avatar
thomas committed
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
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)
thomas's avatar
thomas committed
246
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
247

thomas's avatar
thomas committed
248
#------- LOAD REEF MAPS --------------------------------------------------------
thomas's avatar
thomas committed
249
250
reefsShallowDC = dgDofContainer(groups,1)
reefsDeepDC = dgDofContainer(groups,1)
thomas's avatar
thomas committed
251
252
reefsShallowDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapShallowFilenameTrunk+"_"+simMesh
reefsDeepDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapDeepFilenameTrunk+"_"+simMesh
thomas's avatar
thomas committed
253

thomas's avatar
thomas committed
254
#1. Load reef maps from data if they don't already exist
255
256
257
UTMProjector = None
if usePROJ == 1:
  UTMProjector = slimProjector(mesh_projection, lonlat_projection)
thomas's avatar
thomas committed
258
259
260
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)
261
262
263
264
  if usePROJ == 0:
    reefsObj = reefsRaw_oldProj("./Data/"+reefMapShallowFilename, EarthRadius)
  else:
    reefsObj = reefsRaw("./Data/"+reefMapShallowFilename, UTMProjector)
thomas's avatar
thomas committed
265
266
267
268
269
270
271
272
273
274
275
276
  # 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)
277
278
279
280
  if usePROJ == 0:
    reefsObj = reefsRaw_oldProj("./Data/"+reefMapDeepFilename, EarthRadius)
  else:
    reefsObj = reefsRaw("./Data/"+reefMapDeepFilename, UTMProjector)
thomas's avatar
thomas committed
281
282
283
284
285
286
287
288
289
  # 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)
  
thomas's avatar
thomas committed
290
#2. Load reefs data
thomas's avatar
thomas committed
291
292
293
294
295
296
297
298
299
300
301
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, ")"

# Record simulation info to file
thomas's avatar
thomas committed
302
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']
thomas's avatar
thomas committed
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
LP_SimInfoFile = open(LP_OutputDir+'/LP_runInfo.txt','w')
LP_SimInfoFile.writelines(LP_SimInfo)
LP_SimInfoFile.close()

# 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'

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)

#dir2RF = directionToReefsDC.getFunction()
#dirFunction = functionC(glib,"current2",3,[directionToReefsDC.getFunction()])
#directionToReefsDC.exportMshNodeData("./DIRECTION_TO_REEFS",0,0, '')

nbOfShallowReefsFConst = functionConstant(nbOfShallowReefs)
reefsDualDC = dgDofContainer(groups, 1)
thomas's avatar
thomas committed
353
reefsDualFC = functionC(glib, "overwrite2", 1, [bathDC.getFunction(), reefsShallowDC.getFunction(), reefsDeepDC.getFunction(), nbOfShallowReefsFConst])
thomas's avatar
thomas committed
354
355
reefsDualDC.interpolate(reefsDualFC)
reefsDualDC.exportMsh("LPT_Functions/REEFS_RAW_LPT_DUAL_"+simMesh,0,0)
thomas's avatar
thomas committed
356
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
357

thomas's avatar
thomas committed
358
#--------- SEED PARTICLES ------------------------------------------------------
thomas's avatar
thomas committed
359
print(colored('Seeding particles ...',"red"))
thomas's avatar
thomas committed
360
361
362
363
364
365
366
367
368
369
370

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

thomas's avatar
thomas committed
371
# ***SEED OVER ALL REEFS***
thomas's avatar
thomas committed
372
#initCM = fullMatrixDouble( nbOfReefs, nbOfReefs )
thomas's avatar
thomas committed
373
374
#particles = particleArray()
#particles.addParticlesWithReefMap(groups, reefsDualDC, initCM, reefsToSeed, concentration, seedOverWhichReefs, nbOfShallowReefs)
thomas's avatar
thomas committed
375
376
377
378
379
380
#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' )
thomas's avatar
thomas committed
381
# ******
thomas's avatar
thomas committed
382
383

# ***SEED AT A POINT*** (NB: Look at method used to print results if changing seeding method)
thomas's avatar
thomas committed
384
UTMProjectorLonLatToUTM = slimProjector(lonlat_projection, mesh_projection)
thomas's avatar
thomas committed
385
reefPointsLatLon = [ [-14.7,145.4,0], [-14.7,145.4,1] ]
thomas's avatar
thomas committed
386
387
388
389
390
391
392
393
394
395
396
397
398
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)

thomas's avatar
thomas committed
399
particles = particleArray()
thomas's avatar
thomas committed
400
for location in range(0,len(reefPointsCart)):
thomas's avatar
thomas committed
401
  particles.addParticlesAtPoint(groups, 100, reefPointsCart[location][0], reefPointsCart[location][1], 0, reefPointsCart[location][3])
thomas's avatar
thomas committed
402
particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
thomas's avatar
thomas committed
403
404
405
particles.printPositions(LP_OutputDir+'/seededParticles','dat',0)
# ******
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
406

thomas's avatar
thomas committed
407
#-------- PRINT SIMULATION INFO ------------------------------------------------
thomas's avatar
thomas committed
408
409
print ''
print '-----------------------'
thomas's avatar
thomas committed
410
print 'HYDRODYNAMICS PARAMETERS:'
thomas's avatar
thomas committed
411
412
413
414
415
416
417
418
419
420
421
422
423
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 '-----------------------'
thomas's avatar
thomas committed
424
425
426
print ''
print '-----------------------'
print 'BIOLOGICAL PARAMETERS:'
thomas's avatar
thomas committed
427
428
print 'Mortality flag:', mortParams(0,0), '  (0: OFF, 1: CONST, 2: WEIBULL)'
if (int(mortParams(0,0)) == 1):
thomas's avatar
thomas committed
429
  print 'Mortality rate (const):', mortParams(0,1), '(per day)'
thomas's avatar
thomas committed
430
431
432
elif (int(mortParams(0,0)) == 2):
  print 'Weibull lambda:', mortParams(0,2), ', Weibull nu:', mortParams(0,3)
print '-----------------------'
thomas's avatar
thomas committed
433
print 'Settling flag: ', settleParams(0,0), '  (0: OFF, 1: CONST, 2: TRIANGLE, 3: PROGRESSIVE)'
thomas's avatar
thomas committed
434
435
436
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)'
thomas's avatar
thomas committed
437
438
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)'
thomas's avatar
thomas committed
439
440
441
442
443
444
445
446
447
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 '-----------------------'
thomas's avatar
thomas committed
448
print ''
thomas's avatar
thomas committed
449
450
451
452
453
454
455
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 ''
thomas's avatar
thomas committed
456
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
457

thomas's avatar
thomas committed
458
#------- SET UP HYDRO DOFCONTAINERS --------------------------------------------
thomas's avatar
thomas committed
459
460
461
simSolution = dgDofContainer(groups, 3)
simSolutionUnder = dgDofContainer(groups, 3)
simSolutionOver = dgDofContainer(groups, 3)
thomas's avatar
thomas committed
462
463
#Keep track of the solutions loaded into the three DCs:
loadedSolutions = [-1] * 3 #[0]:simSolution, [1]:simSolutionOver, [2]:simSolutionUnder
thomas's avatar
thomas committed
464
465
466

connectivityMatrix = fullMatrixDouble( nbOfReefs, nbOfReefs )
dummyPrinter = connectivityMatrixPrint()
thomas's avatar
thomas committed
467
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
468

thomas's avatar
thomas committed
469
#------ SET UP PARTICLE TRACKER OBJECT -----------------------------------------
thomas's avatar
thomas committed
470
471
#Create particle tracker object
print 'Initialising particle tracker ...'
thomas's avatar
thomas committed
472
473
474
475
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)
thomas's avatar
thomas committed
476
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
477

thomas's avatar
thomas committed
478
#-------- PARTICLE LOOP --------------------------------------------------------
thomas's avatar
thomas committed
479
print 'Beginning loop ...'
thomas's avatar
thomas committed
480
481
482
483
nextIterToSeed=itersBetweenSeeds
seedNbTimes = seedNbTimes-1 #seedNbTimes now counts down how many more times we need to seed
startcpu = time.clock()

thomas's avatar
thomas committed
484
for n in range(1,LP_TotIter+1):
thomas's avatar
thomas committed
485
486
487
488
489
490
491
  #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 ***'
thomas's avatar
thomas committed
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
  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.readMsh(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.readMsh(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.readMsh(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)
thomas's avatar
thomas committed
583
    if (int(settleParams(0,0)) != 0):
thomas's avatar
thomas committed
584
      particles.printPositions(LP_OutputDir+'/particlesSettled_%06d'%( n ), 'pos', 1)
thomas's avatar
thomas committed
585
    if (int(mortParams(0,0)) != 0):
thomas's avatar
thomas committed
586
587
588
589
590
591
592
593
594
595
596
597
      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 ) )
thomas's avatar
thomas committed
598
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
599
600
601

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