dg: LPT fix

git-svn-id: https://geuz.org/svn/gmsh/trunk@16299 59c8cabf-10d6-4a28-b43f-3936b85956c4

benchmarks/GreatBarrierReef/LPT_Caller_DualReefMap.py

@@ -40,8 +40,8 @@ if ( len(sys.argv)>1 ):
 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 = '/scratch/cthomas/temp'#simOutputDir +'/LPT_Output/SettProg_Mort007'
+LP_OutputDir = 'temp'
 
 
 ##
@@ -55,22 +55,19 @@ reefMapDeepFilename = reefMapDeepFilenameTrunk+".dat"
 
 # Radius: 6.37101e6 for old meshes (sma_zoomed), else 6.371e6 (used for reefsObj and reefsDistanceObj)
 EarthRadius=6.37101e6
-# Set whether 1: to use PROJ projeciton library (set parameters in GBR.cc) or 0: to use old projection
+# 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:
 lonlat_projection = "+proj=latlong +ellps=WGS84"
-mesh_projection = "+proj=utm +ellps=WGS84 +zone=55"
+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 = 1
+seedOverWhichReefs = 2
 settleOverWhichReefs = 2
 ##
 
-if(not os.path.exists(LP_OutputDir)):
-  try : os.mkdir(LP_OutputDir)
-  except: 0;
 LP_OutputDir = LP_OutputDir + cpuExt
 if(not os.path.exists(LP_OutputDir)):
   try : os.mkdir(LP_OutputDir)
@@ -347,23 +344,35 @@ reefsDualDC.exportMsh("LPT_Functions/REEFS_RAW_LPT_DUAL_"+simMesh,0,0)
 
 print(colored('Seeding particles ...',"red"))
 # ***SEED OVER ALL REEFS***
-initCM = fullMatrixDouble( nbOfReefs, nbOfReefs )
-particles = particleArray(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(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)
-#cartCoords = [0]*3
-#geoCoordTo3DCartesian(-19.9400,149.1100,cartCoords,EarthRadius)
-#projectPoint(cartCoords)
-#particles = particleArray(groups, 0, cartCoords[0], cartCoords[1], 0, 0)
-#for location in range(0,10):
-  #particles.addParticles(groups, 1000, cartCoords[0], cartCoords[1], 0, location)
-#particles.printPositions(LP_OutputDir+'/seededParticles','pos',0)
+
+UTMProjectorLonLatToUTM = slimProjector(lonlat_projection, mesh_projection)
+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(groups, 0, 0, 0, 0, 0)
+for location in range(0,len(reefPointsCart)):
+  particles.addParticles(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)
 
 
@@ -425,7 +434,10 @@ dummyPrinter = connectivityMatrixPrint()
 
 #Create particle tracker object
 print 'Initialising particle tracker ...'
-particleTracker = dgParticleTracker2D(groups, particles, bathDC, diffDC, reefsDualDC, distToReefsDC, directionToReefsDC, settleParams, mortParams, swimParams, nbOfShallowReefs, LP_TotIter, LP_dt, LP_OutputDir)
+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:

benchmarks/GreatBarrierReef/Libraries/GBR.cc

@@ -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) sol.set(i, 0, shallowDC(i,0));
+    if (shallowDC(i,0) > 0.1 && (bath(i,0) < 20.0)) sol.set(i, 0, shallowDC(i,0));
     else {
-      if ((deepDC(i,0) > 0.1) && (bath(i,0) > 6.0)) sol.set(i, 0, deepDC(i,0)+nbOfShallowReefs(0,0));
+      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));
       else sol.set(i, 0, 0.0);
     }
   }
@@ -440,8 +440,23 @@ void bottomDrag_const_n008(dataCacheMap *, fullMatrix<double> &val, fullMatrix<d
 
 void bottomDrag_const_n016(dataCacheMap *, fullMatrix<double> &val, fullMatrix<double> &sol, fullMatrix<double> &bath){
   for(int i=0;  i<val.size1(); i++){
+    double depthFactor=1.0;
     double H = sol(i, 0)+bath(i, 0);
-    val.set(i, 0, 9.81*0.016*0.016/(pow(H, 1.333333333333)));
+//     if H > 200 depthFactor = 10.0;
+//     if H > 500 depthFactor = 100.0;
+    val.set(i, 0, 9.81*depthFactor*0.016*0.016/(pow(H, 1.333333333333)));
+  }
+}
+
+void bottomDrag_reefs_n016(dataCacheMap *, fullMatrix<double> &val, fullMatrix<double> &sol, fullMatrix<double> &bath, fullMatrix<double> &reefs){
+  for(int i=0;  i<val.size1(); i++){
+    double reef = 1.0;
+    double H = sol(i, 0)+bath(i, 0);
+    //Verify if there is a reef or not.
+    //  Reef --> multiply normal bottom drag by factor 10x
+    int reefYesOrNo = int(reefs(i, 0));
+    if (reefYesOrNo > 0 && bath(i,0) < 25.0) reef = 10.0;
+    val.set(i, 0, 9.81*reef*0.016*0.016/(pow(H, 1.333333333333)));
   }
 }
 
@@ -623,7 +638,7 @@ void functionCurrentAsFWind(dataCacheMap *, fullMatrix<double> &sol, fullMatrix<
 
 void functionCurrentAsFWind_calcDepthFactor(dataCacheMap *, fullMatrix<double> &sol, fullMatrix<double> &bath) {
   for (size_t i=0; i<sol.size1(); i++) {
-    double factor = 200.0/bath(i,0);
+    double factor = 100.0/bath(i,0);
     if (factor > 2.0) factor = 2.0;
     sol.set(i, 0, factor);
   }

dgParticleTracker2D.h

@@ -364,6 +364,8 @@ public:
   /** Activate particle swimming and load parameters */
   void addSwimming(fullMatrix<double>* swimParams, dgDofContainer* reefs, dgDofContainer* distToReefs, dgDofContainer* directionToReefs) {
     _swimFlag = true;
+    //TEMP: disactivate swimming if swimParams(0,0)=0.0
+    if (swimParams < 0.1) _swimFlag=false;
     //Load swimming parameters:
     _swimType = static_cast<int>(swimParams->get(0,0)+0.1);
     _swimSpeed = swimParams->get(0,1);