LPT_Caller_DualReefMap.py 23.7 KB
Newer Older
thomas's avatar
thomas committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
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
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
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
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
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
522
523
524
525
526
527
528
529
530
531
532
533
534
535
# 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
# 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 ''

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"

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

# Define simOutput directory & particle concentration

#LOCAL:
scratchFolder = 'Saved_Data/'
#scratchFolder = '/media/BackupDrive/Data/GBR_Hydrodynamics/'
simOutputDirName = 'Output10K_4wks_windsOn_allTides_16112000'
#LP_OutputDir = scratchFolder + simOutputDirName

LP_OutputDir = 'temp'#scratchFolder + simOutputDirName +'/LPT_Output/NoSett_NoMort'

#EVERYONE:
simOutputDir = scratchFolder + simOutputDirName
##
concentration = 2.0  #No. of particles to seed per kmsq
reefMapShallowFilename = "SGBR_ReefMap_2121_Solid9200.dat"
reefMapDeepFilename = "SGBR_DeepReefMap_2121_Solid9200.dat"

stereoFlag = 1
if stereoFlag == 1:
  EarthRadius=6371e3
else:
  EarthRadius=0

##

# SELECT REEF TYPES TO SEED AND SETTLE OVER
#0: shallow reefs, 1: deep reefs, 2: both
seedOverWhichReefs = 0
settleOverWhichReefs = 0
##

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()
simMesh = "gbr20K_S_sgbr"#"gbr10K"# #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)]

print(colored('Loading particle-tracking parameters ...',"red"))
LP_dt = 120.                                            #(double) in seconds
LP_TotTime = 3*7*24.0                                  #(double) in hours
LP_Export = 1#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 PARAM
#  -> 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, 0.0)
#[1]: (max if Weibull) mortality rate (PER DAY)
  #NB: 0.09 d-1 for G. retiformis (Connolly 2k11), 0.07 for A. Gemmifera
  #    Reef fish: 0.18 (James et al. 2002)
mortParams.set(0,1, 0.09)
#[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)
#[1]: t_preComp (seconds) = t_timeLag (for [3])
settleParams.set(0,1, 0.0)
#settleParams.set(0,1, 5.0*24.0*3600.0)
#[2]: t_postComp
settleParams.set(0,2, 25.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)
settleParams.set(0,5, 0.5)
#[6]: (for type-3 only) rate of loss of competence (d-1)
settleParams.set(0,6, 0.0)
#[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)

#TODO: convert swimming speeds to stereo

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")):
  print(colored('Importing bathymetry.',"red"))
  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)

# Load reef map(s).
reefsShallowDC = dgDofContainer(groups,1)
reefsDeepDC = dgDofContainer(groups,1)
reefsShallowDCfilename = "LPT_Functions/REEFS_RAW_LPT_SHALLOW_"+simMesh
reefsDeepDCfilename = "LPT_Functions/REEFS_RAW_LPT_DEEP_"+simMesh

#1. Load data if it doesn't already exist
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:
    reefsObj = reefs("./Data/"+reefMapShallowFilename, lonLatDegrees, 1)
  else:
    reefsObj = reefs("./Data/"+reefMapShallowFilename, stereoToLonLatDegrees, 1, EarthRadius)
  # 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:
    reefsObj = reefs("./Data/"+reefMapDeepFilename, lonLatDegrees, 1)
  else:
    reefsObj = reefs("./Data/"+reefMapDeepFilename, stereoToLonLatDegrees, 1, EarthRadius)
  # 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 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, ")"

# Record simulation info to file
LP_SimInfo = ['LPT starts '+str((LP_StartTime-simStartTime)/3600)+' hrs after beginning of hydrodynamics simulation\n', 'LPT simulation length: '+str(LP_TotTime)+' hrs\n', 'Exporting every '+str(LP_Export)+' LPT steps\n', 'LPT dt: '+ str(LP_dt)+' s\n', 'No. of simulation iterations per LP iteration: '+str(simItersPerLPIter)+'\n', 'Mortality type:'+str(mortParams(0,0))+'\n', '(max) Mortality rate:'+str(mortParams(0,1)/(24.0*3600.0))+'(per day)\n','Settling type:'+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()

# 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)
reefsDualFC = functionC(glib, "overwrite2", 1, [reefsShallowDC.getFunction(), reefsDeepDC.getFunction(), nbOfShallowReefsFConst])
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
#geoCoordToCartesian(-19.9400,149.1100,cartCoords)
#projectPoint(cartCoords)
#particles = particleArray(groups, 1000, cartCoords[0], cartCoords[1], 0)
#particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
#particles.printPositions(LP_OutputDir+'/seededParticles','dat',0)


print ''
print '-----------------------'
print 'SIMULATION 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 'Mortality flag:', mortParams(0,0), '  (0: OFF, 1: CONST, 2: WEIBULL)'
if (int(mortParams(0,0)) == 1):
  print 'Mortality rate (max):', 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)'
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)'
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 '-----------------------'

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

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)
    if (int(settleParams(0,0)) == 1):
      particles.printPositions(LP_OutputDir+'/particlesSettled_%06d'%( n ), 'pos', 1)
    if (int(mortParams(0,0)) == 1):
      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)