Skip to content
GitLab
Commit 75443e33 authored by thomas's avatar thomas
Browse files

dg: LPT: update benchmark 

git-svn-id: https://geuz.org/svn/gmsh/trunk@16388 59c8cabf-10d6-4a28-b43f-3936b85956c4
parent 18b28c19
Hide whitespace changes
Inline Side-by-side
benchmarks/GreatBarrierReef/LPT_Caller_DualReefMap.py
View file @ 75443e33
......@@ -76,7 +76,7 @@ simInfo = simInfoFile.readlines()
simDt = float(simInfo[0])
simTotIter = int(simInfo[1])
simExportGap = int(simInfo[2])
simStartTime = int(simInfo[3])
simStartTime = int(float(simInfo[3]))
#simOutputDir = simInfo[4].strip()
simMesh = simInfo[5].strip()
simInfoFile.close()
......@@ -110,10 +110,12 @@ simItersPerLPIter = LP_dt/simDt
#--------- SEEDING, MORTALITY, SETTLING AND SWIMMING PARAMETERS ----------------
# SELECT REEF TYPES TO SEED AND SETTLE OVER
#0: shallow reefs, 1: deep reefs, 2: both
seedOverWhichReefs = 2
settleOverWhichReefs = 2
seedOverWhichReefs = 0
settleOverWhichReefs = 0
# SET SEEDING CONCENTRATION (nb particles/kmsq)
concentration = 2.0
# SET MINIMUM NUMBER OF PARTICLES TO SEED ON EACH REEF AT EACH SEEDING EVENT
minNbParticlesPerReef = 200
# Do we seed all in one go (0) or over a time interval (1)
seedOverTime = 0
# If seedOverTime==1, specify the total time interval (seconds) over which to seed, and how many times to seed within this interval
......@@ -191,8 +193,29 @@ XYZ_DC = dgDofContainer(groups, 3)
XYZ_DC.interpolate(XYZ)
XYZ_DC.exportMsh("XYZ_DC",0,0)
bathDC = dgDofContainer(groups, 1);
# Print size of smallest element:
maxEdge = 100000.0
maxMaxEdge = 0.0
minSurface = 1000000.0
maxSurface = 0.0
nbGroups = groups.getNbElementGroups()
for group in range( 0, nbGroups):
groupOfElements = groups.getElementGroup(group)
nbElements = groupOfElements.getNbElements()
for element in range( 0, nbElements):
el = groupOfElements.getElement(element)
if el.getVolume() < minSurface:
minSurface = el.getVolume()
if el.getVolume() > maxSurface:
maxSurface = el.getVolume()
if el.maxEdge() < maxEdge:
maxEdge = el.maxEdge()
if el.maxEdge() > maxMaxEdge:
maxMaxEdge = el.maxEdge()
print 'Smallest element\'s maxEdge is', maxEdge, 'm, Largest element\'s maxEdge is', maxMaxEdge, 'm'
print 'Smallest element\'s surface is', minSurface, 'msq, Largest element\'s surface is', maxSurface, 'msq'
bathDC = dgDofContainer(groups, 1);
if(os.path.exists("./Bath/"+simMesh+"_bath_smooth/"+simMesh+"_bath_smooth.idx")):
print(colored('Importing bathymetry ...',"red"))
bathDC.readMsh("./Bath/"+simMesh+"_bath_smooth/"+simMesh+"_bath_smooth.idx")
......@@ -297,56 +320,38 @@ reefsDeepDC.minmax(minReefNbFM,maxReefNbFM)
nbOfDeepReefs = int(maxReefNbFM(0,0))
nbOfReefs = nbOfShallowReefs + nbOfDeepReefs
print "Maximum possible number of reefs: ",nbOfReefs, "(shallow:", nbOfShallowReefs, ", deep:", nbOfDeepReefs, ")"
#-------------------------------------------------------------------------------
# Record simulation info to file
#------- WRITE SIMULATION INFO TO FILE -----------------------------------------
LP_SimInfo = [str(nbOfDeepReefs),'\n'+str(LP_dt),'\n',str((LP_StartTime-simStartTime)/3600),'\n', str(LP_TotTime),'\n',str(LP_Export),'\n', str(simItersPerLPIter),'\n',str(mortParams(0,0)),'\n' ,str(mortParams(0,1)/(24.0*3600.0)),'\n',str(settleParams(0,0)),'\n','Pre-comptetency period:'+str(settleParams(0,1)/(24.0*3600.0))+' (days)\n', 'Post-competent at:'+str(settleParams(0,2)/(24.0*3600.0))+' days)\n', '(max) settling rate:'+str(settleParams(0,4))+' per dt\nConcentration:'+str(concentration)+'/kmsq\n', 'Seeding over: '+str(seedOverWhichReefs)+'\n', 'Settling over: '+str(settleOverWhichReefs)+'\n', 'Highest shallow reef Id: '+str(nbOfShallowReefs)+'\n']
LP_SimInfoFile = open(LP_OutputDir+'/LP_runInfo.txt','w')
LP_SimInfoFile.writelines(LP_SimInfo)
LP_SimInfoFile.close()
#-------------------------------------------------------------------------------
# Print size of smallest element:
maxEdge = 100000.0
maxMaxEdge = 0.0
minSurface = 1000000.0
maxSurface = 0.0
nbGroups = groups.getNbElementGroups()
for group in range( 0, nbGroups):
groupOfElements = groups.getElementGroup(group)
nbElements = groupOfElements.getNbElements()
for element in range( 0, nbElements):
el = groupOfElements.getElement(element)
if el.getVolume() < minSurface:
minSurface = el.getVolume()
if el.getVolume() > maxSurface:
maxSurface = el.getVolume()
if el.maxEdge() < maxEdge:
maxEdge = el.maxEdge()
if el.maxEdge() > maxMaxEdge:
maxMaxEdge = el.maxEdge()
print 'Smallest element\'s maxEdge is', maxEdge, 'm, Largest element\'s maxEdge is', maxMaxEdge, 'm'
print 'Smallest element\'s surface is', minSurface, 'msq, Largest element\'s surface is', maxSurface, 'msq'
print(colored('Calculating distances to nearest shallow reefs ...', "red"))
distToReefsDC = dgDofContainer(groups, 1)
directionToReefsDC = dgDofContainer(groups, 3)
if(os.path.exists("./LPT_Functions/DISTANCE_TO_REEFS_"+simMesh+"_COMP_0.msh") and os.path.exists("./LPT_Functions/DIRECTION_TO_REEFS_"+simMesh+"_COMP_0.msh")):
print'./LPT_Functions/DISTANCE_TO_REEFS_'+simMesh+'_COMP_0.msh and ./LPT_Functions/DIRECTION_TO_REEFS_'+simMesh+'_COMP_0.msh files already exist: importing them.'
distToReefsDC.importMsh("./LPT_Functions/DISTANCE_TO_REEFS_"+simMesh)
directionToReefsDC.importMsh("./LPT_Functions/DIRECTION_TO_REEFS_"+simMesh)
else:
#TODO: Fix distance_to_reef generation in stereo
print'./LPT_Functions/DISTANCE_TO_REEFS_'+simMesh+'_COMP_0.msh and/or ./LPT_Functions/DIRECTION_TO_REEFS_'+simMesh+'_COMP_0.msh files do not exist. Generating dof container with distance and directions to nearest reefs ...'
reefDistanceObj = reefDistance(groups, reefsShallowDC, XYZ_DC, EarthRadius)
reefDistanceObj.getNearestReefMap(distToReefsDC, directionToReefsDC)
distToReefsDC.exportMsh("./LPT_Functions/DISTANCE_TO_REEFS_"+simMesh,0,0)
directionToReefsDC.exportMsh("./LPT_Functions/DIRECTION_TO_REEFS_"+simMesh,0,0)
#Also save DIR-TO-REEFS as function, for visualisation
directionToReefsDC.exportFunctionMsh(directionToReefsDC.getFunction(), "./LPT_Functions/DIRECTION_TO_REEFS_vector_"+simMesh,0,0)
#dir2RF = directionToReefsDC.getFunction()
#dirFunction = functionC(glib,"current2",3,[directionToReefsDC.getFunction()])
#directionToReefsDC.exportMshNodeData("./DIRECTION_TO_REEFS",0,0, '')
#------- LOAD MAPS OF DISTANCE AND DIRECTION TO NEAREST REEF -------------------
distToReefsDC = None
directionToReefsDC = None
reefsDualDC = None
nbOfShallowReefsFConst = None
if swimParams(0,0) > 0.1:
print(colored('Calculating distances to nearest shallow reefs ...', "red"))
distToReefsDC = dgDofContainer(groups, 1)
directionToReefsDC = dgDofContainer(groups, 3)
if(os.path.exists("./LPT_Functions/DISTANCE_TO_REEFS_"+simMesh+"_COMP_0.msh") and os.path.exists("./LPT_Functions/DIRECTION_TO_REEFS_"+simMesh+"_COMP_0.msh")):
print'./LPT_Functions/DISTANCE_TO_REEFS_'+simMesh+'_COMP_0.msh and ./LPT_Functions/DIRECTION_TO_REEFS_'+simMesh+'_COMP_0.msh files already exist: importing them.'
distToReefsDC.importMsh("./LPT_Functions/DISTANCE_TO_REEFS_"+simMesh)
directionToReefsDC.importMsh("./LPT_Functions/DIRECTION_TO_REEFS_"+simMesh)
else:
#TODO: Fix distance_to_reef generation in stereo
print'./LPT_Functions/DISTANCE_TO_REEFS_'+simMesh+'_COMP_0.msh and/or ./LPT_Functions/DIRECTION_TO_REEFS_'+simMesh+'_COMP_0.msh files do not exist. Generating dof container with distance and directions to nearest reefs ...'
reefDistanceObj = reefDistance(groups, reefsShallowDC, XYZ_DC, EarthRadius)
reefDistanceObj.getNearestReefMap(distToReefsDC, directionToReefsDC)
distToReefsDC.exportMsh("./LPT_Functions/DISTANCE_TO_REEFS_"+simMesh,0,0)
directionToReefsDC.exportMsh("./LPT_Functions/DIRECTION_TO_REEFS_"+simMesh,0,0)
#Also save DIR-TO-REEFS as function, for visualisation
directionToReefsDC.exportFunctionMsh(directionToReefsDC.getFunction(), "./LPT_Functions/DIRECTION_TO_REEFS_vector_"+simMesh,0,0)
nbOfShallowReefsFConst = functionConstant(nbOfShallowReefs)
reefsDualDC = dgDofContainer(groups, 1)
......@@ -369,38 +374,39 @@ if seedOverTime == 1:
print' -> Total concentration:', concentration*seedNbTimes, ', concentration per seeding event:', concentration
# ***SEED OVER ALL REEFS***
#initCM = fullMatrixDouble( nbOfReefs, nbOfReefs )
#particles = particleArray()
#particles.addParticlesWithReefMap(groups, reefsDualDC, initCM, reefsToSeed, concentration, seedOverWhichReefs, nbOfShallowReefs)
#particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
#dummyPrinterICM = connectivityMatrixPrint()
#if settleParams(0,0) == 0:
#dummyPrinterICM.CMPrinter(initCM, LP_OutputDir+'/connectivityMatrix_initial' )
#else:
#dummyPrinterICM.CMPrinter(initCM, LP_OutputDir+'/settlementMatrix_initial' )
initCM = fullMatrixDouble( nbOfReefs, nbOfReefs )
particles = particleArray()
particles.addParticlesWithReefMap(groups, reefsDualDC, initCM, reefsToSeed, concentration, minNbParticlesPerReef, seedOverWhichReefs, nbOfShallowReefs)
#particles.addParticlesWithDofContainer(groups, reefsDualDC, concentration, 200)
particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
dummyPrinterICM = connectivityMatrixPrint()
if settleParams(0,0) == 0:
dummyPrinterICM.CMPrinter(initCM, LP_OutputDir+'/connectivityMatrix_initial' )
else:
dummyPrinterICM.CMPrinter(initCM, LP_OutputDir+'/settlementMatrix_initial' )
# ******
# ***SEED AT A POINT*** (NB: Look at method used to print results if changing seeding method)
UTMProjectorLonLatToUTM = slimProjector(lonlat_projection, mesh_projection)
reefPointsLatLon = [ [-14.7,145.4,0], [-14.7,145.4,1] ]
reefPointsCart = [ [0]*4 for i in range(len(reefPointsLatLon)) ]
pj_latlong = pj_init_plus(lonlat_projection)
pj_merc = pj_init_plus(mesh_projection)
for i in range(0,len(reefPointsLatLon)):
reefPointsCart[i][3] = reefPointsLatLon[i][2]
src_points = fullMatrixDouble(2,2)
dest_points = fullMatrixDouble(2,2)
src_points.set(0,0,reefPointsLatLon[i][1]/180.0*pi)
src_points.set(1,1,reefPointsLatLon[i][0]/180.0*pi)
UTMProjectorLonLatToUTM.projectPoint(src_points, dest_points)
reefPointsCart[i][0] = dest_points(0,0)
reefPointsCart[i][1] = dest_points(1,1)
#UTMProjectorLonLatToUTM = slimProjector(lonlat_projection, mesh_projection)
#reefPointsLatLon = [ [-14.7,145.4,0], [-14.7,145.4,1] ]
#reefPointsCart = [ [0]*4 for i in range(len(reefPointsLatLon)) ]
#pj_latlong = pj_init_plus(lonlat_projection)
#pj_merc = pj_init_plus(mesh_projection)
#for i in range(0,len(reefPointsLatLon)):
#reefPointsCart[i][3] = reefPointsLatLon[i][2]
#src_points = fullMatrixDouble(2,2)
#dest_points = fullMatrixDouble(2,2)
#src_points.set(0,0,reefPointsLatLon[i][1]/180.0*pi)
#src_points.set(1,1,reefPointsLatLon[i][0]/180.0*pi)
#UTMProjectorLonLatToUTM.projectPoint(src_points, dest_points)
#reefPointsCart[i][0] = dest_points(0,0)
#reefPointsCart[i][1] = dest_points(1,1)
particles = particleArray()
for location in range(0,len(reefPointsCart)):
particles.addParticlesAtPoint(groups, 100, reefPointsCart[location][0], reefPointsCart[location][1], 0, reefPointsCart[location][3])
particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
particles.printPositions(LP_OutputDir+'/seededParticles','dat',0)
#particles = particleArray()
#for location in range(0,len(reefPointsCart)):
#particles.addParticlesAtPoint(groups, 100, reefPointsCart[location][0], reefPointsCart[location][1], 0, reefPointsCart[location][3])
#particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
#particles.printPositions(LP_OutputDir+'/seededParticles','dat',0)
# ******
#-------------------------------------------------------------------------------
......@@ -472,7 +478,7 @@ print 'Initialising particle tracker ...'
particleTracker = dgParticleTracker2D(groups, particles, bathDC, diffDC, LP_TotIter, LP_dt, LP_OutputDir)
particleTracker.addMortality(mortParams)
particleTracker.addSettling(settleParams, reefsDualDC, nbOfShallowReefs)
particleTracker.addSwimming(swimParams, reefsDualDC, distToReefsDC, directionToReefsDC)
#particleTracker.addSwimming(swimParams, reefsDualDC, distToReefsDC, directionToReefsDC)
#-------------------------------------------------------------------------------
#-------- PARTICLE LOOP --------------------------------------------------------
......
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