Commit 32d4e32f authored by thomas's avatar thomas
Browse files

dg: dgParticleTracker clean-up & GBR update

git-svn-id: https://geuz.org/svn/gmsh/trunk@15817 59c8cabf-10d6-4a28-b43f-3936b85956c4
parent 9b4ca250
......@@ -32,7 +32,8 @@ glib="Libraries/lib_gbr.so"
cpuExt=""
if ( len(sys.argv)>1 ):
if sys.argv[1] == "-c":
cpuExt = "/"+sys.argv[2]
cpuExt = "/CPU"+sys.argv[2]
time.sleep(float(sys.argv[2])*2) #Avoid simultaneous simulations crashing each other
# Define simOutput directory & particle concentration
......@@ -45,16 +46,20 @@ LP_OutputDir = '/scratch/cthomas/temp'#simOutputDir +'/LPT_Output/SettProg_Mort0
##
concentration = 2.0 #No. of particles to seed per kmsq
#reefMapShallowFilename = "SGBR_ReefMap_2121_Solid9200.dat"
#reefMapShallowFilenameTrunk = "SGBR_ReefMap_2121_Solid9200.dat"
#reefMapDeepFilenameTrunk = "SGBR_DeepReefMap_2121_Solid9200.dat"
reefMapShallowFilenameTrunk = "CGBR_ReefMap_1093_Cropped_Sma_Solid9200"
reefMapDeepFilenameTrunk = "CGBR_DeepReefMap_1093_Cropped_Sma_Solid9200"
reefMapShallowFilename = reefMapShallowFilenameTrunk+".dat"
reefMapDeepFilename = reefMapDeepFilenameTrunk+".dat"
#Radius: 6.37101e6 for old meshes (sma_zoomed), else 6.371e6
# Radius: 6.37101e6 for old meshes (sma_zoomed), else 6.371e6 (used for reefsObj and reefsDistanceObj)
EarthRadius=6.37101e6
stereoFlag = 0
# Set whether 1: to use PROJ projeciton library (set parameters in GBR.cc) or 0: to use old projection
usePROJ = 1
# If using PROJ, set projection systems:
lonlat_projection = "+proj=latlong +ellps=WGS84"
mesh_projection = "+proj=utm +ellps=WGS84 +zone=55"
##
# SELECT REEF TYPES TO SEED AND SETTLE OVER
......@@ -72,8 +77,6 @@ if(not os.path.exists(LP_OutputDir)):
except: 0;
XYZ = function.getCoordinates()
lonLatDegrees = functionC(glib,"lonLatDegrees",3,[XYZ])
stereoToLonLatDegrees = functionC(glib,"stereoToLonLatDegrees",3,[XYZ])
print(colored('Loading simulation parameters ...',"red"))
simInfoFile = open(simOutputDir+'/simulationInfo.txt','r')
......@@ -120,14 +123,14 @@ mortParams = fullMatrixDouble(4,4)
#[0]: 0: mortality off, 1: const. mortality, 2: Weibull mortality
mortParams.set(0,0, 1.0)
#[1]: (max if Weibull) mortality rate (PER DAY)
#(Connolly 2k11): 0.09 d-1 for G. retiformis, 0.067 for A. Gemmifera
#G. retiformis: 0.09, A. gemmifera: 0.067
# Reef fish: 0.18 (James et al. 2002)
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
#A. millepora: 0.043, A. valida: 0.019, P. daedalea: 0.060
mortParams.set(0,2, 0.043)
#[3]: Weibull only (else unused): nu parameter <-- (Connolly & Baird 2010)nu = c(Pinder et al. 1978)
#A. Millepora: 0.57
#A. millepora: 0.57, A. valida: 0.46, P. daedalea: 0.72
mortParams.set(0,3, 0.57)
# SETTLING PARAMETERS
......@@ -137,7 +140,7 @@ settleParams = fullMatrixDouble(8,8)
#[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) == for type-3: t_timeLag
# G. Retiformis: 0.0, A. Gemmifera: 3.47
# G. retiformis: 0.0, A. gemmifera: 3.47, A. millepora: 3.239, A. valida: 0.0, P. daedalea: 2.937
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)
......@@ -146,10 +149,10 @@ 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)
# G. Retiformis: 0.58, A. Gemmifera: 0.39
# G. retiformis: 0.58, A. gemmifera: 0.39, A. millepora: 0.18, A. valida: 0.22, P. daedalea: 0.39
settleParams.set(0,5, 0.39)
#[6]: (for type-3 only) rate of loss of competence (d-1)
# G. Retiformis: 0.096, A. Gemmifera: 0.145
# G. retiformis: 0.096, A. gemmifera: 0.145, A. millepora: 0.05, A. valida: 0.031, P. daedalea: 0.099
settleParams.set(0,6, 0.145)
#[7]: settle over shallow, deep or both reef types
settleParams.set(0,7, settleOverWhichReefs)
......@@ -165,16 +168,14 @@ swimParams.set(0,2, 2.0)
#[3]: swimming starts after .. hours
swimParams.set(0,3, 0.0)
#STEREO: still need to convert swimming speeds to stereo
print(colored('Loading model ...',"red"))
model = GModel()
if stereoFlag == 0:
print " (Cartesian co-ordinates)"
if usePROJ == 0:
print " (Using old projection)"
model.load("./Meshes/"+simMesh+"_tan.msh")
else:
print " (Stereographic co-ordinates)"
model.load("./Meshes/"+simMesh+".msh")
print " (Using PROJ for projections)"
model.load("./Meshes/"+simMesh+"_utm.msh")
groups = dgGroupCollection(model, 2, 1) #model, dimension, order
#groups.buildGroupsOfInterfaces()
XYZ = function.getCoordinates()
......@@ -217,7 +218,7 @@ extern "C" {
}
"""
tmpCLib = "tempCLib.dylib"
tmpCLib = "./tempCLib.dylib"
if (Msg.GetCommRank() == 0 ) :
functionC.buildLibrary (CCode, tmpCLib)
......@@ -241,13 +242,16 @@ reefsShallowDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapShallowFilenameTr
reefsDeepDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapDeepFilenameTrunk+"_"+simMesh
#1. Load reef maps from data if they don't already exist
UTMProjector = None
if usePROJ == 1:
UTMProjector = slimProjector(mesh_projection, lonlat_projection)
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)
if stereoFlag == 0:
reefsObj = reefs("./Data/"+reefMapShallowFilename, lonLatDegrees, 1, EarthRadius)
#else:
#reefsObj = reefs("./Data/"+reefMapShallowFilename, stereoToLonLatDegrees, 1, EarthRadius)
if usePROJ == 0:
reefsObj = reefsRaw_oldProj("./Data/"+reefMapShallowFilename, EarthRadius)
else:
reefsObj = reefsRaw("./Data/"+reefMapShallowFilename, UTMProjector)
# 2. Interpolate the function onto DC
reefsShallowDC.interpolate(reefsObj)
reefsShallowDC.exportMsh(reefsShallowDCfilename,0,0)
......@@ -260,10 +264,10 @@ else:
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)
if stereoFlag == 0:
reefsObj = reefs("./Data/"+reefMapDeepFilename, lonLatDegrees, 1, EarthRadius)
#else:
#reefsObj = reefs("./Data/"+reefMapDeepFilename, stereoToLonLatDegrees, 1, EarthRadius)
if usePROJ == 0:
reefsObj = reefsRaw_oldProj("./Data/"+reefMapDeepFilename, EarthRadius)
else:
reefsObj = reefsRaw("./Data/"+reefMapDeepFilename, UTMProjector)
# 2. Interpolate the function onto DC
reefsDeepDC.interpolate(reefsObj)
reefsDeepDC.exportMsh(reefsDeepDCfilename,0,0)
......
......@@ -376,6 +376,15 @@ void cutOffMin(dataCacheMap *, fullMatrix<double> &sol, fullMatrix<double> &ourF
//Physical functions:
/** Calculate water depth */
void getWaterDepth(dataCacheMap *, fullMatrix<double> &fct, fullMatrix<double> &sol, fullMatrix<double> &bath) {
for (size_t i = 0; i< fct.size1(); i++) {
double eta = sol.get(i, 0);
double depth = bath.get(i, 0);
fct.set(i, 0, eta+depth);
}
}
/** Constant wind */
void wind (dataCacheMap *, fullMatrix<double> &sol) {
for (size_t i = 0; i< sol.size1(); i++) {
......
......@@ -40,7 +40,7 @@ if(Msg.GetCommRank() == 0):
print(colored('SLIM: Multirate GBR simulation starting. Loading 1:Mesh, 2:Bathymetry 3:Reef Map 4:Conservation law.',"yellow"))
print''
print(colored(' ________ ___ ___ ___ ___ ',"red"))
print(colored(' ________ ___ ___ ___ ___ ',"red"))
print(colored(' / _____// / / / / \/ \ ',"red"))
print(colored(' ____\____ // /___/ /_/__/\__/\ \_____ ',"red"))
print(colored('/__________/ \_____/__/___________\__\____\ [ Great Barrier Reef ]',"red"))
......@@ -382,7 +382,7 @@ elif gbrResidual == 3:
boundaryUV_modulated_CCH = functionC(glib,"currentsAlongVecDotNormals", 3, [function.getNormals(), data_South, boundaryUV_CCH])
boundaryUVElev_CCH = functionC(glib, "merge", 3, [boundaryElev_CCH, boundaryUV_modulated_CCH])
# .... Boundary currents as funciton of wind
# .... Boundary currents as function of wind
elif gbrResidual == 4:
if gbrWind == 0:
print 'Error: gbrResidual type is',gbrResidual,'but gbrWind is',gbrWind,'. Cannot force currents using wind if no wind exists!'
......@@ -570,7 +570,13 @@ if(Msg.GetCommRank() == 0):
exporterSol=dgIdxExporter(solution, outputDir+'/mr-gbr')
#exporterSol=dgIdxExporter(currentFunction, outputDir+'/mr-gbr-func')
norm=solution.norm()
norm = solution.norm()
waterDepth = functionC(glib, "getWaterDepth", 1, [solution.getFunction(), bathDC.getFunction()])
integrator = dgFunctionIntegrator(groups, waterDepth, solution)
totalMass=fullMatrixDouble(1, 1)
totalMassRef=fullMatrixDouble(1, 1)
integrator.compute(totalMassRef, "")
if (Msg.GetCommRank() == 0):
print ''
......@@ -585,6 +591,7 @@ if (Msg.GetCommRank() == 0):
print ' - Number of Iterations:', nbSteps
print ' - Theoretical Speedup:', rk.speedUp()
print ' - Norm of Solution at Ti:', norm
print ' - Mass of Solution at Ti:', totalMassRef(0, 0)
print ''
print ' - Tides:', gbrTideLabels[gbrTide], ' || Wind:', gbrWindLabels[gbrWind], ' || Residual current:', gbrResidualLabels[gbrResidual]
print ''
......@@ -599,8 +606,7 @@ if (Msg.GetCommRank() == 0):
simInfoFile.writelines(simInfo)
simInfoFile.close()
#Function Evaluators
#Initialise Function Evaluators
FunctionEval = dgFunctionEvaluator(groups, currentFunction)
BathEval = dgFunctionEvaluator(groups, bathDC.getFunction())
ElevationEval = dgFunctionEvaluator(groups, solution.getFunction())
......@@ -623,17 +629,23 @@ for point in range(0,len(evalPointsGBROOSLatLon)):
outBathGBROOSEval.write(str(res(0,0))+'\n')
outBathGBROOSEval.close()
#Create total mass file
if (Msg.GetCommRank() == 0):
outMass = open(outputDir+'/totalMass.dat',"w")
outMass.write(str(t/3600-Ti/3600)+'\t'+str(totalMassRef(0,0))+'\n')
outMass.close()
# **********************************************
exporterSol.exportIdx(0, t)
# ****************** Iterate ******************
t_exportStart = Ti + t_exportStart
startcpu = time.clock()
#TEST: OUTPUT VALUE OF WIND (PART 1 OF 2)
if (Msg.GetCommRank() == 0):
for siteId in range(0,len(evalPointsGBROOSLatLon)):
outWind = open(outputDir+'/Evaluator/WINDUV_%02d'%(evalPointsGBROOSLatLon[siteId][3])+'.dat',"w")
outWind.close()
##TEST: OUTPUT VALUE OF WIND (PART 1 OF 2)
#if (Msg.GetCommRank() == 0):
#for siteId in range(0,len(evalPointsGBROOSLatLon)):
#outWind = open(outputDir+'/Evaluator/WINDUV_%02d'%(evalPointsGBROOSLatLon[siteId][3])+'.dat',"w")
#outWind.close()
# OPTIONAL 3/3: if starting from an initial solution, adjust start & end steps accordingly
firstStep = 1
......@@ -669,8 +681,6 @@ for i in range(firstStep, lastStep):
print(colored('------------------------------------------------------', "red"))
print''
if ( i%export == 0 and t > t_exportStart ):
#solution.exportFunctionMsh(currentFunction, outputDir+'/Current/current-%06d'%(i), t, nbExport, '')
#solution.exportMsh(outputDir+"/Solution/sol-%06d"%(i),t,nbExport)
exporterSol.exportIdx(i, t)
if ( i%export_extraTerms == 0 ):
#exporterDissipationTerm.exportIdx(i, t)
......@@ -689,7 +699,12 @@ for i in range(firstStep, lastStep):
exporterGBROOS_CCHSol.exportIdx(i, t)
if ( i%export_fnEval == 0 ):
#Export Function Evaluator results:
nbExport = nbExport + 1
integrator.compute(totalMass, "")
relMass=totalMass(0, 0)
if (Msg.GetCommRank() == 0):
outMass = open(outputDir+'/totalMass.dat',"a")
outMass.write(str(t/3600-Ti/3600)+'\t'+str(relMass)+'\n')
outMass.close()
#*** Current evaluator: ***
for point in range(0,len(evalPointsCart)):
#Only write to file if we are on correct cpu:
......
//Lagrangian random walk advection-diffusion code. Offline particle-tracker. Called from /GreatBarrierReef/LPT_Caller.py
//Contact: C. Thomas
#include <fstream>
#include <iostream>
......@@ -123,42 +122,41 @@ void dgParticleTracker2D::particleMove(dgDofContainer* solution, fullMatrix<doub
//Pre-calculation checks:
//1. Check if particle dies
if (_mortType != 0) {
if (_mortFlag == true) {
if (isDead(t) == 1) {
dead++;
_ParticleArray->getParticle(m)->setState(2);
// outdead << m << "\t" << t << "\t" << x_n << "\t" << y_n << endl;
continue;
}
}
//2. Check if particle competency changes
if (_settleType == 3) {
if (t > _t_preComp) {
//2&3. Check if particle competency changes and if particle settles
if (_settleFlag == true) {
if (_settleType == 3 && t > _t_preComp) {
int particleComp = _ParticleArray->getParticle(m)->getCompetency();
_ParticleArray->getParticle(m)->setCompetency( assessCompetency(particleComp, t) );
}
}
//3. Check if particle settles
//3a. Check if settling is turned on
if (_settleType != 0) {
//3b. Check if we're on a reef
fullMatrix<double> reefFM;
ReefEval.computeElem(iGroup, iElem, 0.25, 0.25, 0.0, reefFM);
int reefId = static_cast<int>(reefFM(0,0));
if ( (reefId > 0 && _settleOverWhichReefType == 2)
|| (reefId > _nbOfShallowReefs && _settleOverWhichReefType == 1)
|| (reefId > 0 && reefId <= _nbOfShallowReefs && _settleOverWhichReefType == 0) ) //We are on REEF!
{
if (isSettled(t,m) == 1) { //NB: could change this to settleMass(double) if needed
settled++;
_ParticleArray->getParticle(m)->setState(1);
//Update connectivity matrix
//NB: subtract 1 from reef IDs so that the lowest ReefId (reefId=1) goes to first row/col of CMatrix (row/col 0)
double previousConnectivity = connectivityMatrix->get(_ParticleArray->getParticle(m)->getSource()-1, reefId-1);
connectivityMatrix->set(_ParticleArray->getParticle(m)->getSource()-1, reefId-1, previousConnectivity+1);
continue;
//3. Check if particle settles
if (_settleType != 0) { //NB: double-check that settling is also enabled in settleParams
//3b. Check if we're on a reef
fullMatrix<double> reefFM;
ReefEval.computeElem(iGroup, iElem, 0.25, 0.25, 0.0, reefFM);
int reefId = static_cast<int>(reefFM(0,0));
if ( (reefId > 0 && _settleOverWhichReefType == 2)
|| (reefId > _nbOfShallowReefs && _settleOverWhichReefType == 1)
|| (reefId > 0 && reefId <= _nbOfShallowReefs && _settleOverWhichReefType == 0) ) //We are on REEF!
{
if (isSettled(t,m) == 1) { //NB: could change this to settleMass(double) if needed
settled++;
_ParticleArray->getParticle(m)->setState(1);
//Update connectivity matrix
//NB: subtract 1 from reef IDs so that the lowest ReefId (reefId=1) goes to first row/col of CMatrix (row/col 0)
double previousConnectivity = connectivityMatrix->get(_ParticleArray->getParticle(m)->getSource()-1, reefId-1);
connectivityMatrix->set(_ParticleArray->getParticle(m)->getSource()-1, reefId-1, previousConnectivity+1);
continue;
}
}
}
}
......@@ -182,24 +180,26 @@ void dgParticleTracker2D::particleMove(dgDofContainer* solution, fullMatrix<doub
alive++;
//5. Calculate swimming velocity
if (_swimType > 0 && t > _swimStart*3600.0) {
//Check if we are close enough to a reef to swim but not yet on a reef
// -> nb: use different criteria for these two tests, as reefsDistFM can=0 even when not on reef - cf. notes
fullMatrix<double> reefsDistFM;
fullMatrix<double> reefFM;
ReefsDistanceEval.computeElem(iGroup, iElem, uvw[0], uvw[1], uvw[2], reefsDistFM);
ReefEval.computeElem(iGroup, iElem, 0.25, 0.25, 0.0, reefFM);
if (reefsDistFM(0,0) < _swimDist*1000.0 && reefFM(0,0) < 0.1) {
//Get vector pointing to nearest reef
fullMatrix<double> reefDir;
ReefsDirectionEval.computeElem(iGroup, iElem, uvw[0], uvw[1], uvw[2], reefDir);
//Normalise vector
double vectorMag = sqrt(reefDir(0,0)*reefDir(0,0) + reefDir(0,1)*reefDir(0,1));
//Calculate swimming vector
pullu = (reefDir(0,0)/vectorMag) * _swimSpeed;
pullv = (reefDir(0,1)/vectorMag) * _swimSpeed;
//Test:
//if (n%1000==0) outSwimming << "VP ("<<x_n<<","<<y_n<<",0.0) {"<<pullu<<","<<pullv<<",0.0};"<<endl;
if (_swimFlag == true) {
if (t > _swimStart*3600.0) {
//Check if we are close enough to a reef to swim but not yet on a reef
// -> nb: use different criteria for these two tests, as reefsDistFM can=0 even when not on reef - cf. notes
fullMatrix<double> reefsDistFM;
fullMatrix<double> reefFM;
ReefsDistanceEval.computeElem(iGroup, iElem, uvw[0], uvw[1], uvw[2], reefsDistFM);
ReefEval.computeElem(iGroup, iElem, 0.25, 0.25, 0.0, reefFM);
if (reefsDistFM(0,0) < _swimDist*1000.0 && reefFM(0,0) < 0.1) {
//Get vector pointing to nearest reef
fullMatrix<double> reefDir;
ReefsDirectionEval.computeElem(iGroup, iElem, uvw[0], uvw[1], uvw[2], reefDir);
//Normalise vector
double vectorMag = sqrt(reefDir(0,0)*reefDir(0,0) + reefDir(0,1)*reefDir(0,1));
//Calculate swimming vector
pullu = (reefDir(0,0)/vectorMag) * _swimSpeed;
pullv = (reefDir(0,1)/vectorMag) * _swimSpeed;
//Test:
//if (n%1000==0) outSwimming << "VP ("<<x_n<<","<<y_n<<",0.0) {"<<pullu<<","<<pullv<<",0.0};"<<endl;
}
}
}
......@@ -307,7 +307,7 @@ void dgParticleTracker2D::particleMove(dgDofContainer* solution, fullMatrix<doub
uvw[0] = uvw[0] - (duvw[0]*remaining);
uvw[1] = uvw[1] - (duvw[1]*remaining);
/*Update remaining vector*/
dxyz[0] *= remaining; //Used to recalculate duvw at beginning of every cycle
dxyz[0] *= remaining; //Used to recalculate duvw at beginning of every cycle
dxyz[1] *= remaining;
/*Retrieve and record current co-ordinates*/
e->primaryPnt(uvw[0],uvw[1],uvw[2], xyz);
......@@ -322,7 +322,7 @@ void dgParticleTracker2D::particleMove(dgDofContainer* solution, fullMatrix<doub
int faceGroupId, faceId, connectionId;
_neighbourMap->getFace(iGroup, iElem, edge_id, faceGroupId, faceId, connectionId);
std::string tag = _group->getFaceGroup(faceGroupId)->physicalTag();
if (tag == "Open Sea North" || tag == "Open Sea Central" || tag == "Open Sea South" || tag == "OpenSeaNorth" || tag == "OpenSeaCentral" || tag == "OpenSeaSouth") {
if (tag == "Open Sea North" || tag == "Open Sea Central" || tag == "Open Sea South" || tag == "Open Sea Neutral") {
_ParticleArray->getParticle(m)->setState(3);
break;
}
......@@ -407,27 +407,25 @@ void dgParticleTracker2D::particleMove(dgDofContainer* solution, fullMatrix<doub
_ParticleArray->getParticle(m)->setPosition(x_nplus1, y_nplus1);
//If t > t_precomp AND settling is disabled, then update connectivity matrix using exposure time method
if (t > _t_preComp && t < _t_postComp && _settleType == 0) {
//Are we on REEF?
fullMatrix<double> reefFM;
ReefEval.computeElem(iGroup, iElem, 0.25, 0.25, 0.0, reefFM);
int reefId = static_cast<int>(reefFM(0,0));
// //Test:
// if (_ParticleArray->getParticle(m)->getSource() < 8) printf("Particle %d, from source reef %d is currently on reef %d\n",m, _ParticleArray->getParticle(m)->getSource(), reefId);
if ( reefId > 0 ) //We are on REEF!
{
//NB: subtract 1 from reef IDs so that the lowest ReefId (reefId=1) goes to first row/col of CMatrix (row/col 0)
double previousConnectivity = connectivityMatrix->get(_ParticleArray->getParticle(m)->getSource()-1, reefId-1);
connectivityMatrix->set(_ParticleArray->getParticle(m)->getSource()-1, reefId-1, previousConnectivity+_dt);
if (_connectivityExposureTimeFlag == true && _settleFlag == false) {
if (t > _t_preComp && t < _t_postComp) {
//Are we on REEF?
fullMatrix<double> reefFM;
ReefEval.computeElem(iGroup, iElem, 0.25, 0.25, 0.0, reefFM);
int reefId = static_cast<int>(reefFM(0,0));
if ( reefId > 0 ) //We are on REEF!
{
//NB: subtract 1 from reef IDs so that the lowest ReefId (reefId=1) goes to first row/col of CMatrix (row/col 0)
double previousConnectivity = connectivityMatrix->get(_ParticleArray->getParticle(m)->getSource()-1, reefId-1);
connectivityMatrix->set(_ParticleArray->getParticle(m)->getSource()-1, reefId-1, previousConnectivity+_dt);
}
}
}
//Reset Reef-proximity advection terms
pullu = 0;
pullv = 0;
}//END: m-loop
//Record no. of particles in each state at the end of this time-steps
......@@ -435,16 +433,17 @@ void dgParticleTracker2D::particleMove(dgDofContainer* solution, fullMatrix<doub
// 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[5] = {0};
//_ParticleArray->getStats(particleStats);
for (int m=0; m<M; m++) {
int state = _ParticleArray->getParticle(m)->getState();
particleStats[state]++;
if (state == 0 && _settleType == 3) {
int aliveAndCompetent = _ParticleArray->getParticle(m)->getCompetency();
if (aliveAndCompetent == 1) particleStats[4]++;
if (_settleFlag == true) {
for (int m=0; m<M; m++) {
int state = _ParticleArray->getParticle(m)->getState();
particleStats[state]++;
if (state == 0 && _settleType == 3) {
int aliveAndCompetent = _ParticleArray->getParticle(m)->getCompetency();
if (aliveAndCompetent == 1) particleStats[4]++;
}
}
}
_ParticleArray->setStats(particleStats);
outStats << t << " " << particleStats[0] << " " << particleStats[1] << " " << particleStats[2] << " " << particleStats[3] << " " << particleStats[4] << endl;
......@@ -715,9 +714,26 @@ void velocity3D (double stereoCoords[2], double stereoVector[2], double xyzVecto
xyzVector[2] = uz;
}
/** Method to print contents of a connectivity matrix */
void connectivityMatrixPrint::CMPrinter(fullMatrix<double>* FM, const std::string filename) {
std::ofstream outFile;
std::ostringstream fname;
//Define filename
fname << filename << ".dat";
//Open file (replace if already exists)
std::string temp = fname.str();
outFile.open( temp.c_str(), std::ios::out | std::ios::trunc );
for (int row=0; row<FM->size1(); row++) {
for (int col=0; col<FM->size2(); col++) {
outFile << FM->get(row,col) << " ";
}
outFile << std::endl;
}
outFile.close();
}
/**Seed particles using reef map DC*/
particleArray::particleArray(dgGroupCollection* groups, dgDofContainer* reefsDC, fullMatrix<double>* initCM, fullMatrix<double>* reefsToSeed, double concentration, int seedOverWhichReefs, int nbOfShallowReefs)
{
particleArray::particleArray(dgGroupCollection* groups, dgDofContainer* reefsDC, fullMatrix<double>* initCM, fullMatrix<double>* reefsToSeed, double concentration, int seedOverWhichReefs, int nbOfShallowReefs) {
/*Plan: loop through every element, evaluate reefsDC in each and seed particles where it is non-zero.
* Number of particles to seed depends on element's surface area
* Particles to be seeded randomly in element
......@@ -838,6 +854,114 @@ particleArray::particleArray(dgGroupCollection* groups, dgDofContainer* reefsDC,
}
void particleArray::printPositions(const std::string filename, const std::string filetype, int status) {
std::ofstream outFile;
std::ostringstream fname, res;
//Define filename
fname << "-";//<< std::setw( 5 ) << std::setfill( '0' ) << m; //<--Creates seperate file for every particle (if want to track particles)
res << filename << fname.str() << "." << filetype;
//Open file (replace if already exists)
std::string temp = res.str();
outFile.open( temp.c_str(), std::ios::out | std::ios::trunc );
//outFile.open( temp.c_str(), std::ios::out | std::ios::app );
if (filetype == "dat") {
//Write to file as list
for (size_t m=0; m<_particle_array.size(); m++) {
//Only register if particle is of state 'status'
if (_particle_array[m].getState() == status) {
outFile << _particle_array[m].x() << " " << _particle_array[m].y() << " " << _particle_array[m].getSource() << std::endl;
}
}
outFile.close();
}
else if (filetype == "pos") {
//Write to file in .pos format
outFile << "View \"particlePositions\" {"<<std::endl;
for (size_t m=0; m<_particle_array.size(); m++) {
//Only register if particle is of state 'status'
if (_particle_array[m].getState() == status) {
outFile << "SP (" << _particle_array[m].x() <<","<< _particle_array[m].y() <<",0) {"<<_particle_array[m].getSource()<<"};" << std::endl;
}
}
outFile << "};"<<std::endl;
outFile.close();
}
else {
std::cout << "ERROR: printPositions filetype argument is invalid." << std::endl;
outFile.close();
}
}
void particleArray::printPositionsFollow(const std::string filename, const std::string filetype, int n, int N, double dt) {
for (size_t m=0; m<_particle_array.size(); m++) {
std::ofstream outFile;
std::ostringstream fname, res;
//Define filename
fname << "-" << std::setw( 5 ) << std::setfill( '0' ) << m; //<--Creates seperate file for every particle
res << filename << fname.str() << "." << filetype;
std::string temp = res.str();
if (n==0) {
outFile.open( temp.c_str(), std::ios::out | std::ios::trunc );
if (filetype == "pos") outFile << "View \"particlePositions\" {"<<std::endl;
}
else outFile.open( temp.c_str(), std::ios::out | std::ios::app );
if (filetype == "dat") {
//Write to file as list
outFile << _particle_array[m].x() << " " << _particle_array[m].y() << std::endl;
outFile.close();
}
else if (filetype == "pos") {
//Write to file in .pos format
outFile << "SP (" << _particle_array[m].x() <<","<< _particle_array[m].y() <<",0) {" << static_cast<double>(n)*dt/3600.0 << "};" << std::endl;
if (n==N) outFile << "};"<<std::endl;
outFile.close();
}
else {
std::cout << "ERROR: printPositions filetype argument is invalid." << std::endl;
outFile.close();
}
}
}
void particleArray::printPositionsFollowSingle(const std::string filename, const std::string filetype, int n, int N, double dt, int particleId) {
std::ofstream outFile;
std::ostringstream fname, res;
//Define filename
fname << "-" << std::setw( 5 ) << std::setfill( '0' ) << particleId; //<--Creates file labelled by particleId
res << filename << fname.str() << "." << filetype;
std::string temp = res.str();
if (n==0) {
outFile.open( temp.c_str(), std::ios::out | std::ios::trunc );
if (filetype == "pos") { outFile << "View \"particlePositions\" {"<<std::endl; }
}