LPT_Caller_DualReefMap.py 29.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
#      --> 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
##

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

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

thomas's avatar
thomas committed
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
# 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'
thomas's avatar
thomas committed
217

thomas's avatar
thomas committed
218
bathDC = dgDofContainer(groups, 1);
thomas's avatar
thomas committed
219
if(os.path.exists("./Bath/"+simMesh+"_bath_smooth/"+simMesh+"_bath_smooth.idx")):
thomas's avatar
thomas committed
220
  print(colored('Importing bathymetry ...',"red"))
221
  bathDC.importIdx("./Bath/"+simMesh+"_bath_smooth/"+simMesh+"_bath_smooth.idx")
thomas's avatar
thomas committed
222
223
else:
  print(colored('**** ERROR: Bathymetry file does not exist! ****',"red"))
thomas's avatar
thomas committed
224
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
225

thomas's avatar
thomas committed
226
#------- SET UP DIFFUSIVITY ----------------------------------------------------
thomas's avatar
thomas committed
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
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);
    }
  }
}

"""
253
tmpCLib = "./tempCLib.dylib"
thomas's avatar
thomas committed
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
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
269
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
270

thomas's avatar
thomas committed
271
#------- LOAD REEF MAPS --------------------------------------------------------
thomas's avatar
thomas committed
272
273
reefsShallowDC = dgDofContainer(groups,1)
reefsDeepDC = dgDofContainer(groups,1)
thomas's avatar
thomas committed
274
275
reefsShallowDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapShallowFilenameTrunk+"_"+simMesh
reefsDeepDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapDeepFilenameTrunk+"_"+simMesh
thomas's avatar
thomas committed
276

thomas's avatar
thomas committed
277
#1. Load reef maps from data if they don't already exist
278
279
280
UTMProjector = None
if usePROJ == 1:
  UTMProjector = slimProjector(mesh_projection, lonlat_projection)
thomas's avatar
thomas committed
281
282
283
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)
284
285
286
287
  if usePROJ == 0:
    reefsObj = reefsRaw_oldProj("./Data/"+reefMapShallowFilename, EarthRadius)
  else:
    reefsObj = reefsRaw("./Data/"+reefMapShallowFilename, UTMProjector)
thomas's avatar
thomas committed
288
289
290
291
292
293
294
295
296
297
298
299
  # 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)
300
301
302
303
  if usePROJ == 0:
    reefsObj = reefsRaw_oldProj("./Data/"+reefMapDeepFilename, EarthRadius)
  else:
    reefsObj = reefsRaw("./Data/"+reefMapDeepFilename, UTMProjector)
thomas's avatar
thomas committed
304
305
306
307
308
309
310
311
312
  # 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
313
#2. Load reefs data
thomas's avatar
thomas committed
314
315
316
317
318
319
320
321
322
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, ")"
thomas's avatar
thomas committed
323
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
324

thomas's avatar
thomas committed
325
#------- WRITE SIMULATION INFO TO FILE -----------------------------------------
thomas's avatar
thomas committed
326
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
327
328
329
LP_SimInfoFile = open(LP_OutputDir+'/LP_runInfo.txt','w')
LP_SimInfoFile.writelines(LP_SimInfo)
LP_SimInfoFile.close()
thomas's avatar
thomas committed
330
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
331

thomas's avatar
thomas committed
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
#------- LOAD MAPS OF DISTANCE AND DIRECTION TO NEAREST REEF -------------------
distToReefsDC = None
directionToReefsDC = None
reefsDualDC = None
nbOfShallowReefsFConst = None
if swimParams(0,0) > 0.1:
  print(colored('Calculating distances to nearest shallow reefs ...', "red"))
  distToReefsDC = dgDofContainer(groups, 1)
  directionToReefsDC = dgDofContainer(groups, 3)

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

nbOfShallowReefsFConst = functionConstant(nbOfShallowReefs)
reefsDualDC = dgDofContainer(groups, 1)
thomas's avatar
thomas committed
358
reefsDualFC = functionC(glib, "overwrite2", 1, [bathDC.getFunction(), reefsShallowDC.getFunction(), reefsDeepDC.getFunction(), nbOfShallowReefsFConst])
thomas's avatar
thomas committed
359
360
reefsDualDC.interpolate(reefsDualFC)
reefsDualDC.exportMsh("LPT_Functions/REEFS_RAW_LPT_DUAL_"+simMesh,0,0)
thomas's avatar
thomas committed
361
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
362

thomas's avatar
thomas committed
363
#--------- SEED PARTICLES ------------------------------------------------------
thomas's avatar
thomas committed
364
print(colored('Seeding particles ...',"red"))
thomas's avatar
thomas committed
365
366
367
368
369
370
371
372
373
374
375

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
376
# ***SEED OVER ALL REEFS***
thomas's avatar
thomas committed
377
378
379
380
381
382
383
384
385
386
initCM = fullMatrixDouble( nbOfReefs, nbOfReefs )
particles = particleArray()
particles.addParticlesWithReefMap(groups, reefsDualDC, initCM, reefsToSeed, concentration, minNbParticlesPerReef, seedOverWhichReefs, nbOfShallowReefs)
#particles.addParticlesWithDofContainer(groups, reefsDualDC, concentration, 200)
particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
dummyPrinterICM = connectivityMatrixPrint()
if settleParams(0,0) == 0:
  dummyPrinterICM.CMPrinter(initCM, LP_OutputDir+'/connectivityMatrix_initial' )
else:
  dummyPrinterICM.CMPrinter(initCM, LP_OutputDir+'/settlementMatrix_initial' )
thomas's avatar
thomas committed
387
# ******
thomas's avatar
thomas committed
388
389

# ***SEED AT A POINT*** (NB: Look at method used to print results if changing seeding method)
thomas's avatar
thomas committed
390
391
392
393
394
395
396
397
398
399
400
401
402
403
#UTMProjectorLonLatToUTM = slimProjector(lonlat_projection, mesh_projection)
#reefPointsLatLon = [ [-14.7,145.4,0], [-14.7,145.4,1] ]
#reefPointsCart = [ [0]*4 for i in range(len(reefPointsLatLon)) ]
#pj_latlong = pj_init_plus(lonlat_projection)
#pj_merc = pj_init_plus(mesh_projection)
#for i in range(0,len(reefPointsLatLon)):
  #reefPointsCart[i][3] = reefPointsLatLon[i][2]
  #src_points = fullMatrixDouble(2,2)
  #dest_points = fullMatrixDouble(2,2)
  #src_points.set(0,0,reefPointsLatLon[i][1]/180.0*pi)
  #src_points.set(1,1,reefPointsLatLon[i][0]/180.0*pi)
  #UTMProjectorLonLatToUTM.projectPoint(src_points, dest_points)
  #reefPointsCart[i][0] = dest_points(0,0)
  #reefPointsCart[i][1] = dest_points(1,1)
thomas's avatar
thomas committed
404

thomas's avatar
thomas committed
405
406
407
408
409
#particles = particleArray()
#for location in range(0,len(reefPointsCart)):
  #particles.addParticlesAtPoint(groups, 100, reefPointsCart[location][0], reefPointsCart[location][1], 0, reefPointsCart[location][3])
#particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
#particles.printPositions(LP_OutputDir+'/seededParticles','dat',0)
thomas's avatar
thomas committed
410
411
# ******
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
412

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

thomas's avatar
thomas committed
464
#------- SET UP HYDRO DOFCONTAINERS --------------------------------------------
thomas's avatar
thomas committed
465
466
467
simSolution = dgDofContainer(groups, 3)
simSolutionUnder = dgDofContainer(groups, 3)
simSolutionOver = dgDofContainer(groups, 3)
thomas's avatar
thomas committed
468
469
#Keep track of the solutions loaded into the three DCs:
loadedSolutions = [-1] * 3 #[0]:simSolution, [1]:simSolutionOver, [2]:simSolutionUnder
thomas's avatar
thomas committed
470
471
472

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

thomas's avatar
thomas committed
475
#------ SET UP PARTICLE TRACKER OBJECT -----------------------------------------
thomas's avatar
thomas committed
476
477
#Create particle tracker object
print 'Initialising particle tracker ...'
thomas's avatar
thomas committed
478
479
480
particleTracker = dgParticleTracker2D(groups, particles, bathDC, diffDC, LP_TotIter, LP_dt, LP_OutputDir)
particleTracker.addMortality(mortParams)
particleTracker.addSettling(settleParams, reefsDualDC, nbOfShallowReefs)
thomas's avatar
thomas committed
481
#particleTracker.addSwimming(swimParams, reefsDualDC, distToReefsDC, directionToReefsDC)
thomas's avatar
thomas committed
482
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
483

thomas's avatar
thomas committed
484
#-------- PARTICLE LOOP --------------------------------------------------------
thomas's avatar
thomas committed
485
print 'Beginning loop ...'
thomas's avatar
thomas committed
486
487
488
489
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
490
for n in range(1,LP_TotIter+1):
thomas's avatar
thomas committed
491
492
493
494
495
496
497
  #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
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
  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 )+').'
521
      simSolution.importIdx(simOutputDir+'/mr-gbr/mr-gbr-%06d'%( iterNumberWanted )+'.idx')
thomas's avatar
thomas committed
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
      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
548
      simSolutionUnder.importIdx(simOutputDir+'/mr-gbr/mr-gbr-%06d'%( iterNumberWanted-l )+'.idx')
thomas's avatar
thomas committed
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
      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
573
      simSolutionOver.importIdx(simOutputDir+'/mr-gbr/mr-gbr-%06d'%( iterNumberWanted+m )+'.idx')
thomas's avatar
thomas committed
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
      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
589
    if (int(settleParams(0,0)) != 0):
thomas's avatar
thomas committed
590
      particles.printPositions(LP_OutputDir+'/particlesSettled_%06d'%( n ), 'pos', 1)
thomas's avatar
thomas committed
591
    if (int(mortParams(0,0)) != 0):
thomas's avatar
thomas committed
592
593
594
595
596
597
598
599
600
601
602
603
      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
604
#-------------------------------------------------------------------------------
thomas's avatar
thomas committed
605
606
607

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