LPT_Caller_DualReefMap.py 24.3 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
23
24
#      --> 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 ''

if (Msg.GetCommRank()==0):
  print ''
thomas's avatar
thomas committed
25
  print(colored('Compiling Libraries ...', "red"))
thomas's avatar
thomas committed
26
27
28
29
30
31
32
33
34
35
36
37
38
  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"

cpuExt=""
if ( len(sys.argv)>1 ):
  if sys.argv[1] == "-c":
    cpuExt = "/"+sys.argv[2]

# Define simOutput directory & particle concentration

thomas's avatar
thomas committed
39
40
41
simOutputDir = '/scratch/cthomas/Hydro_Realistic/gbr500K_sma_1svIn_bDf10_allTides_windOn_28112007_5wks/'
#simOutputDir = 'Output10K_4wks_windsOn_allTides_16112000'

thomas's avatar
thomas committed
42
LP_OutputDir = '/scratch/cthomas/temp'#simOutputDir +'/LPT_Output/SettProg_Mort007'
thomas's avatar
thomas committed
43
#LP_OutputDir = 'temp'
thomas's avatar
thomas committed
44
45
46
47


##
concentration = 2.0  #No. of particles to seed per kmsq
thomas's avatar
thomas committed
48
49
50
51
52
#reefMapShallowFilename = "SGBR_ReefMap_2121_Solid9200.dat"
reefMapShallowFilenameTrunk = "CGBR_ReefMap_1093_Cropped_Sma_Solid9200"
reefMapDeepFilenameTrunk = "CGBR_DeepReefMap_1093_Cropped_Sma_Solid9200"
reefMapShallowFilename = reefMapShallowFilenameTrunk+".dat"
reefMapDeepFilename = reefMapDeepFilenameTrunk+".dat"
thomas's avatar
thomas committed
53

thomas's avatar
thomas committed
54
55
#Radius: 6.37101e6 for old meshes (sma_zoomed), else 6.371e6
EarthRadius=6.37101e6
thomas's avatar
thomas committed
56
stereoFlag = 0
thomas's avatar
thomas committed
57
58
59
60
61

##

# SELECT REEF TYPES TO SEED AND SETTLE OVER
#0: shallow reefs, 1: deep reefs, 2: both
thomas's avatar
thomas committed
62
63
seedOverWhichReefs = 1
settleOverWhichReefs = 2
thomas's avatar
thomas committed
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
##

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

XYZ = function.getCoordinates()
lonLatDegrees = functionC(glib,"lonLatDegrees",3,[XYZ])
stereoToLonLatDegrees = functionC(glib,"stereoToLonLatDegrees",3,[XYZ])

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
86
simMesh = simInfo[5].strip()
thomas's avatar
thomas committed
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
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)]

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

thomas's avatar
thomas committed
112
# REEF SEEDING PARAM
thomas's avatar
thomas committed
113
114
115
116
117
118
119
120
#  -> 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
121
mortParams.set(0,0, 1.0)
thomas's avatar
thomas committed
122
#[1]: (max if Weibull) mortality rate (PER DAY)
thomas's avatar
thomas committed
123
  #(Connolly 2k11): 0.09 d-1 for G. retiformis, 0.067 for A. Gemmifera
thomas's avatar
thomas committed
124
  #    Reef fish: 0.18 (James et al. 2002)
thomas's avatar
thomas committed
125
mortParams.set(0,1, 0.067)
thomas's avatar
thomas committed
126
127
128
129
130
131
132
133
134
135
136
137
138
#[2]: Weibull only (else unused): lambda parameter <-- (Connolly & Baird 2010)lambda = 1/b(Pinder et al. 1978)
#A. Millepora: 0.043
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
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
139
#[1]: t_preComp (seconds) == for type-3: t_timeLag
thomas's avatar
thomas committed
140
# G. Retiformis: 0.0, A. Gemmifera: 3.47
thomas's avatar
thomas committed
141
142
143
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
144
145
146
147
148
#[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)
thomas's avatar
thomas committed
149
# G. Retiformis: 0.58, A. Gemmifera: 0.39
thomas's avatar
thomas committed
150
settleParams.set(0,5, 0.39)
thomas's avatar
thomas committed
151
#[6]: (for type-3 only) rate of loss of competence (d-1)
thomas's avatar
thomas committed
152
# G. Retiformis: 0.096, A. Gemmifera: 0.145
thomas's avatar
thomas committed
153
settleParams.set(0,6, 0.145)
thomas's avatar
thomas committed
154
155
156
157
158
159
160
161
162
163
164
165
166
167
#[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
168
#STEREO: still need to convert swimming speeds to stereo
thomas's avatar
thomas committed
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187

print(colored('Loading model ...',"red"))
model = GModel()
if stereoFlag == 0:
  print " (Cartesian co-ordinates)"
  model.load("./Meshes/"+simMesh+"_tan.msh")
else:
  print " (Stereographic co-ordinates)"
  model.load("./Meshes/"+simMesh+".msh")
groups = dgGroupCollection(model, 2, 1)   #model, dimension, order
#groups.buildGroupsOfInterfaces()
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
188
  print(colored('Importing bathymetry ...',"red"))
thomas's avatar
thomas committed
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
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
230
231
232
233
234
235
236
  bathDC.readMsh("./Bath/"+simMesh+"_bath_smooth/"+simMesh+"_bath_smooth.idx")
else:
  print(colored('**** ERROR: Bathymetry file does not exist! ****',"red"))


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)

thomas's avatar
thomas committed
237
# Load reef maps
thomas's avatar
thomas committed
238
239
reefsShallowDC = dgDofContainer(groups,1)
reefsDeepDC = dgDofContainer(groups,1)
thomas's avatar
thomas committed
240
241
reefsShallowDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapShallowFilenameTrunk+"_"+simMesh
reefsDeepDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapDeepFilenameTrunk+"_"+simMesh
thomas's avatar
thomas committed
242

thomas's avatar
thomas committed
243
#1. Load reef maps from data if they don't already exist
thomas's avatar
thomas committed
244
245
246
247
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 stereoFlag == 0:
thomas's avatar
thomas committed
248
249
250
    reefsObj = reefs("./Data/"+reefMapShallowFilename, lonLatDegrees, 1, EarthRadius)
  #else:
    #reefsObj = reefs("./Data/"+reefMapShallowFilename, stereoToLonLatDegrees, 1, EarthRadius)
thomas's avatar
thomas committed
251
252
253
254
255
256
257
258
259
260
261
262
263
  # 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 stereoFlag == 0:
thomas's avatar
thomas committed
264
265
266
    reefsObj = reefs("./Data/"+reefMapDeepFilename, lonLatDegrees, 1, EarthRadius)
  #else:
    #reefsObj = reefs("./Data/"+reefMapDeepFilename, stereoToLonLatDegrees, 1, EarthRadius)
thomas's avatar
thomas committed
267
268
269
270
271
272
273
274
275
  # 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
276
#2. Load reefs data
thomas's avatar
thomas committed
277
278
279
280
281
282
283
284
285
286
287
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
288
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
289
290
291
292
293
294
295
296
297
298
299
300
301
302
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
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
340
reefsDualFC = functionC(glib, "overwrite2", 1, [bathDC.getFunction(), reefsShallowDC.getFunction(), reefsDeepDC.getFunction(), nbOfShallowReefsFConst])
thomas's avatar
thomas committed
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
reefsDualDC.interpolate(reefsDualFC)
reefsDualDC.exportMsh("LPT_Functions/REEFS_RAW_LPT_DUAL_"+simMesh,0,0)

print(colored('Seeding particles ...',"red"))
# ***SEED OVER ALL REEFS***
initCM = fullMatrixDouble( nbOfReefs, nbOfReefs )
particles = particleArray(groups, reefsDualDC, initCM, reefsToSeed, concentration, seedOverWhichReefs, nbOfShallowReefs)
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)
#cartCoords = [0]*3
thomas's avatar
thomas committed
357
#geoCoordTo3DCartesian(-19.9400,149.1100,cartCoords,EarthRadius)
thomas's avatar
thomas committed
358
#projectPoint(cartCoords)
thomas's avatar
thomas committed
359
360
361
#particles = particleArray(groups, 0, cartCoords[0], cartCoords[1], 0, 0)
#for location in range(0,10):
  #particles.addParticles(groups, 1000, cartCoords[0], cartCoords[1], 0, location)
thomas's avatar
thomas committed
362
363
364
365
366
367
#particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
#particles.printPositions(LP_OutputDir+'/seededParticles','dat',0)


print ''
print '-----------------------'
thomas's avatar
thomas committed
368
print 'HYDRODYNAMICS PARAMETERS:'
thomas's avatar
thomas committed
369
370
371
372
373
374
375
376
377
378
379
380
381
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
382
383
384
print ''
print '-----------------------'
print 'BIOLOGICAL PARAMETERS:'
thomas's avatar
thomas committed
385
386
print 'Mortality flag:', mortParams(0,0), '  (0: OFF, 1: CONST, 2: WEIBULL)'
if (int(mortParams(0,0)) == 1):
thomas's avatar
thomas committed
387
  print 'Mortality rate (const):', mortParams(0,1), '(per day)'
thomas's avatar
thomas committed
388
389
390
elif (int(mortParams(0,0)) == 2):
  print 'Weibull lambda:', mortParams(0,2), ', Weibull nu:', mortParams(0,3)
print '-----------------------'
thomas's avatar
thomas committed
391
print 'Settling flag: ', settleParams(0,0), '  (0: OFF, 1: CONST, 2: TRIANGLE, 3: PROGRESSIVE)'
thomas's avatar
thomas committed
392
393
394
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
395
396
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
397
398
399
400
401
402
403
404
405
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
406
print ''
thomas's avatar
thomas committed
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
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 ''

simSolution = dgDofContainer(groups, 3)
simSolutionUnder = dgDofContainer(groups, 3)
simSolutionOver = dgDofContainer(groups, 3)

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

#Create particle tracker object
print 'Initialising particle tracker ...'
particleTracker = dgParticleTracker2D(groups, particles, bathDC, diffDC, reefsDualDC, distToReefsDC, directionToReefsDC, settleParams, mortParams, swimParams, nbOfShallowReefs, LP_TotIter, LP_dt, LP_OutputDir)

startcpu = time.clock()
#Keep track of the solutions loaded into the three DCs:
loadedSolutions = [-1] * 3 #[0]:simSolution, [1]:simSolutionOver, [2]:simSolutionUnder
thomas's avatar
thomas committed
429
print 'Beginning loop ...'
thomas's avatar
thomas committed
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
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
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) | 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
522
    if (int(settleParams(0,0)) != 0):
thomas's avatar
thomas committed
523
      particles.printPositions(LP_OutputDir+'/particlesSettled_%06d'%( n ), 'pos', 1)
thomas's avatar
thomas committed
524
    if (int(mortParams(0,0)) != 0):
thomas's avatar
thomas committed
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
      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)