Develop

SRC/EXTENSIONS/GAD/CUDA/cuda_GADDevice.cu

@@ -236,7 +236,7 @@ extern "C" int cuda_GADDeviceCleanup(){
      cudaFree(GAD_turbineVseries_d);
      cudaFree(u_sampAvg_d);
      cudaFree(v_sampAvg_d);
-     cudaFree(GAD_yawError);
+     cudaFree(GAD_yawError_d);
      cudaFree(GAD_anFactor_d);
      cudaFree(GAD_Xcoords_d);
      cudaFree(GAD_Ycoords_d);

SRC/EXTENSIONS/GAD/GAD.c

@@ -665,11 +665,12 @@ int GADConstructor(){
 int GADInitTurbineRefChars(float dt){
   int errorCode = GAD_SUCCESS;
   int iturb,i,j,k;
-  int ijk;
+  int ijk,ijkp1,ijkm1;
   int ij;
   float rVec;
   float rVec0;
   float deltaz, deltaz0;
+  float deltaz_p1, deltaz_m1;
 
   /*Initialize requisite parameters for the RefMag and RefDir calculations*/
   GADsamplingAvgLength = (int) floor(GADrefSampleWindow/dt);    //Determine the number of model timesteps in a sample window (high frequencies filter)
@@ -702,30 +703,39 @@ int GADInitTurbineRefChars(float dt){
     }//end if inFile == NULL
     for(i=iMin-Nh; i < iMax+Nh; i++){
       for(j=jMin-Nh; j < jMax+Nh; j++){
-        for(k=kMin-Nh; k < kMax+Nh; k++){
+        for(k=kMin-Nh+1; k < kMax+Nh-1; k++){
            ijk = i*(Nyp+2*Nh)*(Nzp+2*Nh)+j*(Nzp+2*Nh)+k;
            ij = i*(Nyp+2*Nh)+j;
+           ijkp1 = i*(Nyp+2*Nh)*(Nzp+2*Nh)+j*(Nzp+2*Nh)+(k+1);
+           ijkm1 = i*(Nyp+2*Nh)*(Nzp+2*Nh)+j*(Nzp+2*Nh)+(k-1);
            rVec = sqrt( pow((GAD_Xcoords[iturb]-xPos[ijk]),2.0)
                        +pow((GAD_Ycoords[iturb]-yPos[ijk]),2.0));
-           if(rVec <= sqrt(pow(dX,2.0)+pow(dY,2.0))){ //Should be a candiate gridcell for (nacelle center) reference location
-	     if(rVec <= rVec0){
-	       if(rVec < rVec0){
+           if((rVec <= sqrt(pow(dX,2.0)+pow(dY,2.0))) && (rVec <= rVec0)){ //Should be a candiate gridcell for (nacelle center) reference location
+	     if(rVec < rVec0){
+	         printf("%d/%d: iturb = %d, rVec = %.9f, rVec0 = %.9f\n",
+                      mpi_rank_world,mpi_size_world, iturb, rVec, rVec0);
 	         rVec0 = rVec;
 	         GAD_turbineRank[iturb] = mpi_rank_world; 
                  GAD_turbineRefi[iturb] = i;
 	         GAD_turbineRefj[iturb] = j;
-	       }
-	       deltaz = sqrt(pow((GAD_hubHeights[GAD_turbineType[iturb]]-(zPos[ijk]-topoPos[ij])),2.0));
+	     }
+	     if ((GAD_turbineRank[iturb]==mpi_rank_world) && (GAD_turbineRefi[iturb]==i) && (GAD_turbineRefj[iturb]==j)){
+	       deltaz = fabsf(GAD_hubHeights[GAD_turbineType[iturb]]-(zPos[ijk]-topoPos[ij]));
+	       deltaz_p1 = GAD_hubHeights[GAD_turbineType[iturb]]-(zPos[ijkp1]-topoPos[ij]);
+	       deltaz_m1 = GAD_hubHeights[GAD_turbineType[iturb]]-(zPos[ijkm1]-topoPos[ij]);
 #ifdef DEBUG_TURBCHAR
-               printf("%d/%d: deltaz = %f, 0.5/(J33[ijk]*dZi)) = %f\n",
-                      mpi_rank_world,mpi_size_world, deltaz, 0.5/(J33[ijk]*dZi));
+               printf("%d/%d: iturb = %d, deltaz = %.9f, deltaz0 = %.9f, deltaz_m1 = %.9f, deltaz_p1 = %.9f\n",
+                      mpi_rank_world,mpi_size_world, iturb, deltaz, deltaz0, deltaz_m1, deltaz_p1);
 #endif
-               if(deltaz <= 0.5/(J33[ijk]*dZi)){         // 1/(J33[ijk]*dZi)) = dz of cell
+               if((deltaz_m1 > 0.0) && (deltaz_p1 < 0.0) && (deltaz < deltaz0)){
 		 deltaz0 = deltaz;      
 		 GAD_turbineRefk[iturb] = k;      
+#ifdef DEBUG_TURBCHAR
+		 printf("%d/%d: iturb = %d, k = %d, j=%d, i=%d, deltaz = %.9f, deltaz0 = %.9f\n", mpi_rank_world,mpi_size_world, iturb, k, j, i, deltaz, deltaz0); 
+#endif
                }//end if vertical delta < dz ...
-             }//end if rVec < rVec0...
-           }//end if rVec...
+	     }// end if GAD_turbine...
+           }//end if rVec<=...
         } //end for(k...
       } // end for(j...
     } // end for(i...

SRC/EXTENSIONS/URBAN/CUDA/cuda_urbanDevice.cu

@@ -41,6 +41,7 @@ extern "C" int cuda_urbanDeviceSetup(){
 
    cudaMemcpyToSymbol(delta_aware_bdg_d, &delta_aware_bdg, sizeof(float));
 
+   cudaMemcpyToSymbol(urban_heatRedis_d, &urban_heatRedis, sizeof(int));
    if(urban_heatRedis > 0){
      Nelems = (Nxp+2*Nh)*(Nyp+2*Nh);
      fecuda_DeviceMalloc(Nelems, &urban_heat_redis_d);
@@ -157,7 +158,7 @@ __global__ void cudaDevice_URBANfinal(float* hydroFlds_d, float* hydroFldsFrhs_d
       ijk = i*iStride + j*jStride + k*kStride;
 
       cudaDevice_UrbanDragMethod(&hydroFlds_d[fldStride*RHO_INDX+ijk],&hydroFlds_d[fldStride*U_INDX+ijk],&hydroFlds_d[fldStride*V_INDX+ijk],&hydroFlds_d[fldStride*W_INDX+ijk],
-                                 &hydroFlds_d[fldStride*THETA_INDX+ijk],&hydroBaseStateFlds_d[fldStride*THETA_INDX+ijk],&hydroBaseStateFlds_d[fldStride*RHO_INDX+ijk],
+                                 &hydroFlds_d[fldStride*THETA_INDX+ijk],&hydroBaseStateFlds_d[fldStride*THETA_INDX_BS+ijk],&hydroBaseStateFlds_d[fldStride*RHO_INDX_BS+ijk],
                                  &hydroFldsFrhs_d[fldStride*U_INDX+ijk],&hydroFldsFrhs_d[fldStride*V_INDX+ijk],&hydroFldsFrhs_d[fldStride*W_INDX+ijk],
                                  &hydroFldsFrhs_d[fldStride*THETA_INDX+ijk],&hydroFldsFrhs_d[fldStride*RHO_INDX+ijk],&building_mask_d[ijk]);
       if(NhydroAuxScalars_d > 0){

SRC/FEMAIN/FastEddy.c

@@ -63,13 +63,12 @@ int main(int argc, char **argv){
      /* Parse the command line arguments */
      if(argc != 2){
         printf("usage: %s paramFile \n",argv[0]);
-        fflush(stdout);
         exit(0);
      }else{
         sscanf(argv[1],"%s", paramFile);
         printf("Obtaining parameters from %s\n", paramFile);
-        fflush(stdout);
      }
+     fflush(stdout);
   } //end if(mpi_rank == 0 )
 
   /*** ---------------------------------------------------------------------------------------------- ***/
@@ -260,6 +259,56 @@ int main(int argc, char **argv){
     errorCode = hydro_corePrepareFromInitialConditions(simTime_itRestart, dt);
   }//end if inFile !=NULL
 
+
+  /*** ---------------------------------------------------------------------------------------------- ***/
+  /*** ----------------- Initialize/configure any specified Profile/Planes IO functionality ----------***/
+  /*** ---------------------------------------------------------------------------------------------- ***/
+  if(towerIOSelector > 0){
+    int iprofile;
+    int tmp_rank;
+    errorCode = ioProfilePreparations();
+    for(iprofile = 0; iprofile < nProfs; iprofile++){
+      if(towerProfiles.coordType == 0){
+        tmp_rank = gridGetRankFromLatLonPosition(towerProfiles.coordsLon[iprofile],towerProfiles.coordsLat[iprofile]);
+        if(tmp_rank >= 0){
+          towerProfiles.mpi_ranks[iprofile] = tmp_rank;
+          printf("Rank %d/%d: profile ID = %d at position (lat,lon) = (%f,%f), found in mpi_rank = %d subdomain!\n",
+                 mpi_rank_world, mpi_size_world, towerProfiles.profIDs[iprofile],towerProfiles.coordsLat[iprofile],towerProfiles.coordsLon[iprofile],
+ 	         towerProfiles.mpi_ranks[iprofile]);
+        }else{
+          printf("Rank %d/%d: profile ID = %d at position (lat,lon) = (%f,%f), not in simulation domain!\n",
+                 mpi_rank_world, mpi_size_world, towerProfiles.profIDs[iprofile],towerProfiles.coordsLat[iprofile],towerProfiles.coordsLon[iprofile]);
+        }
+      }else{
+        tmp_rank = gridGetRankFromXYPosition(towerProfiles.coordsWE[iprofile],towerProfiles.coordsSN[iprofile]);
+        if(tmp_rank >= 0){
+          towerProfiles.mpi_ranks[iprofile] = tmp_rank;
+          printf("Rank %d/%d: profile ID = %d at position (x,y) = (%f,%f), found in mpi_rank = %d subdomain!\n",
+                 mpi_rank_world, mpi_size_world, towerProfiles.profIDs[iprofile],towerProfiles.coordsWE[iprofile],towerProfiles.coordsSN[iprofile],
+ 	         towerProfiles.mpi_ranks[iprofile]);
+        }else{
+          printf("Rank %d/%d: profile ID = %d at position (x,y) = (%f,%f), not in simulation domain!\n",
+                 mpi_rank_world, mpi_size_world, towerProfiles.profIDs[iprofile],towerProfiles.coordsWE[iprofile],towerProfiles.coordsSN[iprofile]);
+        }
+      }
+    }
+    fflush(stdout);
+    errorCode = hydro_coreAllocateTowersDataStructure(nProfs, towerProfiles, NtBatch);
+    MPI_Barrier(MPI_COMM_WORLD); 
+
+    //Write Towers static/initial conditions files
+    ioWriteBinaryTowerInitialFile(dt, simTime_itRestart, Nxp, Nyp, Nzp, Nh, 
+                                  towersData, towersSurfData,
+                                  towerIDs, rank_nTowers, tower_iInds, tower_jInds, 
+	    		          towerProfiles.coordType, tower_xOffsets, tower_yOffsets, tower_LonOffsets, tower_LatOffsets,
+				  NtBatch, towerInstanceSize, towerSurfInstanceSize,
+				  zPos, yPos, xPos, topoPos, surflayer_offshore, sea_mask);
+  }// endif towerIOSelector > 0
+  MPI_Barrier(MPI_COMM_WORLD); 
+  printf("Rank %d/%d: Profile preparations complete!\n",mpi_rank_world, mpi_size_world);
+  fflush(stdout);
+  MPI_Barrier(MPI_COMM_WORLD); 
+
   /*** ---------------------------------------------------------------------------------------------- ***/
   /*** ----------------- Initialize the CUDA-layer of each model-component module --------------------***/
   /*** ---------------------------------------------------------------------------------------------- ***/
@@ -299,8 +348,8 @@ int main(int argc, char **argv){
 
 #endif /* ifndef NOTCUDA: THIS SECTION PREPARED FOR CUDA FASTEDDY SIMULATION */
 
+  fflush(stdout);
   MPI_Barrier(MPI_COMM_WORLD); 
-
   /*** ---------------------------------------------------------------------------------------------- ***/
   /*** ------- Final pre-check logging and initialization before entering the main time-loop ---------***/
   /*** ---------------------------------------------------------------------------------------------- ***/
@@ -359,7 +408,7 @@ int main(int argc, char **argv){
      /*If appropriate timing to do so, update the nesting boundary conditions*/ 
      if(hydroBCs == 1){
        if((it%((int)roundf(dtBdyPlaneBCs/dt))==0)&&(it > simTime_itRestart)){    //If due for an update and after the simulation start
-         printf("FastEddy MAin timestepping loop: Reading new BdyPlanes at it=%d...\n",it);
+         printf("FastEddy Main timestepping loop: Reading new BdyPlanes at it=%d...\n",it);
          fflush(stdout);
          errorCode = timeIntBdyPlaneUpdates();
          if((cellpertSelector==1)&&(cellpert_tvcp==1)){ // update CP parameters with dynamic LBCs
@@ -368,7 +417,9 @@ int main(int argc, char **argv){
        }//end if hydroBCs == 1
      }
      MPI_Barrier(MPI_COMM_WORLD);
+     fflush(stdout);
 
+     mpi_t3 = MPI_Wtime();    //Mark the walltime to measure IO/logging duration.
      if(it%frqOutput == 0){
        MPI_Barrier(MPI_COMM_WORLD); 
        if(mpi_rank_world == 0){
@@ -378,6 +429,7 @@ int main(int argc, char **argv){
        } //if mpi_rank_world
 
        MPI_Barrier(MPI_COMM_WORLD); 
+       fflush(stdout);
 
        /*Every rank calls the StateLogDump*/
        hydro_coreStateLogDump();
@@ -389,7 +441,6 @@ int main(int argc, char **argv){
          fflush(stdout);
        } //if mpi_rank_world
 
-       mpi_t3 = MPI_Wtime();    //Mark the walltime to measure IO duration.
        /* Dump the root output file. */
 #ifndef IO_OFF
        if(ioOutputMode==0){
@@ -405,13 +456,22 @@ int main(int argc, char **argv){
          errorCode = ioWriteBinaryoutFileSingleTime(it, Nxp, Nyp, Nzp, Nh);
 #endif
        }
-#endif
-       mpi_t4 = MPI_Wtime();    //Mark the walltime to measure IO duration
+#endif 
        if(mpi_rank_world == 0){
          printf("Dumped state at timestep = %d...\n",it);
          fflush(stdout);
        } //if mpi_rank_world
      } //end if (it%frqOutput == 0) ....   (We log summary info and dump outputs)
+     //Dump tower data if appropriate
+     if((towerIOSelector > 0) && (it > simTime_itRestart)){
+       ioWriteBinaryTowerFileSingleBatch(it, NtBatch, Nz, simTimeBatch, towersData, towersSurfData, 
+		                         towerIDs, rank_nTowers, towerInstanceSize, towerSurfInstanceSize);
+       if(mpi_rank_world == 0){
+         printf("Dumped batch tower data at timestep = %d...\n",it);
+         fflush(stdout);
+       } //if mpi_rank_world
+     }
+     mpi_t4 = MPI_Wtime();    //Mark the walltime to measure IO/logging duration
 #ifdef NOTCUDA 
      /* OBSELETE!!!!! There is longer any CPU model integration functionality */
 #else  /* ---------------  CUDA FASTEDDY !!!!! -------------------------  */
@@ -429,7 +489,7 @@ int main(int argc, char **argv){
      /*Kernel return*/
      itTmp = itTmp+NtBatch;  
 #endif
-     mpi_t2 = MPI_Wtime();    //Mark the walltime to measure duration of  a batch of timesteps.
+     mpi_t2 = MPI_Wtime();    //Mark the walltime to measure duration of a batch of timesteps (including any IO/logging).
      if(mpi_rank_world == 0){
         printf("\n\t\t\t!!!!!\t  TIMESTEP PERFORMANCE  \t !!!!! \n");
         printf(" Total Time (s) | Batch Steps \t| Time/step (s) | Comp./step (s) | IO Time (s)\n");
@@ -460,6 +520,37 @@ int main(int argc, char **argv){
   hydro_coreStateLogDump();
   MPI_Barrier(MPI_COMM_WORLD); 
 
+  MPI_Barrier(MPI_COMM_WORLD); 
+  if(mpi_rank_world == 0){
+    printf("\n_____________________#######_________  TOWER-SUMMARY @ it = %d _________#######____________________ \n", it);
+    fflush(stdout);
+  } //if mpi_rank_world  
+  MPI_Barrier(MPI_COMM_WORLD); 
+#ifdef DEBUG_TOWER
+  for(int mrank=0; mrank < mpi_size_world; mrank++){
+     MPI_Barrier(MPI_COMM_WORLD); 
+     if(mrank == mpi_rank_world){
+       for(int towerCount = 0; towerCount < rank_nTowers; towerCount++){
+          printf("=========================================  TOWER-ID %d =========================================== \n", towerIDs[towerCount]);
+          for(int k=0; k < Nz; k++){
+	    printf("%d: ",k);
+            for(int towfld=0; towfld < 15; towfld++){
+	     printf("%f, ",towersData[(NtBatch)*towerCount*towerInstanceSize + (NtBatch-1)*towerInstanceSize+towfld*Nz+k]);
+	    }
+	    printf("\n");
+	  }
+          printf("************ surface values ****************\n");
+          for(int surfld=0; surfld < 6; surfld++){
+	     printf("%f, ",towersSurfData[(NtBatch)*towerCount*towerSurfInstanceSize + (NtBatch-1)*towerSurfInstanceSize+surfld]);
+	  }
+	  printf("\n");
+       }
+       fflush(stdout);
+     }
+     MPI_Barrier(MPI_COMM_WORLD); 
+  }//end for mrank	  
+  MPI_Barrier(MPI_COMM_WORLD); 
+#endif
   if(mpi_rank_world == 0){
     printf("Dumping state at timestep = %d...\n",it);
     fflush(stdout);
@@ -481,6 +572,11 @@ int main(int argc, char **argv){
 #endif
   }
 #endif
+  MPI_Barrier(MPI_COMM_WORLD); 
+  //Dump tower data if appropriate
+  if(towerIOSelector > 0){
+    ioWriteBinaryTowerFileSingleBatch(it, NtBatch, Nz, simTimeBatch, towersData, towersSurfData, towerIDs, rank_nTowers, towerInstanceSize, towerSurfInstanceSize);
+  }
   MPI_Barrier(MPI_COMM_WORLD); 
   mpi_t4 = MPI_Wtime();    //Mark the walltime to measure IO duration
   mpi_t2 = MPI_Wtime();    //Mark the walltime to measure final timestep summary and performance.
@@ -492,7 +588,6 @@ int main(int argc, char **argv){
     printf("   %8.4f \t| %8d \t|  %8.4f \t|  %8.4f \t |  %9.6f \n", (mpi_t2-mpi_t1), 0, 
             (mpi_t2-mpi_t1)/NtBatch, (mpi_t2-mpi_t1-(mpi_t4-mpi_t3))/NtBatch, (mpi_t4-mpi_t3));
     printf("\n********************************************************************************************************\n");
-    fflush(stdout);
     printf("Your FastEddy simulation is complete!\n");
     printf("Cleaning up...\n");
     fflush(stdout);
@@ -532,6 +627,7 @@ int main(int argc, char **argv){
   /* Finalize the FEMPI environment */
   if(mpi_rank_world == 0){
     printf("Shutting down MPI...\n Goodbye!\n");
+    fflush(stdout);
   } //if mpi_rank_world == 0
   MPI_Barrier(MPI_COMM_WORLD);
   errorCode = fempi_FinalizeMPI();

SRC/FEMAIN/Makefile

@@ -36,7 +36,7 @@ DEBUG_CFLAGS = -g
 DEFINES = -DCUB_IGNORE_DEPRECATED_CPP_DIALECT -DTHRUST_IGNORE_DEPRECATED_CPP_DIALECT
 
 TEST_CFLAGS = -Wall -m64 ${DEFINES} ${INCLUDES} ${OTHER_INCLUDES}
-ARCH_CU_FLAGS = -arch=sm_70
+ARCH_CU_FLAGS = -arch=sm_80
 TEST_CU_CFLAGS = ${ARCH_CU_FLAGS} -m64 -std=c++11 ${DEFINES} ${INCLUDES} ${OTHER_INCLUDES}
 
 L_CPPFLAGS =
@@ -189,7 +189,8 @@ all: FastEddy
         ../HYDRO_CORE/CUDA/cuda_largeScaleForcingsDevice.cu \
         ../HYDRO_CORE/CUDA/cuda_moistureDevice.cu \
         ../HYDRO_CORE/CUDA/cuda_filtersDevice.cu \
-	../HYDRO_CORE/CUDA/cuda_cellpertDevice.cu
+	../HYDRO_CORE/CUDA/cuda_cellpertDevice.cu \
+	../HYDRO_CORE/CUDA/cuda_towersDevice.cu
 	$(TEST_CU_CC) $(TEST_CU_CFLAGS) -dc $< -o $@
 ################################################################################
 # Generic Executable