LPT_Caller_DualReefMap.py 25.4 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
  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":
35
36
    cpuExt = "/CPU"+sys.argv[2]
    time.sleep(float(sys.argv[2])*2) #Avoid simultaneous simulations crashing each other
thomas's avatar
thomas committed
37
38
39

# Define simOutput directory & particle concentration

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

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


##
concentration = 2.0  #No. of particles to seed per kmsq
49
50
#reefMapShallowFilenameTrunk = "SGBR_ReefMap_2121_Solid9200.dat"
#reefMapDeepFilenameTrunk = "SGBR_DeepReefMap_2121_Solid9200.dat"
thomas's avatar
thomas committed
51
52
53
54
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
55

56
# Radius: 6.37101e6 for old meshes (sma_zoomed), else 6.371e6 (used for reefsObj and reefsDistanceObj)
thomas's avatar
thomas committed
57
EarthRadius=6.37101e6
thomas's avatar
thomas committed
58
# Set whether 1: to use PROJ projection library (set parameters in GBR.cc) or 0: to use old projection
59
60
61
usePROJ = 1
# If using PROJ, set projection systems:
lonlat_projection = "+proj=latlong +ellps=WGS84"
thomas's avatar
thomas committed
62
mesh_projection = "+proj=utm +ellps=WGS84 +zone=55 +south"
thomas's avatar
thomas committed
63
64
65
66
##

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

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

XYZ = function.getCoordinates()

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)
123
  #G. retiformis: 0.09, A. gemmifera: 0.067
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
#[2]: Weibull only (else unused): lambda parameter <-- (Connolly & Baird 2010)lambda = 1/b(Pinder et al. 1978)
127
#A. millepora: 0.043, A. valida: 0.019, P. daedalea: 0.060
thomas's avatar
thomas committed
128
129
mortParams.set(0,2, 0.043)
#[3]: Weibull only (else unused): nu parameter <-- (Connolly & Baird 2010)nu = c(Pinder et al. 1978)
130
#A. millepora: 0.57, A. valida: 0.46, P. daedalea: 0.72
thomas's avatar
thomas committed
131
132
133
134
135
136
137
138
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
140
# 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
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)
149
# 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
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)
152
# 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
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
168
169
#[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)

print(colored('Loading model ...',"red"))
model = GModel()
170
171
if usePROJ == 0:
  print " (Using old projection)"
thomas's avatar
thomas committed
172
173
  model.load("./Meshes/"+simMesh+"_tan.msh")
else:
174
175
  print " (Using PROJ for projections)"
  model.load("./Meshes/"+simMesh+"_utm.msh")
thomas's avatar
thomas committed
176
177
178
179
180
181
182
183
184
185
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
186
  print(colored('Importing bathymetry ...',"red"))
thomas's avatar
thomas committed
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
  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);
    }
  }
}

"""
218
tmpCLib = "./tempCLib.dylib"
thomas's avatar
thomas committed
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
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
235
# Load reef maps
thomas's avatar
thomas committed
236
237
reefsShallowDC = dgDofContainer(groups,1)
reefsDeepDC = dgDofContainer(groups,1)
thomas's avatar
thomas committed
238
239
reefsShallowDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapShallowFilenameTrunk+"_"+simMesh
reefsDeepDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapDeepFilenameTrunk+"_"+simMesh
thomas's avatar
thomas committed
240

thomas's avatar
thomas committed
241
#1. Load reef maps from data if they don't already exist
242
243
244
UTMProjector = None
if usePROJ == 1:
  UTMProjector = slimProjector(mesh_projection, lonlat_projection)
thomas's avatar
thomas committed
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)
248
249
250
251
  if usePROJ == 0:
    reefsObj = reefsRaw_oldProj("./Data/"+reefMapShallowFilename, EarthRadius)
  else:
    reefsObj = reefsRaw("./Data/"+reefMapShallowFilename, UTMProjector)
thomas's avatar
thomas committed
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)
264
265
266
267
  if usePROJ == 0:
    reefsObj = reefsRaw_oldProj("./Data/"+reefMapDeepFilename, EarthRadius)
  else:
    reefsObj = reefsRaw("./Data/"+reefMapDeepFilename, UTMProjector)
thomas's avatar
thomas committed
268
269
270
271
272
273
274
275
276
  # 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
277
#2. Load reefs data
thomas's avatar
thomas committed
278
279
280
281
282
283
284
285
286
287
288
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
289
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
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
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
341
reefsDualFC = functionC(glib, "overwrite2", 1, [bathDC.getFunction(), reefsShallowDC.getFunction(), reefsDeepDC.getFunction(), nbOfShallowReefsFConst])
thomas's avatar
thomas committed
342
343
344
345
346
reefsDualDC.interpolate(reefsDualFC)
reefsDualDC.exportMsh("LPT_Functions/REEFS_RAW_LPT_DUAL_"+simMesh,0,0)

print(colored('Seeding particles ...',"red"))
# ***SEED OVER ALL REEFS***
thomas's avatar
thomas committed
347
348
349
350
351
352
353
354
#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' )
thomas's avatar
thomas committed
355
356

# ***SEED AT A POINT*** (NB: Look at method used to print results if changing seeding method)
thomas's avatar
thomas committed
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375

UTMProjectorLonLatToUTM = slimProjector(lonlat_projection, mesh_projection)
reefPointsCart = [ [0]*4 for i in range(len(reefPointsLatLon)) ]
pj_latlong = pj_init_plus(lonlat_projection)
pj_merc = pj_init_plus(mesh_projection)
for i in range(0,len(reefPointsLatLon)):
  reefPointsCart[i][3] = reefPointsLatLon[i][2]
  src_points = fullMatrixDouble(2,2)
  dest_points = fullMatrixDouble(2,2)
  src_points.set(0,0,reefPointsLatLon[i][1]/180.0*pi)
  src_points.set(1,1,reefPointsLatLon[i][0]/180.0*pi)
  UTMProjectorLonLatToUTM.projectPoint(src_points, dest_points)
  reefPointsCart[i][0] = dest_points(0,0)
  reefPointsCart[i][1] = dest_points(1,1)

particles = particleArray(groups, 0, 0, 0, 0, 0)
for location in range(0,len(reefPointsCart)):
  particles.addParticles(groups, 100, reefPointsCart[location][0], reefPointsCart[location][1], 0, reefPointsCart[location][3])
particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
thomas's avatar
thomas committed
376
377
378
379
380
#particles.printPositions(LP_OutputDir+'/seededParticles','dat',0)


print ''
print '-----------------------'
thomas's avatar
thomas committed
381
print 'HYDRODYNAMICS PARAMETERS:'
thomas's avatar
thomas committed
382
383
384
385
386
387
388
389
390
391
392
393
394
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
395
396
397
print ''
print '-----------------------'
print 'BIOLOGICAL PARAMETERS:'
thomas's avatar
thomas committed
398
399
print 'Mortality flag:', mortParams(0,0), '  (0: OFF, 1: CONST, 2: WEIBULL)'
if (int(mortParams(0,0)) == 1):
thomas's avatar
thomas committed
400
  print 'Mortality rate (const):', mortParams(0,1), '(per day)'
thomas's avatar
thomas committed
401
402
403
elif (int(mortParams(0,0)) == 2):
  print 'Weibull lambda:', mortParams(0,2), ', Weibull nu:', mortParams(0,3)
print '-----------------------'
thomas's avatar
thomas committed
404
print 'Settling flag: ', settleParams(0,0), '  (0: OFF, 1: CONST, 2: TRIANGLE, 3: PROGRESSIVE)'
thomas's avatar
thomas committed
405
406
407
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
408
409
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
410
411
412
413
414
415
416
417
418
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
419
print ''
thomas's avatar
thomas committed
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
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 ...'
thomas's avatar
thomas committed
437
438
439
440
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
441
442
443
444

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
445
print 'Beginning loop ...'
thomas's avatar
thomas committed
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
536
537
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
538
    if (int(settleParams(0,0)) != 0):
thomas's avatar
thomas committed
539
      particles.printPositions(LP_OutputDir+'/particlesSettled_%06d'%( n ), 'pos', 1)
thomas's avatar
thomas committed
540
    if (int(mortParams(0,0)) != 0):
thomas's avatar
thomas committed
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
      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)