Commit 8a2aa46f authored by thomas's avatar thomas
Browse files

dg: LPT: script update

git-svn-id: https://geuz.org/svn/gmsh/trunk@16383 59c8cabf-10d6-4a28-b43f-3936b85956c4
parent ac321918
......@@ -20,6 +20,7 @@ print ''
print(colored('Lagrangian Advection-Diffusion particle tracking module for SLIM. Working offline.', "yellow"))
print ''
#--------- SETUP LIBRARIES ----------------------------------------------------
if (Msg.GetCommRank()==0):
print ''
print(colored('Compiling Libraries ...', "red"))
......@@ -28,33 +29,40 @@ if (Msg.GetCommRank()==0):
glib="Libraries/lib_gbr.so"
#rlib="Libraries/lib_reef.so"
#-------------------------------------------------------------------------------
#---------- SET INPUT & OUTPUT FOLDERS -----------------------------------------
#simOutputDir = '/scratch/cthomas/Hydro_Realistic/gbr500K_sma_1svIn_bDf10_allTides_windOn_28112007_5wks/'
simOutputDir = '/home/cthomas/GBR_slim/Saved_Data/gbr460KTests/460K_tOnwCSFRbNCJ0p8SvSth_FWind1o50Every_bDn016Reefs10'
#LP_OutputDir = '/scratch/cthomas/temp'#simOutputDir +'/LPT_Output/SettProg_Mort007'
LP_OutputDir = 'temp'
##
# Add /CPUX to LP_OutputDir if -c option used
cpuExt=""
if ( len(sys.argv)>1 ):
if sys.argv[1] == "-c":
cpuExt = "/CPU"+sys.argv[2]
time.sleep(float(sys.argv[2])*2) #Avoid simultaneous simulations crashing each other
# Define simOutput directory & particle concentration
simOutputDir = '/scratch/cthomas/Hydro_Realistic/gbr500K_sma_1svIn_bDf10_allTides_windOn_28112007_5wks/'
#simOutputDir = 'Output10K_4wks_windsOn_allTides_16112000'
#LP_OutputDir = '/scratch/cthomas/temp'#simOutputDir +'/LPT_Output/SettProg_Mort007'
LP_OutputDir = 'temp'
LP_OutputDir = LP_OutputDir + cpuExt
if(not os.path.exists(LP_OutputDir)):
try : os.mkdir(LP_OutputDir)
except: 0;
##
#-------------------------------------------------------------------------------
##
concentration = 2.0 #No. of particles to seed per kmsq
#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"
reefMapShallowFilenameTrunk = "GBR_ReefMap_2966_Solid9200"
reefMapDeepFilenameTrunk = "GBR_DeepReefMap_2966_Solid9200"
reefMapShallowFilename = reefMapShallowFilenameTrunk+".dat"
reefMapDeepFilename = reefMapDeepFilenameTrunk+".dat"
# Radius: 6.37101e6 for old meshes (sma_zoomed), else 6.371e6 (used for reefsObj and reefsDistanceObj)
EarthRadius=6.37101e6
# Radius: 6.37101e6 for old meshes (sma_zoomed), else 6.371e6 (used for reefsObj and reefsDistanceObj, the former only if useProj=0)
EarthRadius=6.371e6
# Set whether 1: to use PROJ projection library (set parameters in GBR.cc) or 0: to use old projection
usePROJ = 1
# If using PROJ, set projection systems:
......@@ -62,19 +70,6 @@ lonlat_projection = "+proj=latlong +ellps=WGS84"
mesh_projection = "+proj=utm +ellps=WGS84 +zone=55 +south"
##
# SELECT REEF TYPES TO SEED AND SETTLE OVER
#0: shallow reefs, 1: deep reefs, 2: both
seedOverWhichReefs = 2
settleOverWhichReefs = 2
##
LP_OutputDir = LP_OutputDir + cpuExt
if(not os.path.exists(LP_OutputDir)):
try : os.mkdir(LP_OutputDir)
except: 0;
XYZ = function.getCoordinates()
print(colored('Loading simulation parameters ...',"red"))
simInfoFile = open(simOutputDir+'/simulationInfo.txt','r')
simInfo = simInfoFile.readlines()
......@@ -98,7 +93,9 @@ if(not os.path.exists('./LPT_Functions')):
except: 0;
simSolutions = [simExportGap*solutionNumber for solutionNumber in range(0,simTotIter/simExportGap+1)]
#-------------------------------------------------------------------------------
#--------- SET LPT TIME PARAMETERS ---------------------------------------------
print(colored('Loading particle-tracking parameters ...',"red"))
LP_dt = 90. #(double) in seconds
LP_TotTime = 4*7*24.0 #(double) in hours
......@@ -108,8 +105,21 @@ LP_TotIter = int(LP_TotTime*3600.0/LP_dt)
LP_StartTime = simStartTime + (7.0*24.0*3600.0) #in seconds
LP_EndTime = LP_StartTime + LP_TotTime*3600.0
simItersPerLPIter = LP_dt/simDt
#-------------------------------------------------------------------------------
# REEF SEEDING PARAM
#--------- 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
# SET SEEDING CONCENTRATION (nb particles/kmsq)
concentration = 2.0
# 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
seedOverInterval = 24.0*3600.0
seedNbTimes = 4
# WHICH REEFS SHALL WE SEED?
# -> 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)
......@@ -164,7 +174,9 @@ swimParams.set(0,1, 0.1)
swimParams.set(0,2, 2.0)
#[3]: swimming starts after .. hours
swimParams.set(0,3, 0.0)
#-------------------------------------------------------------------------------
#--------- LOAD MODEL & BATHYMETRY ---------------------------------------------
print(colored('Loading model ...',"red"))
model = GModel()
if usePROJ == 0:
......@@ -174,7 +186,6 @@ else:
print " (Using PROJ for projections)"
model.load("./Meshes/"+simMesh+"_utm.msh")
groups = dgGroupCollection(model, 2, 1) #model, dimension, order
#groups.buildGroupsOfInterfaces()
XYZ = function.getCoordinates()
XYZ_DC = dgDofContainer(groups, 3)
XYZ_DC.interpolate(XYZ)
......@@ -187,8 +198,9 @@ if(os.path.exists("./Bath/"+simMesh+"_bath_smooth/"+simMesh+"_bath_smooth.idx"))
bathDC.readMsh("./Bath/"+simMesh+"_bath_smooth/"+simMesh+"_bath_smooth.idx")
else:
print(colored('**** ERROR: Bathymetry file does not exist! ****',"red"))
#-------------------------------------------------------------------------------
#------- SET UP DIFFUSIVITY ----------------------------------------------------
print(colored('Defining diffusivity ...',"red"))
CCode = """
......@@ -231,8 +243,9 @@ dofMan = dgDofManager.newCG (groups, 1, sys)
diffProjector = L2ProjectionContinuous(dofMan)
diffProjector.apply(diffDC,diffFC)
diffDC.exportMsh("./LPT_Functions/DIFFUSIVITY_smooooth_"+simMesh,0,0)
#-------------------------------------------------------------------------------
# Load reef maps
#------- LOAD REEF MAPS --------------------------------------------------------
reefsShallowDC = dgDofContainer(groups,1)
reefsDeepDC = dgDofContainer(groups,1)
reefsShallowDCfilename = "LPT_Functions/REEFS_RAW_LPT_"+reefMapShallowFilenameTrunk+"_"+simMesh
......@@ -333,7 +346,6 @@ else:
#dir2RF = directionToReefsDC.getFunction()
#dirFunction = functionC(glib,"current2",3,[directionToReefsDC.getFunction()])
#directionToReefsDC.exportMshNodeData("./DIRECTION_TO_REEFS",0,0, '')
nbOfShallowReefsFConst = functionConstant(nbOfShallowReefs)
......@@ -341,21 +353,36 @@ reefsDualDC = dgDofContainer(groups, 1)
reefsDualFC = functionC(glib, "overwrite2", 1, [bathDC.getFunction(), reefsShallowDC.getFunction(), reefsDeepDC.getFunction(), nbOfShallowReefsFConst])
reefsDualDC.interpolate(reefsDualFC)
reefsDualDC.exportMsh("LPT_Functions/REEFS_RAW_LPT_DUAL_"+simMesh,0,0)
#-------------------------------------------------------------------------------
#--------- SEED PARTICLES ------------------------------------------------------
print(colored('Seeding particles ...',"red"))
secondsBetweenSeeds = 0
itersBetweenSeeds = 0
if seedOverTime == 1:
secondsBetweenSeeds = seedOverInterval/float(seedNbTimes)
itersBetweenSeeds = secondsBetweenSeeds/LP_dt
concentration = concentration/float(seedNbTimes)
print' We are seeding particles PROGRESSIVELY OVER TIME:'
print' -> Seeding',seedNbTimes,'times, at intervals of', secondsBetweenSeeds/3600.0, 'hours ('+str(itersBetweenSeeds)+' iters)'
print' -> Total concentration:', concentration*seedNbTimes, ', concentration per seeding event:', concentration
# ***SEED OVER ALL REEFS***
#initCM = fullMatrixDouble( nbOfReefs, nbOfReefs )
#particles = particleArray(groups, reefsDualDC, initCM, reefsToSeed, concentration, seedOverWhichReefs, nbOfShallowReefs)
#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' )
# ******
# ***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)
......@@ -369,13 +396,15 @@ for i in range(0,len(reefPointsLatLon)):
reefPointsCart[i][0] = dest_points(0,0)
reefPointsCart[i][1] = dest_points(1,1)
particles = particleArray(groups, 0, 0, 0, 0, 0)
particles = particleArray()
for location in range(0,len(reefPointsCart)):
particles.addParticles(groups, 100, reefPointsCart[location][0], reefPointsCart[location][1], 0, reefPointsCart[location][3])
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.printPositions(LP_OutputDir+'/seededParticles','dat',0)
# ******
#-------------------------------------------------------------------------------
#-------- PRINT SIMULATION INFO ------------------------------------------------
print ''
print '-----------------------'
print 'HYDRODYNAMICS PARAMETERS:'
......@@ -424,26 +453,42 @@ print '-----------------------'
print ''
print 'Output directory:', LP_OutputDir
print ''
#-------------------------------------------------------------------------------
#------- SET UP HYDRO DOFCONTAINERS --------------------------------------------
simSolution = dgDofContainer(groups, 3)
simSolutionUnder = dgDofContainer(groups, 3)
simSolutionOver = dgDofContainer(groups, 3)
#Keep track of the solutions loaded into the three DCs:
loadedSolutions = [-1] * 3 #[0]:simSolution, [1]:simSolutionOver, [2]:simSolutionUnder
connectivityMatrix = fullMatrixDouble( nbOfReefs, nbOfReefs )
dummyPrinter = connectivityMatrixPrint()
#-------------------------------------------------------------------------------
#------ SET UP PARTICLE TRACKER OBJECT -----------------------------------------
#Create particle tracker object
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)
#-------------------------------------------------------------------------------
startcpu = time.clock()
#Keep track of the solutions loaded into the three DCs:
loadedSolutions = [-1] * 3 #[0]:simSolution, [1]:simSolutionOver, [2]:simSolutionUnder
#-------- PARTICLE LOOP --------------------------------------------------------
print 'Beginning loop ...'
nextIterToSeed=itersBetweenSeeds
seedNbTimes = seedNbTimes-1 #seedNbTimes now counts down how many more times we need to seed
startcpu = time.clock()
for n in range(1,LP_TotIter+1):
#Check whether to seed more particles
if seedOverTime == 1 and seedNbTimes > 0 and n >= nextIterToSeed:
print '*** Seeding more particles ... ***'
particles.addParticlesWithReefMap(groups, reefsDualDC, initCM, reefsToSeed, concentration, seedOverWhichReefs, nbOfShallowReefs)
nextIterToSeed = nextIterToSeed + itersBetweenSeeds
seedNbTimes = seedNbTimes - 1
print '***', seedNbTimes,'seeding events left ***'
t = simStartTime + LPsim_TimeOffset + float(n)*LP_dt
iterNumberWanted = int(LPsim_TimeOffset/simDt + n*simItersPerLPIter)
print ''
......@@ -550,7 +595,7 @@ for n in range(1,LP_TotIter+1):
dummyPrinter.CMPrinter(connectivityMatrix, LP_OutputDir+'/connectivityMatrix%06d'%( n ) )
else:
dummyPrinter.CMPrinter(connectivityMatrix, LP_OutputDir+'/settlementMatrix%06d'%( n ) )
#-------------------------------------------------------------------------------
print(colored('Particle tracking finished. Closing module.',"yellow"))
Msg.Exit(0)
......@@ -357,9 +357,9 @@ void merge3D_with_UV(dataCacheMap *,fullMatrix<double> &sol, fullMatrix<double>
//Function to merge shallow (hard) and deep (soft) DCs
void overwrite2(dataCacheMap *,fullMatrix<double> &sol, fullMatrix<double> &bath, fullMatrix<double> &shallowDC, fullMatrix<double> &deepDC, fullMatrix<double> &nbOfShallowReefs){
for (size_t i = 0; i < sol.size1(); i++) {
if (shallowDC(i,0) > 0.1 && (bath(i,0) < 20.0)) sol.set(i, 0, shallowDC(i,0));
if (shallowDC(i,0) > 0.1) sol.set(i, 0, shallowDC(i,0));
else {
if ((deepDC(i,0) > 0.1) && (bath(i,0) > 20.0) && (bath(i,0) < 60.0)) sol.set(i, 0, deepDC(i,0)+nbOfShallowReefs(0,0));
if (deepDC(i,0) > 0.1) sol.set(i, 0, deepDC(i,0)+nbOfShallowReefs(0,0));
else sol.set(i, 0, 0.0);
}
}
......@@ -375,7 +375,6 @@ 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++) {
......
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