Commit c6c7d5f5 authored by thomas's avatar thomas
Browse files

dg: LPT bug fix

git-svn-id: https://geuz.org/svn/gmsh/trunk@13585 59c8cabf-10d6-4a28-b43f-3936b85956c4
parent 6e2768f0
# 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
# Mesh must be in planar coords in Meshes/ folder
# Mesh must be in planar coords in Meshes/ folder. Stereo-coords: WIP
# 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
#
# !!! NB: MUST UPDATE t_preComp in settleParams, EVEN IF settleType=0 !!!
# !!! NB: Must update t_preComp in settleParams, EVEN IF settleType=0 !!!
# --> Because _t_preComp is defined here, even for exposure time-connectivity method
from dgpy import *
......@@ -48,10 +48,13 @@ LP_OutputDir = 'temp'#scratchFolder + simOutputDirName +'/LPT_Output/NoSett_NoMo
simOutputDir = scratchFolder + simOutputDirName
##
concentration = 2.0 #No. of particles to seed per kmsq
reefMapShallowFilename = "SGBR_ReefMap_2121_Solid9200.dat"
reefMapDeepFilename = "SGBR_DeepReefMap_2121_Solid9200.dat"
#reefMapShallowFilename = "SGBR_ReefMap_2121_Solid9200.dat"
reefMapShallowFilenameTrunk = "CGBR_ReefMap_1093_Cropped_Sma_Solid9200"
reefMapDeepFilenameTrunk = "CGBR_DeepReefMap_1093_Cropped_Sma_Solid9200"
reefMapShallowFilename = reefMapShallowFilenameTrunk+".dat"
reefMapDeepFilename = reefMapDeepFilenameTrunk+".dat"
stereoFlag = 1
stereoFlag = 0
if stereoFlag == 1:
EarthRadius=6371e3
else:
......@@ -61,8 +64,8 @@ else:
# SELECT REEF TYPES TO SEED AND SETTLE OVER
#0: shallow reefs, 1: deep reefs, 2: both
seedOverWhichReefs = 0
settleOverWhichReefs = 0
seedOverWhichReefs = 1
settleOverWhichReefs = 2
##
if(not os.path.exists(LP_OutputDir)):
......@@ -85,7 +88,7 @@ simTotIter = int(simInfo[1])
simExportGap = int(simInfo[2])
simStartTime = int(simInfo[3])
#simOutputDir = simInfo[4].strip()
simMesh = "gbr20K_S_sgbr"#"gbr10K"# #simInfo[5].strip()
simMesh = simInfo[5].strip()
simInfoFile.close()
# Lauch Makefile to generate meshes and download forcings
......@@ -104,14 +107,14 @@ simSolutions = [simExportGap*solutionNumber for solutionNumber in range(0,simTot
print(colored('Loading particle-tracking parameters ...',"red"))
LP_dt = 120. #(double) in seconds
LP_TotTime = 3*7*24.0 #(double) in hours
LP_Export = 1#int( (12.0*3600.0)/LP_dt ) #(int) export particle positions every ... steps
LP_Export = int( (12.0*3600.0)/LP_dt ) #(int) export particle positions every ... steps
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)
LP_StartTime = simStartTime # + (7.0*24.0*3600.0) #in seconds
LP_StartTime = simStartTime + (7.0*24.0*3600.0) #in seconds
LP_EndTime = LP_StartTime + LP_TotTime*3600.0
simItersPerLPIter = LP_dt/simDt
# SEEDING PARAM
# REEF SEEDING PARAM
# -> 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)
......@@ -120,11 +123,11 @@ reefsToSeed.set(0,0,0)
# MORTALITY PARAMETERS
mortParams = fullMatrixDouble(4,4)
#[0]: 0: mortality off, 1: const. mortality, 2: Weibull mortality
mortParams.set(0,0, 0.0)
mortParams.set(0,0, 1.0)
#[1]: (max if Weibull) mortality rate (PER DAY)
#NB: 0.09 d-1 for G. retiformis (Connolly 2k11), 0.07 for A. Gemmifera
#(Connolly 2k11): 0.09 d-1 for G. retiformis, 0.067 for A. Gemmifera
# Reef fish: 0.18 (James et al. 2002)
mortParams.set(0,1, 0.09)
mortParams.set(0,1, 0.067)
#[2]: Weibull only (else unused): lambda parameter <-- (Connolly & Baird 2010)lambda = 1/b(Pinder et al. 1978)
#A. Millepora: 0.043
mortParams.set(0,2, 0.043)
......@@ -138,19 +141,21 @@ 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)
#[1]: t_preComp (seconds) = t_timeLag (for [3])
settleParams.set(0,1, 0.0)
#settleParams.set(0,1, 5.0*24.0*3600.0)
#[2]: t_postComp
settleParams.set(0,2, 25.0*24.0*3600.0)
#[1]: t_preComp (seconds) == for type-3: t_timeLag
# A. Gemmifera: 3.47
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)
#[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)
settleParams.set(0,5, 0.5)
# A. Gemmifera: 0.39
settleParams.set(0,5, 0.39)
#[6]: (for type-3 only) rate of loss of competence (d-1)
settleParams.set(0,6, 0.0)
# A. Gemmifera: 0.145
settleParams.set(0,6, 0.145)
#[7]: settle over shallow, deep or both reef types
settleParams.set(0,7, settleOverWhichReefs)
......@@ -165,7 +170,7 @@ swimParams.set(0,2, 2.0)
#[3]: swimming starts after .. hours
swimParams.set(0,3, 0.0)
#TODO: convert swimming speeds to stereo
#STEREO: still need to convert swimming speeds to stereo
print(colored('Loading model ...',"red"))
model = GModel()
......@@ -234,13 +239,13 @@ diffProjector = L2ProjectionContinuous(dofMan)
diffProjector.apply(diffDC,diffFC)
diffDC.exportMsh("./LPT_Functions/DIFFUSIVITY_smooooth_"+simMesh,0,0)
# Load reef map(s).
# Load reef maps
reefsShallowDC = dgDofContainer(groups,1)
reefsDeepDC = dgDofContainer(groups,1)
reefsShallowDCfilename = "LPT_Functions/REEFS_RAW_LPT_SHALLOW_"+simMesh
reefsDeepDCfilename = "LPT_Functions/REEFS_RAW_LPT_DEEP_"+simMesh
reefsShallowDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapShallowFilenameTrunk+"_"+simMesh
reefsDeepDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapDeepFilenameTrunk+"_"+simMesh
#1. Load data if it doesn't already exist
#1. Load reef maps from data if they don't already exist
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)
......@@ -273,8 +278,7 @@ else:
print(colored('Loading deep reef map from file ...',"red"))
reefsDeepDC.importMsh(reefsDeepDCfilename)
#2. Load data
#2. Load reefs data
minReefNbFM = fullMatrixDouble()
maxReefNbFM = fullMatrixDouble()
reefsShallowDC.minmax(minReefNbFM,maxReefNbFM)
......@@ -366,7 +370,7 @@ else:
print ''
print '-----------------------'
print 'SIMULATION PARAMETERS:'
print 'HYDRODYNAMICS PARAMETERS:'
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)'
......@@ -379,18 +383,22 @@ 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 '-----------------------'
print ''
print '-----------------------'
print 'BIOLOGICAL PARAMETERS:'
print 'Mortality flag:', mortParams(0,0), ' (0: OFF, 1: CONST, 2: WEIBULL)'
if (int(mortParams(0,0)) == 1):
print 'Mortality rate (max):', mortParams(0,1), '(per day)'
print 'Mortality rate (const):', mortParams(0,1), '(per day)'
elif (int(mortParams(0,0)) == 2):
print 'Weibull lambda:', mortParams(0,2), ', Weibull nu:', mortParams(0,3)
print '-----------------------'
print 'Settling flag: ', settleParams(0,0), ' (0: OFF, 1: CONST, 2: TRIANGLE)'
print 'Settling flag: ', settleParams(0,0), ' (0: OFF, 1: CONST, 2: TRIANGLE, 3: PROGRESSIVE)'
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)'
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)'
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 '-----------------------'
......@@ -400,7 +408,7 @@ if (swimParams(0,0) > 0):
print 'Swimming max distance: ',swimParams(0,2), 'km'
print 'Swimming starts after: ', swimParams(0,3), 'hrs'
print '-----------------------'
print ''
LPsim_TimeOffset = LP_StartTime - simStartTime
print '-----------------------'
print 'Particle-tracker / Simulation time offset is:', LPsim_TimeOffset/(24.0*3600.0),'(days)'
......@@ -423,7 +431,7 @@ particleTracker = dgParticleTracker2D(groups, particles, bathDC, diffDC, reefsDu
startcpu = time.clock()
#Keep track of the solutions loaded into the three DCs:
loadedSolutions = [-1] * 3 #[0]:simSolution, [1]:simSolutionOver, [2]:simSolutionUnder
print 'Beginning loop ...'
for n in range(1,LP_TotIter+1):
t = simStartTime + LPsim_TimeOffset + float(n)*LP_dt
iterNumberWanted = int(LPsim_TimeOffset/simDt + n*simItersPerLPIter)
......@@ -516,9 +524,9 @@ for n in range(1,LP_TotIter+1):
# Print output:
if ( n%(LP_Export) == 0 ):
particles.printPositions(LP_OutputDir+'/particlesAlive_%06d'%( n ), 'pos', 0)
if (int(settleParams(0,0)) == 1):
if (int(settleParams(0,0)) != 0):
particles.printPositions(LP_OutputDir+'/particlesSettled_%06d'%( n ), 'pos', 1)
if (int(mortParams(0,0)) == 1):
if (int(mortParams(0,0)) != 0):
particles.printPositions(LP_OutputDir+'/particlesDead_%06d'%( n ), 'pos', 2)
# *** More useful only for SEEDING AT SINGLE POINT ***
# if ( n%LP_Export == 0 ):
......
......@@ -466,13 +466,11 @@ if (startFromExistingSolution == 1):
else:
firstStep = 1
lastStep = nbSteps + 1
print "Beginning loop ..."
for i in range(firstStep, lastStep):
if(i == lastStep-1):
dt=Tf-t
print"got here"
norm = rk.iterate (solution, dt, t)
print"not here"
t = t +dt
timeFunction.set(t)
if ( i%iter == 0 ):
......
......@@ -240,11 +240,6 @@ void dgParticleTracker2D::particleMove(dgDofContainer* solution, fullMatrix<doub
if (u == 0 && v == 0) cout<<"Warning! Advective components are null! n is "<<n<<endl;
//10. Compute next step of path (ie. pos.n at t=t+_dt)
//TODO: convert two components to stereo
dx=0;
dy=0;
randoox[m]=0;
randooy[m]=0;
x_nplus1 = x_n + dx + randoox[m] * sqrt( 2.0 * K_xyplusHalfdxdy * _dt );
y_nplus1 = y_n + dy + randooy[m] * sqrt( 2.0 * K_xyplusHalfdxdy * _dt );
......@@ -434,11 +429,15 @@ void dgParticleTracker2D::particleMove(dgDofContainer* solution, fullMatrix<doub
//Record no. of particles in each state at the end of this time-steps
// -> NB: the above loop 'sees' all particles that are still alive,
// but not all the particles that are dead or settled or out of domain (only the ones settling/dying/leaving domain in this loop)
int particleStats[4] = {0};
int particleStats[5] = {0};
//_ParticleArray->getStats(particleStats);
for (int m=0; m<M; m++) {
int state = _ParticleArray->getParticle(m)->getState();
particleStats[state]++;
if (_settleType != 0) {
int aliveAndCompetent = _ParticleArray->getParticle(m)->getCompetency();
if (aliveAndCompetent == 1) particleStats[4]++;
}
}
// particleStats[0] = alive;
// particleStats[1] += settled;
......@@ -446,8 +445,8 @@ void dgParticleTracker2D::particleMove(dgDofContainer* solution, fullMatrix<doub
_ParticleArray->setStats(particleStats);
outStats << t << " " << particleStats[0] << " " << particleStats[1] << " " << particleStats[2] << " " << particleStats[3] << endl;
printf("End of loop cycle %d of %d. Particle stats: %d alive, %d settled, %d dead, %d out of domain. TOTAL: %d\n", n, _N, particleStats[0], particleStats[1], particleStats[2], particleStats[3], particleStats[0]+particleStats[1]+particleStats[2]+particleStats[3]);
outStats << t << " " << particleStats[0] << " " << particleStats[1] << " " << particleStats[2] << " " << particleStats[3] << " " << particleStats[4] << endl;
printf("End of loop cycle %d of %d. Particle stats: %d alive, %d settled, %d dead, %d out of domain. %d alive & competent. TOTAL: %d\n", n, _N, particleStats[0], particleStats[1], particleStats[2], particleStats[3], particleStats[4], particleStats[0]+particleStats[1]+particleStats[2]+particleStats[3]);
// cout << "Cycle finished." << endl;
// cout << "Alive: " << alive << "; Just dead: " << dead << "; Just settled = " << settled << "." << endl;
// cout << "\t Random numbers: "<< randoraw() <<", "<< randoraw() << endl;
......@@ -969,6 +968,7 @@ particleArray::particleArray(dgGroupCollection* groups, dgDofContainer* reefsDC,
_alive = static_cast<int>( _particle_array.size() );
_settled = 0;
_dead = 0;
_aliveAndCompetent = static_cast<int>( _particle_array.size() );
printf("**** particleArray:: SEEDED %d PARTICLES. ****\n --> ... of which %d (%.2f) over shallow reefs and %d (%.2f) over deep reefs.\n", static_cast<int>( _particle_array.size() ), nbOfShallowParticles, 100.0*static_cast<double>(nbOfShallowParticles)/static_cast<double>( _particle_array.size() ), nbOfDeepParticles, 100.0*static_cast<double>(nbOfDeepParticles)/static_cast<double>( _particle_array.size() ) );
printf(" --> We have %d 'reef' elements, of which %d shallow and %d deep.\n --> We have seeded %d 'reef' elements.\n", nbOfReefElements, nbOfShallowReefElements, nbOfDeepReefElements, nbOfElementsSeeded);
printf(" --> We have seeded %d reefs, or %.2f of total (%.2f of shallow reefs, %.2f of deep reefs).\n", static_cast<int>(reefsSeeded.size()), 100.0*static_cast<double>( reefsSeeded.size() )/static_cast<double>( maxReefId ), 100.0*static_cast<double>(nbOfShallowReefsSeeded)/static_cast<double>(nbOfShallowReefs), 100.0*static_cast<double>(nbOfDeepReefsSeeded)/static_cast<double>(maxReefId-nbOfShallowReefs) );
......
......@@ -152,7 +152,7 @@ public:
class particleArray {
private:
std::vector<particle> _particle_array;
int _alive, _settled, _dead, _outOfDomain;
int _alive, _settled, _dead, _outOfDomain, _aliveAndCompetent;
public:
//Constructors:
......@@ -180,6 +180,7 @@ public:
_settled = 0;
_dead = 0;
_outOfDomain = 0;
_aliveAndCompetent = M;
}
/** Single point constructor (try to avoid using this and use the above instead) */
......@@ -210,14 +211,7 @@ public:
/**Seed particles using reef map DC*/
particleArray(dgGroupCollection* groups, dgDofContainer* reefsDC, fullMatrix<double>* initCM, fullMatrix<double>* reefsToSeed, double concentration, int seedOverWhichReefs, int nbOfShallowReefs);
void setStats(int *stats) {
_alive = stats[0];
_settled = stats[1];
_dead = stats[2];
_outOfDomain = stats[3];
}
void addParticles(dgGroupCollection* _group, int M, double x0, double y0, int status, int source_reef)
{
//Find group & element IDs
......@@ -233,14 +227,24 @@ public:
_particle_array.push_back( particle(x0,y0,iGroup,iElem,source_reef,status) );
}
_alive += M;
_aliveAndCompetent += M;
}
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() {
......
Supports Markdown
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