dgParticleTracker2D.h 15 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts 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
#ifndef _DG_PARTICLE_TRACKER_2D_H_
#define _DG_PARTICLE_TRACKER_2D_H_

#include <fstream>
#include <iomanip>
#include <sstream>
#include <iostream>
#include <vector>
#include"math.h"
#include"dgDofContainer.h"
#include"dgFunctionEvaluator.h"
#include"dgGroupOfElements.h"
#include"dgNeighbourMap.h"
#include "functor.h"

double randoraw();
void rando(double *randoo, int M);
void Inverse(double Matrix[3][3], double Inverse[3][3]);
void Transpose(double Matrix[3][3]);
void projection(double* v1, double* v2, double* v1_proj_on_v2);
double projectionMag(double* v1, double* v2);

class particle {
private:
  double _xcoord, _ycoord;
  int _iter_seeded;     //iteration at which particle enters domain
  int _competency;      //0: pre-competent, 1: competent, 2: post-competent
  int _int_state;       //0: alive, 1: settled, 2: dead, 3: out of domain
  int _source_reef, _sink_reef;
  //Remember element and group Id particle is in:
  int _groupId, _elementId;
  
public:
  // Default constructor
  particle()
  {
    _xcoord = 0.0;
    _ycoord = 0.0;
    _int_state = 0;
    _source_reef = -1;
    _sink_reef = -1;
    _groupId = -1;
    _elementId = -1;
    _competency = 0;
  }
  
  // Cartesian constructor which SETS mesh element Id's
  particle(double x0, double y0, int iGroup, int iElement, int status, int iter)
  {
    _xcoord = x0;
    _ycoord = y0;
    _int_state = status;
    _source_reef = -1;
    _sink_reef = -1;
    _groupId = iGroup;
    _elementId = iElement;
    _competency = 0;
    _iter_seeded = iter;
  }
  
  // Cartesian constructor which DOESN'T set mesh element Id's
  particle(double x0, double y0, int status, int iter)
  {
    _xcoord = x0;
    _ycoord = y0;
    _int_state = status;
    _source_reef = -1;
    _groupId = -1;
    _elementId = -1;
    _competency = 0;
    _iter_seeded = iter;
  }
  
  // Cartesian constructor which sets source reef AND mesh element Id's
  particle(double x0, double y0, int iGroup, int iElement, int srcReef, int status, int iter)
  {
    _xcoord = x0;
    _ycoord = y0;
    _int_state = status;
    _source_reef = srcReef;
    _sink_reef = -1;
    _groupId = iGroup;
    _elementId = iElement;
    _competency = 0;
    _iter_seeded = iter;
  }
  
  // Cartesian constructor which sets source reef but NOT mesh element Id's (try not to use this, doesn't store element we're in)
  particle(double x0, double y0, int srcReef, int status, int iter)
  {
    _xcoord = x0;
    _ycoord = y0;
    _int_state = status;
    _source_reef = srcReef;
    _sink_reef = -1;
    _groupId = -1;
    _elementId = -1;
    _competency = 0;
    _iter_seeded = iter;
  }
  
  //Access methods:
  double x()
  {
    return _xcoord;
  }
  
  double y()
  {
    return _ycoord;
  }
  
  int getState()
  {
    return _int_state;
  }

  int getSource() 
  {
    return _source_reef;
  }
  
  int getSink() 
  {
    return _sink_reef;
  }
  
  void getElementId(int *iGroup, int *iElement)
  {
    *iGroup = _groupId;
    *iElement = _elementId;
  }
  
  int getCompetency()
  {
    return _competency;
  }
  
  int getIterSeeded()
  {
    return _iter_seeded;
  }
  
  void setPosition(double x, double y)
  {
    _xcoord=x;
    _ycoord=y;
  }
  
  void setState(int status)
  {
    _int_state = status;
  }
  
  void setSource(int source) 
  {
    _source_reef = source;
  }
  
  void setSink(int sink) 
  {
    _sink_reef = sink;
  }
  
  void setElementId(int iGroup, int iElement)
  {
    _groupId = iGroup;
    _elementId = iElement;
  }
  
  void setCompetency(int neoComp)
  {
    _competency = neoComp;
  }
  
};

class particleArray {
private:
  std::vector<particle> _particle_array;
  int _alive, _settled, _dead, _outOfDomain, _aliveAndCompetent;
  int _totalSeeded;
183
  std::vector<std::string> _openBoundaries;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
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
  
public:
  /** Default constructor: use this in conjunction with addParticles to seed particles. */
  particleArray()
  {
    _alive = 0;
    _settled = 0;
    _dead = 0;
    _outOfDomain = 0;
    _aliveAndCompetent = 0;
    _totalSeeded = 0;
  }
  
  /** Seed particles at a given point. */
  void addParticlesAtPoint(dgGroupCollection* _group, int M, double x0, double y0, int iter, int status, int source_reef);
  
  /** Seed particles using a DofContainer. 
   *    minNbParticlesPerHabitat gives the minimum number of particles to seed per habitat (different habitats are identified by having different integer values of the dof container in each element) */
  void addParticlesWithDofContainer(dgGroupCollection* groups, dgDofContainer* DC, double concentration, int minNbParticlesPerHabitat, int iter);
  
  /** Seed particles using reef map DofContainer.*/
  void addParticlesWithReefMap(dgGroupCollection* groups, dgDofContainer* reefsDC, fullMatrix<double>* initCM, fullMatrix<double>* reefsToSeed, double concentration, int minNbParticlesPerReef, int seedOverWhichReefs, int nbOfShallowReefs, int iter);
  
  /** Get/set particle stats. */
  void setStats(int *stats) {
    _alive = stats[0];
    _settled = stats[1];
    _dead = stats[2];
    _outOfDomain = stats[3];
    _aliveAndCompetent = stats[4];
  }
  
  void getStats(int *stats) {
    stats[0] = _alive;
    stats[1] = _settled;
    stats[2] = _dead;
    stats[3] = _outOfDomain;
    stats[4] = _aliveAndCompetent;
  }  
  
  int printNbParticles() {
    return static_cast<int>(_particle_array.size());
  }
  
  /** Access a given particle. */
  particle* getParticle(int m)
  {
    return &_particle_array[m];
  }

  double getParticle_x(int m)
  {
    return _particle_array[m].x();
  }
  
  double getParticle_y(int m)
  {
    return _particle_array[m].y();
  } 
  
  /** Print positions of all particles.
   *  Generates 1 file per call with particles colour-coded by source reef.
   *  Only records particles of state 'status'
   *    filetype = dat or pos 
   */
  void printPositions(const std::string filename, const std::string filetype, int status);
  
  /** Print positions of all particles. 
   *  1 file per particle, colour-coded by time-step.
   *  filetype = dat or pos
   */
  void printPositionsFollow(const std::string filename, const std::string filetype, int n, int N, double dt);
  
  /** Print position of one particle. 
   *  Colour-coded by time-step.
   *  filetype = dat or pos
   */
  void printPositionsFollowSingle(const std::string filename, const std::string filetype, int n, int N, double dt, int particleId);

  /** Method to save state of a particle array for later recovery.
  * Output file:
  * Line1: Header info:
  *          _totalSeeded _alive _settled _dead _outOfDomain _aliveAndCompetent
  * Lines2->End: Particle info:
  *          x y state sourceReef competency groupId elementId
  * */
  void saveState(const std::string &filename);
271
272
273

  void setOpenBoundaries(std::vector<std::string> tags) {_openBoundaries = tags;}
  bool isOpenBoundary(std::string tag) const;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
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
  
  /** Method to load state of a particle array from file generated using saveState(). */
  void loadState(const std::string &filename);
  
  /** Method which returns a list of all the mesh elements containing any particles inside */
  fullMatrix<double> getElementsContainingParticles(dgGroupCollection* groups);
  
  /** Method to calculate and return a connectivity matrix considering only particles which have settled on specified elements (element IDs to consider supplied as arg in a fullMatrix).
   * Arg1: Empty CMatrix to put results in (must have the same size as the full CMatrix)
   * Arg2: List of element Ids of elements to consider 
   * */
  void getCMatrixForListOfElements(fullMatrix<double>* CMatrix, const fullMatrix<double>* listOfElements);
};

class connectivityMatrixPrint {
private:
  
public:
  /** Constructor */
  connectivityMatrixPrint() {
  }
  /** Method to print contents of matrix */
  void CMPrinter(fullMatrix<double>* FM, const std::string filename);
};

class dgParticleTracker2D {
private:
  double _dt;
  int _N;
  int _settleType, _mortType, _swimType, _windDragType, _settleOverWhichReefType, _nbOfShallowReefs;
304
  double _t_preComp, _t_postComp, _t_settMax, _settRate_max, _compAcquisitionRate, _compLossRate;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
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
  double _mortRate_max, _mortLambda, _mortNu;
  double _swimSpeed, _swimDist, _swimStart;
  double _windDragStart, _windDragStop, _windDragCoeff;
  dgGroupCollection *_group;
  particleArray *_ParticleArray;
  dgDofContainer* _bathymetry;
  dgDofContainer* _diffusivity;
  dgDofContainer* _reefs;
  dgDofContainer* _distToReefs;
  dgDofContainer* _directionToReefs;
  functor* _windUV;
  GModel* _gm;
  dgNeighbourMap* _neighbourMap;
  std::string _LP_OutputDir;
  bool _mortFlag, _settleFlag, _swimFlag, _connectivityExposureTimeFlag, _windFlag;
  
public:
  /** Constructor */
  dgParticleTracker2D(dgGroupCollection* group, particleArray* ParticleArray, dgDofContainer* bathymetry, dgDofContainer* diffusivity, int N, double dt, const std::string LP_OutputDir) {
    _dt = dt;
    _N = N;
    _group = group;
    _ParticleArray = ParticleArray;
    _bathymetry = bathymetry;
    _diffusivity = diffusivity;
    _LP_OutputDir = LP_OutputDir;
    
    //By default, set mortality, settling, swimming, use of exposure time connectivity matrix, wind forcing to OFF
    _mortFlag = false;
    _settleFlag = false;
    _swimFlag = false;
    _connectivityExposureTimeFlag = false;
    _windFlag = false;
338
    _settleOverWhichReefType = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
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
    
    //Load model:
    _neighbourMap = new dgNeighbourMap(*_group);
    
    //Put a dummy dgDofContainer into containers used for swimming (will be overwritten by addSwimming if used, else will remain unused)
    _reefs = bathymetry;
    _distToReefs = bathymetry;
    _directionToReefs = bathymetry;
    _windUV = bathymetry->getFunction();
    
    //Initialise output files
    std::ostringstream fname;
    std::ofstream outStats;
    //Define filename
    fname << _LP_OutputDir << "/particleStats.dat";
    //Open file (replace if already exists)
    std::string temp = fname.str();
    outStats.open( temp.c_str(), std::ios::out | std::ios::trunc );
    outStats.close();
  }
  
  ~dgParticleTracker2D() {
    delete _neighbourMap;
  }
  
  /** Activate particle mortality and load parameters */
365
366
367
368
369
370
  /*  mortParams :
   *    [0] : mortType 
   *        0 : no mortality
   *        1 : constant mortality rate given by [1] mortRateMax (by day)
   *        2 : weibullMortalityRate with [2] lambda  and [3] nu;
   */
371
  void addMortality(int type, double rateMax, double weibullLambda = 0, double weibullNu = 0) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
372
    _mortFlag = true;
373
    _mortType = type;
374
375
376
    if (type != 1) {
      Msg::Warning("Weibull mortality implementation may be wrong\n");
    }
377
378
379
    _mortRate_max = rateMax;
    _mortLambda = weibullLambda;   //(Connolly & Baird 2010)lambda = 1/b(Pinder et al. 1978)
    _mortNu = weibullNu;       //(Connolly & Baird 2010)nu = c(Pinder et al. 1978)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
380
381
382
  }
  
  /** Activate particle settling and load parameters */
383
384
385
386
387

  /* reefs :
   *    a dof container with the id of the reef by element (0 = no reef)*/
  /* flat settling all particles can settle between tPreComp and tPostComp  with probabylity by day given py settleRate*/
  void setFlatSettling(double settleRate, dgDofContainer* reefs, double tPreComp, double tPostComp){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
388
    _settleFlag = true;
389
390
391
392
    _settleType = 1;
    _settRate_max = settleRate;
    _t_preComp =  tPreComp;
    _t_postComp =  tPostComp;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
393
394
    _reefs = reefs;
  }
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

  /* all particle can settle with a probability that increases linearly between 0 in [1] tPreComp and [4] settleRateMax in [3] tSettleMax then decreases linearly to 0 in [2] tPostComp */
  void setLinearSettling(double maxRate, dgDofContainer* reefs, double tPreComp, double tPostComp, double tSettleMax) {
    _settleFlag = true;
    _settleType = 2;
    _settRate_max = maxRate;
    _t_preComp = tPreComp;
    _t_postComp = tPostComp;
    _t_settMax = tSettleMax;
    _reefs = reefs;
  }

  /* all particles settles when they are competents, competency is aquired and lost with flat rate [5] compAquisitionRate and [6] compLossRate*/
  void setCompetentSettling(dgDofContainer* reefs, double tPreComp, double tPostComp, double compAcquisitionRate, double compLossRate) {
    _settleFlag = true;
    _settleType = 3;
    _compAcquisitionRate = 1-pow(1-compAcquisitionRate, _dt/(24.0*3600.0));
    _compLossRate = 1-pow(1-compLossRate, _dt/(24.0*3600.0));
    _reefs = reefs;
    _t_preComp = tPreComp;
    _t_postComp = tPostComp;
  }

  /* type : 
   *        0: all reefs
   *        1: only shallow reefs (ie with idx < nbShallowReefs)
   *        2: only deep reefs (ie with idx > nbShallowReefs)
   */
  void setSettlingOverWhichReefs(int type, int nbShallowReefs) {
    _settleOverWhichReefType = type;
    _nbOfShallowReefs = nbShallowReefs;
  }

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
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
  /** Calculate connectivity matrix using exposure time method */
  void addConnectivityExposureTime(double t_preComp, double t_postComp) {
    _connectivityExposureTimeFlag = true;
    _settleFlag = false;
    _t_preComp = t_preComp;
    _t_postComp = t_postComp;
    printf("Connectivity by EXPOSURE TIME is enabled. Settling is therefore being disabled.\n");
  }
  
  /** Activate particle swimming and load parameters */
  void addSwimming(fullMatrix<double>* swimParams, dgDofContainer* reefs, dgDofContainer* distToReefs, dgDofContainer* directionToReefs) {
    _swimFlag = true;
    //TEMP: disactivate swimming if swimParams(0,0)=0.0
    if (swimParams->get(0,0) < 0.1) _swimFlag=false;
    //Load swimming parameters:
    _swimType = static_cast<int>(swimParams->get(0,0)+0.1);
    _swimSpeed = swimParams->get(0,1);
    _swimDist = swimParams->get(0,2);
    _swimStart = swimParams->get(0,3);
    if (_swimFlag == true) printf("swimType: %d, swimSpeed: %.2f(m/s), swimDist: %.2f(km), swimStart: %.2f(hrs)\n", _swimType, _swimSpeed, _swimDist, _swimStart);
    _reefs = reefs;
    _distToReefs = distToReefs;
    _directionToReefs = directionToReefs;
  }
  
  /** Activate direct wind forcing on particles (simulate wind drag at surface).
   ARGS: 1, windUV field
         2, windParams (0,0): type of wind forcing: 0: stop at t_stop defined in [2], 1: stop when particle is competent
                       (0,1): time to START direct wind forcing (seconds, relative to LPT_n=0)
                       (0,2): time to STOP direct wind forcing (seconds, relative to LPT_n=0)
                       (0,3): wind drag coefficient */
  void addWindForcing(functor* windUV, fullMatrix<double>* windParams) {
    _windFlag = true;
    _windUV = windUV;
    
    _windDragType = static_cast<int>(windParams->get(0,0) + 0.1);
    _windDragStart = windParams->get(0,1);
    _windDragStop = windParams->get(0,2);
    _windDragCoeff = windParams->get(0,3);
  }
  
  /** Method to move a particleArray on by a time-step */
  void particleMove(dgDofContainer* solution, fullMatrix<double>* connectivityMatrix, int n);
  /** Method to move a particleArray as seeds on by a time-step */
  void particleMoveSeeds(dgDofContainer* solution, dgDofContainer* windUV, fullMatrix<double>* connectivityMatrix, int n, double windcoef);
  /** Method to check if particle changes competency during time dt */
  int assessCompetency(int particleComp, double t);
  /** Method to check if particle settles or not during time dt */
  int isSettled(double t, int m);
  /** Method to check if particle dies or not during time dt */
  int isDead(double t);
};

#endif