Commit 8c396dc2 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

rm duplicated files

parent aeb1ad59
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"
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);
void convertVectorCoords(double target_v[3],double original_v[3], double jac[3][3]);
SPoint2 intersection(SPoint2 p1, SPoint2 p2, SPoint2 p3, SPoint2 p4, int& flag);
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;
function* _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:
_gm = _group->getModel();
_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(function* 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
9.1810514739
329375
180
1229558400
/scratch/cthomas/Hydro_Realistic/460K_tOnwCSFRbNCJ0p8SvSth_FWind1o50Every_bDn016Reefs10
gbr460K_640m_nov2012
\ No newline at end of file
l = 5e5;
Point(1) = {-l, -l, 0};
Point(2) = {l, -l, 0};
Point(3) = {l, l, 0};
Point(4) = {-l, l, 0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(5) = {3, 4, 1, 2};
Field[1] = MathEval;
Field[1].F = "1e5";
Background Field = 1;
//Field[1] = MathEval;
//Field[1].F = "0.5e4+1*(abs(x)+abs(y))";
//Background Field = 1;
Plane Surface(6) = {5};
Physical Surface("Surface")={6};
Physical Line("Wall") = {1,2,3,4};
Mesh.CharacteristicLengthExtendFromBoundary=1;
Mesh.CharacteristicLengthFromPoints=1;
Mesh.SecondOrderIncomplete=0;
# Lagrangian particle-tracker
# Test for dgParticleTracker2D
from dgpy import *
import time
import os.path
import subprocess
from dgpy.scripts import Common
from dgpy.scripts.termcolor import colored
import numpy as np
print ('')
print(colored('Lagrangian Advection-Diffusion particle tracking module for SLIM. Test #1: basic', "yellow"))
print ('')
#---------- SET OUTPUT FOLDER --------------------------------------------------
LP_OutputDir = 'output'
if(not os.path.exists(LP_OutputDir)):
try : os.mkdir(LP_OutputDir)
except: 0
#-------------------------------------------------------------------------------
# -------- SET HYDRO SIMULATION PARAMETERS -------------------------------------
print(colored('Loading simulation parameters ...',"red"))
simInfoFile = open('simulationInfo.txt','r') # File with info on hydro simulation
simInfo = simInfoFile.readlines()
simDt = float(simInfo[0]) # Simulation dt
simTotIter = int(simInfo[1]) # Tot. number of iters
simExportGap = int(simInfo[2]) # Gap (in iters) between exports
simStartTime = int(float(simInfo[3])) # Simulation start time
simMesh = simInfo[5].strip() # Mesh (not used here)
simInfoFile.close()
#-------------------------------------------------------------------------------
#--------- SET LPT TIME PARAMETERS ---------------------------------------------
print(colored('Loading particle-tracking parameters ...',"red"))
LP_dt = 90. #(double) in seconds
LP_TotTime = 1.0 #(double) in hours
LP_Export = int( (1.0*(3600.0-90.0))/LP_dt ) #(int) export particle positions every ... steps
LP_StartTime = simStartTime
LP_TotIter = int(LP_TotTime*3600.0/LP_dt)
LP_EndTime = LP_StartTime + LP_TotTime*3600.0
simItersPerLPIter = LP_dt/simDt
#-------------------------------------------------------------------------------
#---------- GENERATE & LOAD MESH & BATHYMETRY ----------------------------------
print(colored('Generating and loading mesh and bathymetry ...',"red"))
model = GModel()
# The mesh "stommel_square.msh" is generated using the geometry file "stommel_square.geo"
# The generated mesh is 2d (second argument) and of order 1 (third argument)
Common.genMesh('stommel_square', 2, 1)
model.load('stommel_square.msh')
# Polynomial order of the discretization, and number of dimensions
order = 3
dimension = 2
# Groups are used to differentiate different elements of the mesh
# Surface elements and edges belong to different groups
# Group can be split into sub-groups (i.e. split groups of elements that belong to searate physical entities)
groups = dgGroupCollection(model, dimension, order)
# Bathymetry
bathFC = functionConstant(20.0)
bathDC = dgDofContainer(groups, 1)
bathDC.interpolate(bathFC)
#-------------------------------------------------------------------------------
#------- SET UP DIFFUSIVITY ----------------------------------------------------
print(colored('Defining diffusivity ...',"red"))
def diff(cmap, val, alpha):
maxEdges = dgConservationLawFunction.maxEdge(cmap)
val[:] = alpha * 0.00020551 * np.power(maxEdges[:, 0], 1.15)
alpha = functionConstant(1.0)
diffFC = functionNumpy(1, diff, [alpha])
diffDC = dgDofContainer(groups,1)
# Smooth diffusivity:
sys = linearSystemPETScDouble()
dofMan = dgDofManager.newCG (groups, 1, sys)
diffProjector = L2ProjectionContinuous(dofMan)
diffProjector.apply(diffDC,diffFC)
#-------------------------------------------------------------------------------
#--------- SEED PARTICLES ------------------------------------------------------
print(colored('Seeding particles ...',"red"))
#Set max width of domain to seed one group of particles close to walls
LengthOfDomain = 5e5
seedPoints = [[0, 0], [-100000,-100000], [-100000,100000], [100000,-100000], [100000,100000], [LengthOfDomain-1, LengthOfDomain-1] ]
particles = particleArray()
for location in range(0,len(seedPoints)):
particles.addParticlesAtPoint(groups, 4000, seedPoints[location][0], seedPoints[location][1], 0, 0, location)
particles.printPositions('seededParticles','pos',0)
particles.printPositions('seededParticles','dat',0)
#-------------------------------------------------------------------------------
#-------- PRINT SIMULATION INFO ------------------------------------------------
print ('')
print ('-----------------------')
print ('HYDRODYNAMICS PARAMETERS:')
print ('Dt:', simDt,'s')
print ('Simulation length:', simTotIter, '(iter) |', simTotIter*simDt/3600.0, '(h)' )