Commit 12552a59 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

compilation ok

parent e130a038
......@@ -235,6 +235,7 @@ add_subdirectory(modules/perfectGas)
add_subdirectory(modules/shallowWater)
add_subdirectory(modules/slim3d)
add_subdirectory(modules/slimFunction)
add_subdirectory(modules/slimParticles)
add_subdirectory(modules/slimIce)
add_subdirectory(modules/subsurface)
add_subdirectory(modules/supra)
......
set(SRC
dgParticleTracker2D.cpp
)
dg_add_module(slimParticle "${SRC}" "")
#dg_add_swig_module(wave wave.i dgWave)
#dg_add_test_directory(tests axel.modave@gmail.com)
This diff is collapsed.
#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;
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);
/** 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;
double _t_preComp, _t_postComp, _t_settMax, _settRate_max, _compAcquisitionRate, _compAcquisitionProb, _compLossRate;
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;
//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 */
void addMortality(fullMatrix<double>* mortParams) {
_mortFlag = true;
_mortType = static_cast<int>(mortParams->get(0,0)+0.1);
_mortRate_max = mortParams->get(0,1);
_mortLambda = mortParams->get(0,2); //(Connolly & Baird 2010)lambda = 1/b(Pinder et al. 1978)
_mortNu = mortParams->get(0,3); //(Connolly & Baird 2010)nu = c(Pinder et al. 1978)
}
/** Activate particle settling and load parameters */
void addSettling(fullMatrix<double>* settleParams, dgDofContainer* reefs, int nbOfShallowReefs) {
_settleFlag = true;
_settleType = static_cast<int>(settleParams->get(0,0)+0.1);
_t_preComp = settleParams->get(0,1);
_t_postComp = settleParams->get(0,2);
_t_settMax = settleParams->get(0,3);
_settRate_max = settleParams->get(0,4);
_compAcquisitionRate = settleParams->get(0,5); //(Connolly & Baird 2010): a
_compLossRate = settleParams->get(0,6); //(Connolly & Baird 2010): b
_settleOverWhichReefType = settleParams->get(0,7);
//Get settlement probability in time interval dt
_compAcquisitionProb = 1 - exp( -_compAcquisitionRate * (_dt/(24.0*3600)) );
//Convert competency acquisition & loss rates to (dt-1)
_compAcquisitionRate = (_dt/(24.0*3600.0)) * _compAcquisitionRate;
_compLossRate = (_dt/(24.0*3600.0)) * _compLossRate;
std::cout << "\nTEMP COMPARISON: _compAcquisitionRate: " <<_compAcquisitionRate<< ", _compAcquisitionProb: "<<_compAcquisitionProb<<'\n';
_reefs = reefs;
_nbOfShallowReefs = nbOfShallowReefs;
}
/** 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
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment