#!/usr/bin/env python # -*- coding: utf-8 -*- ## # @package particles # Lagrangian particle-tracker # # @author Thomas, Christopher # @author Vallaeys, Valentin # @author Le Bars, Yoann # @version 1.0 # @date 2015/01/19 # @date 2015/01/23 from dgpy import * import time import os.path import subprocess import Common from 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('output/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 = 96.0 #(double) in hours LP_Export = int( (1.0*(3600.0-90.0))/LP_dt ) #(int) export particle positions every ... steps LP_StartTime = simStartTime + 48 * 3600 LP_TotIter = int(LP_TotTime*3600.0/LP_dt) LP_EndTime = LP_StartTime + LP_TotTime*3600.0 simItersPerLPIter = LP_dt/simDt #------------------------------------------------------------------------------- #---------- LOAD MESH & BATHYMETRY ---------------------------------- print(colored('Loading mesh and bathymetry ...',"red")) model = GModel() # Loading the mesh. model.load("../data/mesh/congoCoast-labels-clean_utm.msh") groups = dgGroupCollection(model) groups.splitGroupsByPhysicalTag() # Bathymetry bathDC = dgDofContainer(groups, 1) bathDC.importIdx("Bath/congoCoast-labels-clean_bath_smooth/congoCoast-labels-clean_bath_smooth.idx") #------------------------------------------------------------------------------- #------- 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")) #seedPoints = [[0, 0], [-100000,-100000], [-100000,100000], [100000,-100000], [100000,100000], [LengthOfDomain-1, LengthOfDomain-1] ] seedPointsLL = [[12.7, -6.]]#, [12.25, -6]]#[[ 7., -6.,], [ 8., -6.]] seedPointsll = seedPointsLL for location in range(0,len(seedPointsLL)): seedPointsll[location][0] = seedPointsLL[location][0]*np.pi/180. seedPointsll[location][1] = seedPointsLL[location][1]*np.pi/180. latlon_pj = pj_init_plus("+proj=latlong +ellps=WGS84") utm_pj = pj_init_plus("+proj=utm +ellps=WGS84 +zone=32") seedPointsUTM = pjTransform(latlon_pj, utm_pj, seedPointsLL) seedPoints = seedPointsLL for location in range(0,len(seedPointsLL)): seedPoints[location][0] = float(seedPointsUTM[location][0]) seedPoints[location][1] = float(seedPointsUTM[location][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)' ) print ('Simulation starts:', simStartTime/3600.0, '(hr since origin) | Ends:', (simStartTime+simTotIter*simDt)/3600.0, '(hr since origin)' ) print ('Data exported every', simExportGap, '(iter) |', simExportGap*simDt/60.0, '(mins)') print ('') print ('-----------------------') print ('PARTICLE-TRACKING PARAMETERS:') print ('Dt:', LP_dt, 's') 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 particles:', particles.printNbParticles()) print ('-----------------------') print ('') LPsim_TimeOffset = LP_StartTime - simStartTime print ('-----------------------') print ('Particle-tracker / Simulation time offset is:', LPsim_TimeOffset/(24.0*3600.0),'(days)') print ('-----------------------') print ('') print ('Output directory:', LP_OutputDir) print ('') #------------------------------------------------------------------------------- #------- SET UP HYDRO DOFCONTAINERS -------------------------------------------- simSolution = dgDofContainer(groups, 3) simSolutionUnder = dgDofContainer(groups, 3) simSolutionOver = dgDofContainer(groups, 3) simSolutions = [simExportGap*solutionNumber for solutionNumber in range(0,int(simTotIter/simExportGap)+1)] loadedSolutions = [-1] * 3 #------------------------------------------------------------------------------- #------ 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) #------------------------------------------------------------------------------- #------ SET UP DUMMY CONNECTIVITY MATRIX --------------------------------------- # Required argument of particleMove: NxN matrix connectivityMatrix = fullMatrixDouble(1,1) #------------------------------------------------------------------------------- #-------- PARTICLE LOOP -------------------------------------------------------- print ('Beginning loop ...') startcpu = time.clock() for n in range(1,LP_TotIter+1): t = simStartTime + LPsim_TimeOffset + float(n)*LP_dt iterNumberWanted = int(LPsim_TimeOffset/simDt + n*simItersPerLPIter) print ('') print('LPT iteration:', n, 'of',LP_TotIter,'| simIteration wanted:',iterNumberWanted, 'for t = %.2f'%( (t-simStartTime)/3600.0 )+ ' (hrs) | Loaded solutions:',loadedSolutions[0],loadedSolutions[1],loadedSolutions[2],'| CPU time:', time.clock()-startcpu) if iterNumberWanted in simSolutions: if iterNumberWanted in loadedSolutions: if iterNumberWanted == loadedSolutions[0]: print('-> Exact solution for simSolution (iteration %06d'%( iterNumberWanted )+') already loaded in simSolution. Using this.') elif iterNumberWanted == loadedSolutions[1]: print('-> Exact solution for simSolution (iteration %06d'%( iterNumberWanted )+') already loaded in simSolutionOver. Using this.') simSolution.copy(simSolutionOver) loadedSolutions[0] = iterNumberWanted elif iterNumberWanted == loadedSolutions[2]: print('-> Exact solution for simSolution (iteration %06d'%( iterNumberWanted )+') already loaded in simSolutionUnder. Are we going back in time??? Using this anyway...') simSolution.copy(simSolutionUnder) loadedSolutions[0] = iterNumberWanted else: print('Problem in loadedSolutions!') else: print('-> simSolution not already loaded. Loading it from file. File is (gbr-%06d'%( iterNumberWanted )+').') simSolution.importIdx('output/idx/idx-%06d'%( iterNumberWanted / simExportGap )+'.idx') loadedSolutions[0] = iterNumberWanted else: l=0 while (not iterNumberWanted-l in simSolutions): l=l+1 if iterNumberWanted-l in loadedSolutions: if iterNumberWanted-l == loadedSolutions[0]: print('-> Exact solution for simSolutionUnder (iteration %06d'%( iterNumberWanted-l )+') already loaded in simSolution. Using this.') simSolutionUnder.copy(simSolution) loadedSolutions[2] = iterNumberWanted-l t_solUnder = simStartTime + float(iterNumberWanted-l)*simDt elif iterNumberWanted-l == loadedSolutions[1]: print('-> Exact solution for simSolutionUnder (iteration %06d'%( iterNumberWanted-l )+') already loaded in simSolutionOver. Using this.') simSolutionUnder.copy(simSolutionOver) loadedSolutions[2] = iterNumberWanted-l t_solUnder = simStartTime + float(iterNumberWanted-l)*simDt elif iterNumberWanted-l == loadedSolutions[2]: print('-> Exact solution for simSolutionUnder (iteration %06d'%( iterNumberWanted-l )+') already loaded in simSolutionUnder. Using this.') else: print('Problem in loadedSolutions!') else: print('-> simSolutionUnder not already loaded. Loading it from file. l is:', l, '. File is (gbr-%06d'%( iterNumberWanted-l )+').') t_solUnder = simStartTime + float(iterNumberWanted-l)*simDt simSolutionUnder.importIdx('output/idx/idx-%06d'%( (iterNumberWanted-l) / simExportGap )+'.idx') loadedSolutions[2] = iterNumberWanted-l m=0 while (not iterNumberWanted+m in simSolutions): m=m+1 if iterNumberWanted+m in loadedSolutions: if iterNumberWanted+m == loadedSolutions[0]: print('-> Exact solution for simSolutionOver (iteration %06d'%( iterNumberWanted+m )+') already loaded in simSolution. Using this.') simSolutionOver.copy(simSolution) loadedSolutions[1] = iterNumberWanted+m t_solOver = simStartTime + float(iterNumberWanted+m)*simDt elif iterNumberWanted+m == loadedSolutions[1]: print('-> Exact solution for simSolutionOver (iteration %06d'%( iterNumberWanted+m )+') already loaded in simSolutionOver. Using this.') elif iterNumberWanted+m == loadedSolutions[2]: print('-> Exact solution for simSolutionOver (iteration %06d'%( iterNumberWanted+m )+') already loaded in simSolutionUnder. Are we going back in time??? Using this anyway...') simSolutionOver.copy(simSolutionUnder) loadedSolutions[1] = iterNumberWanted+m t_solOver = simStartTime + float(iterNumberWanted+m)*simDt else: print('Problem in loadedSolutions!') else: print('-> simSolutionOver not already loaded. Loading it from file. m is:', m, '. File is (gbr-%06d'%( iterNumberWanted+m )+').') t_solOver = simStartTime + float(iterNumberWanted+m)*simDt simSolutionOver.importIdx('output/idx/idx-%06d'%( (iterNumberWanted+m) / simExportGap )+'.idx') loadedSolutions[1] = iterNumberWanted+m scaleFactor = 1.0 - (t-t_solUnder)/(t_solOver-t_solUnder) print(' Info: simSolutionUnder scale factor is', scaleFactor) simSolution.copy(simSolutionUnder) #Put simSolUnder in simSolution simSolution.scale(scaleFactor) #Scale simSolUnder (dof * relative dist. in t of sol. to simSolUnder) simSolution.axpy(simSolutionOver,1.0-scaleFactor) #Scale simSolOver and add to scaled simSolUnder loadedSolutions[0] = -1 #Always give it same solution particleTracker.particleMove(simSolution, connectivityMatrix, n) # Print output: if ( n%(LP_Export) == 0 ): # Export .pos with positions particles.printPositions(LP_OutputDir+'/particlesAlive_%06d'%( n ), 'pos', 0) # Export .dat with positions particles.printPositions(LP_OutputDir+'/particleOutput_%06d'%( n ), 'dat', 0) #------------------------------------------------------------------------------- print(colored('Particle tracking finished. Closing module.',"yellow")) Msg.Exit(0)