diff --git a/SRC/EXTENSIONS/GAD/CUDA/cuda_GADDevice.cu b/SRC/EXTENSIONS/GAD/CUDA/cuda_GADDevice.cu index 7d37b553..a38ca9e7 100644 --- a/SRC/EXTENSIONS/GAD/CUDA/cuda_GADDevice.cu +++ b/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); diff --git a/SRC/EXTENSIONS/GAD/GAD.c b/SRC/EXTENSIONS/GAD/GAD.c index 5909d0c7..e3efe91c 100644 --- a/SRC/EXTENSIONS/GAD/GAD.c +++ b/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... diff --git a/SRC/EXTENSIONS/URBAN/CUDA/cuda_urbanDevice.cu b/SRC/EXTENSIONS/URBAN/CUDA/cuda_urbanDevice.cu index 1d02ee18..10438383 100644 --- a/SRC/EXTENSIONS/URBAN/CUDA/cuda_urbanDevice.cu +++ b/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){ diff --git a/SRC/FEMAIN/FastEddy.c b/SRC/FEMAIN/FastEddy.c index 51d85ab2..e716dd7a 100644 --- a/SRC/FEMAIN/FastEddy.c +++ b/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(); diff --git a/SRC/FEMAIN/Makefile b/SRC/FEMAIN/Makefile index b04ff07d..9ac162df 100644 --- a/SRC/FEMAIN/Makefile +++ b/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 diff --git a/SRC/GRID/CUDA/cuda_gridDevice.cu b/SRC/GRID/CUDA/cuda_gridDevice.cu index c0abd8d5..96f83d55 100644 --- a/SRC/GRID/CUDA/cuda_gridDevice.cu +++ b/SRC/GRID/CUDA/cuda_gridDevice.cu @@ -58,6 +58,9 @@ float *J33_d; // dz/d_zeta float *D_Jac_d; //Determinant of the Jacbian (called scale factor i.e. if d_xi=d_eta=d_zeta=1, then cell volume) float *invD_Jac_d; //inverse Determinant of the Jacbian +float* lat_d; /* latitude in degrees north "()" 2-d array (x by y) (m)*/ +float* lon_d; /* longitude in degrees east "()" 2-d array (x by y) (m)*/ + /*#################------------------- CUDA_GRID module function definitions ---------------------#################*/ /*----->>>>> int cuda_gridDeviceSetup(); ---------------------------------------------------------------------- * Used to cudaMalloc and cudaMemcpy parameters and coordinate arrays, and for the GRID_CUDA module. @@ -65,6 +68,7 @@ float *invD_Jac_d; //inverse Determinant of the Jacbian extern "C" int cuda_gridDeviceSetup(){ int errorCode = CUDA_GRID_SUCCESS; size_t Nelems; + size_t Nelems2d; #ifdef DEBUG cudaEvent_t startE, stopE; float elapsedTime; @@ -97,7 +101,6 @@ extern "C" int cuda_gridDeviceSetup(){ cudaMemcpyToSymbol(jMax_d, &jMax, sizeof(int)); cudaMemcpyToSymbol(kMin_d, &kMin, sizeof(int)); cudaMemcpyToSymbol(kMax_d, &kMax, sizeof(int)); - gpuErrchk( cudaPeekAtLastError() ); /*Check for errors in the cudaMemCpy calls*/ /*Set the full memory block number of elements for grid fields*/ Nelems = (size_t)((Nxp+2*Nh)*(Nyp+2*Nh)*(Nzp+2*Nh)); @@ -115,7 +118,6 @@ extern "C" int cuda_gridDeviceSetup(){ fecuda_DeviceMalloc(Nelems, &J33_d); fecuda_DeviceMalloc(Nelems, &D_Jac_d); fecuda_DeviceMalloc(Nelems, &invD_Jac_d); - gpuErrchk( cudaPeekAtLastError() ); /*Check for errors in the cudaMalloc calls*/ /* cudaMemcpy the GRID arrays from Host to Device*/ /* Coordinate Arrays */ @@ -130,8 +132,15 @@ extern "C" int cuda_gridDeviceSetup(){ cudaMemcpy(J33_d, J33, Nelems*sizeof(float), cudaMemcpyHostToDevice); cudaMemcpy(D_Jac_d, D_Jac, Nelems*sizeof(float), cudaMemcpyHostToDevice); cudaMemcpy(invD_Jac_d, invD_Jac, Nelems*sizeof(float), cudaMemcpyHostToDevice); + + Nelems2d = (size_t)((Nxp+2*Nh)*(Nyp+2*Nh)); + fecuda_DeviceMalloc(Nelems2d, &lat_d); + cudaMemcpy(lat_d, lat, Nelems2d*sizeof(float), cudaMemcpyHostToDevice); + fecuda_DeviceMalloc(Nelems2d, &lon_d); + cudaMemcpy(lon_d, lon, Nelems2d*sizeof(float), cudaMemcpyHostToDevice); + gpuErrchk( cudaPeekAtLastError() ); /*Check for errors in the cudaMemCpy calls*/ - + #ifdef DEBUG /*Launch an independent GPU calculation of the GRID arrays*/ /*Synchronize the Device*/ @@ -176,21 +185,20 @@ extern "C" int cuda_gridDeviceCleanup(){ /* metric tensor fields */ cudaFree(J13_d); cudaFree(J23_d); - gpuErrchk( cudaPeekAtLastError() ); /*Check for errors in the cudaMemCpy calls*/ cudaFree(J31_d); cudaFree(J32_d); cudaFree(J33_d); - gpuErrchk( cudaPeekAtLastError() ); /*Check for errors in the cudaMemCpy calls*/ cudaFree(D_Jac_d); cudaFree(invD_Jac_d); /* coordinate fields */ cudaFree(xPos_d); cudaFree(yPos_d); - gpuErrchk( cudaPeekAtLastError() ); /*Check for errors in the cudaMemCpy calls*/ cudaFree(zPos_d); cudaFree(topoPos_d); - gpuErrchk( cudaPeekAtLastError() ); /*Check for errors in the cudaMemCpy calls*/ - + + cudaFree(lat_d); + cudaFree(lon_d); + return(errorCode); }//end cuda_gridDeviceCleanup() diff --git a/SRC/GRID/CUDA/cuda_gridDevice_cu.h b/SRC/GRID/CUDA/cuda_gridDevice_cu.h index 8e85ec67..335a2d0f 100644 --- a/SRC/GRID/CUDA/cuda_gridDevice_cu.h +++ b/SRC/GRID/CUDA/cuda_gridDevice_cu.h @@ -55,6 +55,9 @@ extern float *J33_d; // dz/d_zeta extern float *D_Jac_d; //Determinant of the Jacbian (called scale factor i.e. if d_xi=d_eta=d_zeta=1, then cell volume) extern float *invD_Jac_d; //inverse Determinant of the Jacbian +extern float* lat_d; /* latitude in degrees north "()" 2-d array (x by y) (m)*/ +extern float* lon_d; /* longitude in degrees east "()" 2-d array (x by y) (m)*/ + /*#################------------------- GRID_CUDADEV module function declarations ---------------------##############*/ /*----->>>>> int cuda_gridDeviceSetup(); ---------------------------------------------------------------------- diff --git a/SRC/GRID/grid.c b/SRC/GRID/grid.c index d2f983d6..46a4a3a8 100644 --- a/SRC/GRID/grid.c +++ b/SRC/GRID/grid.c @@ -63,13 +63,25 @@ float *topoPos; /*Terrain elevation (z in meters) at the cell center position in float *J13; // dx/d_zeta float *J23; // dy/d_zeta +//float *J11; // dx/d_xi -- assumed = 1.0 +//float *J12; // dx/d_eta -- assumed = 0.0 + +//float *J21; // dy/d_xi -- assumed = 0.0 +//float *J22; // dy/d_eta -- assumed = 1.0 + +float *J13; // dx/d_zeta +float *J23; // dy/d_zeta + float *J31; // dz/d_xi float *J32; // dz/d_eta float *J33; // dz/d_zeta float *D_Jac; //Determinant of the Jacobian (called scale factor i.e. if d_xi=d_eta=d_zeta=1, then cell volume) float *invD_Jac; //inverse Determinant of the Jacobian - + +float* lat; /* latitude in degrees north "()" 2-d array (x by y) (m)*/ +float* lon; /* longitude in degrees east "()" 2-d array (x by y) (m)*/ + /*######################------------------- GRID module function definitions ---------------------#################*/ /*----->>>>> int gridGetParams(); ---------------------------------------------------------------------- @@ -247,7 +259,7 @@ int gridInit(){ yPos = memAllocateFloat3DField(Nxp, Nyp, Nzp, Nh, "yPos"); zPos = memAllocateFloat3DField(Nxp, Nyp, Nzp, Nh, "zPos"); topoPos = memAllocateFloat2DField(Nxp, Nyp, Nh, "topoPos"); - topoPosGlobal = memAllocateFloat2DField(Nx, Ny, 0, "topoPos"); + topoPosGlobal = memAllocateFloat2DField(Nx, Ny, 0, "topoPosGlobal"); /* Metric Tensors Fields */ J13 = memAllocateFloat3DField(Nxp, Nyp, Nzp, Nh, "J13"); J23 = memAllocateFloat3DField(Nxp, Nyp, Nzp, Nh, "J23"); @@ -259,14 +271,11 @@ int gridInit(){ } // end if errorCode indicates no errors thus far /*Register these fields with the IO module*/ - /********* FOR THE MOMENT THESE SHOULD BE STRICTLY GLOBAL DOMAIN VARIABLE FIELDS ********/ if(errorCode == GRID_SUCCESS){ ioerrorCode = ioRegisterVar("xPos", "float", 4, dims4d, xPos); ioerrorCode = ioRegisterVar("yPos", "float", 4, dims4d, yPos); ioerrorCode = ioRegisterVar("zPos", "float", 4, dims4d, zPos); ioerrorCode = ioRegisterVar("topoPos", "float", 3, dims2dTD, topoPos); - printf("gridInit:topoPos stored at %p, has been registered with IO.\n", - &topoPosGlobal); fflush(stdout); if(ioerrorCode!=0){ printf("Error in registering GRID module coordinate fields with IO.\n"); @@ -309,6 +318,15 @@ int gridInit(){ #endif } // end if errorCode indicates no errors thus far + // Allocate 2d arrays of latitude and longitude + lat = memAllocateFloat2DField(Nxp, Nyp, Nh, "lat"); + lon = memAllocateFloat2DField(Nxp, Nyp, Nh, "lon"); + errorCode = ioRegisterVar("lat", "float", 3, dims2dTD, lat); + errorCode = ioRegisterVar("lon", "float", 3, dims2dTD, lon); + // Add NetCDF attributes for the registered variable + ioerrorCode = ioAddStandardAttrs("lat", "degrees", "latitude", NULL); + ioerrorCode = ioAddStandardAttrs("lon", "degrees", "longitude", NULL); + #ifdef DEBUG //#if 1 printf("mpi_rank_world %d/%d: Finished gridInit()!\n",mpi_rank_world,mpi_size_world); @@ -873,6 +891,154 @@ int calculateJacobians(){ return(errorCode); } //end calculateJacobians +/*----->>>>> int gridGetIJindsFromXYPosition(); ------------------------------------------------------------ +* Used to determine the i,j indices of an mpi_rank subdomain +* coordinate frame of the cell that contains the point xLoc,Yloc. +*/ +int gridGetIJindsFromXYPosition(float xLoc, float yLoc, int *iIndx, int *jIndx){ + int errorCode = GRID_SUCCESS; + int i,j,ijk; + *iIndx = -1; + *jIndx = -1; + for(i=iMin; i < iMax; i++){ + for(j=jMin; j < jMax; j++){ + ijk = i*(Nyp+2*Nh)*(Nzp+2*Nh)+j*(Nzp+2*Nh)+kMin; + if( (xLoc > (xPos[ijk]-0.5*d_xi) && xLoc <= (xPos[ijk]+0.5*d_xi)) + && (yLoc > (yPos[ijk]-0.5*d_eta) && yLoc <= (yPos[ijk]+0.5*d_eta)) ){ + *iIndx = i; + *jIndx = j; + } + } + } + if((*iIndx < iMin) || (*jIndx < jMin)){ + printf("Rank %d/%d: gridGetIJindsFromXYPosition(): Tower indices not found in this rank's subdomain!\n", + mpi_rank_world, mpi_size_world); + fflush(stdout); + } + return (errorCode); +} //end gridGetIJindsFromXYPosition() + +/*----->>>>> int gridGetIJindsFromLatLonPosition(); ------------------------------------------------------------ +* Used to determine the i,j indices of an mpi_rank subdomain +* coordinate frame of the cell that contains the point latLoc,lonloc. +*/ +int gridGetIJindsFromLatLonPosition(float lonLoc, float latLoc, int *iIndx, int *jIndx){ + int errorCode = GRID_SUCCESS; + int i,j,ij; + double dr; + double min_dr; + + min_dr = DBL_MAX; + *iIndx = -1; + *jIndx = -1; + for(i=iMin; i < iMax; i++){ + for(j=jMin; j < jMax; j++){ + ij = i*(Nyp+2*Nh)+j; + dr = sqrt( pow((lon[ij]-lonLoc),2.0)+pow((lat[ij]-latLoc),2.0) ); + if( (dr >= 0.0 ) && (dr < min_dr) ){ + min_dr = dr; + *iIndx = i; + *jIndx = j; + } + } + } + if((*iIndx < iMin) || (*jIndx < jMin)){ + printf("Rank %d/%d: gridGetIJindsFromXYPosition(): Tower indices not found in this rank's subdomain!\n", + mpi_rank_world, mpi_size_world); + fflush(stdout); + } + return (errorCode); +} //end gridGetIJindsFromLatLonPosition() + +/*----->>>>> int gridGetRankFromXYPosition(); ------------------------------------------------------------ +* Used to determine the mpi_rank with a subdomain +* that contains the xLoc,Yloc. +*/ +int gridGetRankFromXYPosition(float xLoc, float yLoc){ + int errorCode = GRID_SUCCESS; + int ret_rank; + int tmp_ret_rank = -1; + int ijk_min; + int ijk_max; + ijk_min = iMin*(Nyp+2*Nh)*(Nzp+2*Nh)+jMin*(Nzp+2*Nh)+kMin; + ijk_max = (iMax-1)*(Nyp+2*Nh)*(Nzp+2*Nh)+(jMax-1)*(Nzp+2*Nh)+kMin; + if( (xPos[ijk_min]-0.5*d_xi < xLoc) && (yPos[ijk_min]-0.5*d_eta < yLoc) ){ + if( (xPos[ijk_max]+0.5*d_xi >= xLoc) && (yPos[ijk_max]+0.5*d_eta >= yLoc) ){ + tmp_ret_rank = mpi_rank_world; + } + } + errorCode = MPI_Allreduce(&tmp_ret_rank, &ret_rank, 1, + MPI_INTEGER, MPI_MAX, MPI_COMM_WORLD); + if(errorCode != MPI_SUCCESS){ + printf("Rank %d/%d gridGetRankFromXYPosition(): MPI_Allreducei returned with MPI_ERROR = %d!\n", + mpi_rank_world, mpi_size_world, errorCode); + fflush(stdout); + } + return (ret_rank); +} //end gridGetRankFromXYPosition() + +/*----->>>>> int gridGetRankFromLatLonPosition(); ------------------------------------------------------------ +* Used to determine the mpi_rank with a subdomain +* that contains the latLoc,lonloc. +*/ +int gridGetRankFromLatLonPosition(double lonLoc, double latLoc){ + int errorCode = GRID_SUCCESS; + int ret_rank; + int tmp_ret_rank = -1; + int ij_min; + int ij_max; + double dlat; + double dlon; + + ij_min = iMin*(Nyp+2*Nh)+jMin; + ij_max = (iMax-1)*(Nyp+2*Nh)+(jMax-1); + + dlat = fabs(lat[ij_min + 1]-lat[ij_min]); + dlon = fabs(lon[ij_min + (Nyp+2*Nh)]-lon[ij_min]); + + if( (lon[ij_min]-0.5*dlon < lonLoc) && (lat[ij_min]-0.5*dlat < latLoc) ){ + if( (lon[ij_max]+0.5*dlon >= lonLoc) && (lat[ij_max]+0.5*dlat >= latLoc) ){ + tmp_ret_rank = mpi_rank_world; + } + } + errorCode = MPI_Allreduce(&tmp_ret_rank, &ret_rank, 1, + MPI_INTEGER, MPI_MAX, MPI_COMM_WORLD); + if(errorCode != MPI_SUCCESS){ + printf("Rank %d/%d gridGetRankFromLatLonPosition(): MPI_Allreduce returned with MPI_ERROR = %d!\n", + mpi_rank_world, mpi_size_world, errorCode); + fflush(stdout); + } + return (ret_rank); +} //end gridGetRankFromLatLonPosition() + +/*----->>>>> int gridGetXYOffsetsFromXYPosition(); ------------------------------------------------------------ +* Used to determine the x,y position offsets from a predetermined i,j-index cell center x,y coordinate +*/ +int gridGetXYOffsetsFromCellIndices(float xLoc, float yLoc, int iIndx, int jIndx, float *xOff, float *yOff){ + int errorCode = GRID_SUCCESS; + int ijk; + + ijk = iIndx*(Nyp+2*Nh)*(Nzp+2*Nh)+jIndx*(Nzp+2*Nh)+kMin; + *xOff = xPos[ijk]-xLoc; + *yOff = yPos[ijk]-yLoc; + + return (errorCode); +} //end gridGetXYOffsetsFromXYPosition() + +/*----->>>>> int gridGetLatLonOffsetsFromXYPosition(); ------------------------------------------------------------ +* Used to determine the lat,lon position offsets from a predetermined i,j-index cell center lat,lon coordinate +*/ +int gridGetLatLonOffsetsFromCellIndices(double lonLoc, double latLoc, int iIndx, int jIndx, double *lonOff, double *latOff){ + int errorCode = GRID_SUCCESS; + int ij; + + ij = iIndx*(Nyp+2*Nh)+jIndx; + *lonOff = lon[ij]-lonLoc; + *latOff = lat[ij]-latLoc; + + return (errorCode); +} //end gridGetLatLonOffsetsFromCellIndices() + /*----->>>>> int singleRankGridHaloInit(); ------------------------------------------------------------ * Used to setup xPos,yPos,zPos halos on all x-y boundaries * when under single-rank setup (i.e. mpi_size_world ==1). diff --git a/SRC/GRID/grid.h b/SRC/GRID/grid.h index 90f6bba9..0096ded0 100644 --- a/SRC/GRID/grid.h +++ b/SRC/GRID/grid.h @@ -58,6 +58,15 @@ extern float *topoPosGlobal; /*Terrain elevation (z in meters) at the cell cente extern float *J13; // dx/d_zeta extern float *J23; // dy/d_zeta +//extern float *J11; // dx/d_xi -- assumed = 1.0 +//extern float *J12; // dx/d_eta -- assumed = 0.0 + +//extern float *J21; // dy/d_xi -- assumed = 0.0 +//extern float *J22; // dy/d_eta -- assumed = 1.0 + +extern float *J13; // dx/d_zeta +extern float *J23; // dy/d_zeta + extern float *J31; // dz/d_xi extern float *J32; // dz/d_eta extern float *J33; // dz/d_zeta @@ -65,6 +74,9 @@ extern float *J33; // dz/d_zeta extern float *D_Jac; //Determinant of the Jacbian (called scale factor i.e. if d_xi=d_eta=d_zeta=1, then cell volume) extern float *invD_Jac; //inverse Determinant of the Jacbian +extern float* lat; /* latitude in degrees north "()" 2-d array (x by y) (m)*/ +extern float* lon; /* longitude in degrees east "()" 2-d array (x by y) (m)*/ + /*######################------------------- GRID module function declarations ---------------------#################*/ /*----->>>>> int gridGetParams(); ---------------------------------------------------------------------- @@ -88,6 +100,40 @@ int gridSecondaryPreparations(); */ int calculateJacobians(); +/*----->>>>> int gridGetRankFromXYPosition(); ------------------------------------------------------------ +* Used to determine the mpi_rank with a subdomain +* that contains the xLoc,Yloc. +*/ +int gridGetIJindsFromXYPosition(float xLoc, float yLoc, int *iIndx, int *jIndx); + +/*----->>>>> int gridGetRankFromXYPosition(); ------------------------------------------------------------ +* Used to determine the i,j indices of an mpi_rank subdomain +* coordinate frame of the cell that contains the point latLoc,lonloc. +*/ +int gridGetIJindsFromLatLonPosition(float lonLoc, float latLoc, int *iIndx, int *jIndx); + +/*----->>>>> int gridGetRankFromXYPosition(); ------------------------------------------------------------ +* Used to determine the mpi_rank with a subdomain +* that contains the xLoc,Yloc. +*/ +int gridGetRankFromXYPosition(float xLoc, float yLoc); + +/*----->>>>> int gridGetRankFromLatLonPosition(); ------------------------------------------------------------ +* Used to determine the mpi_rank with a subdomain +* that contains the latLoc,lonloc. +*/ +int gridGetRankFromLatLonPosition(double lonLoc, double latLoc); + +/*----->>>>> int gridGetXYOffsetsFromXYPosition(); ------------------------------------------------------------ +* Used to determine the x,y position offsets from a predetermined i,j-index cell center x,y coordinate +*/ +int gridGetXYOffsetsFromCellIndices(float xLoc, float yLoc, int iIndx, int jIndx, float *xOff, float *yOff); + +/*----->>>>> int gridGetLatLonOffsetsFromXYPosition(); ------------------------------------------------------------ +* Used to determine the lat,lon position offsets from a predetermined i,j-index cell center lat,lon coordinate +*/ +int gridGetLatLonOffsetsFromCellIndices(double lonLoc, double latLoc, int iIndx, int jIndx, double *lonOff, double *latOff); + /*----->>>>> int singleRankGridHaloInit(); ------------------------------------------------------------ * Used to setup xPos,yPos,zPos halos on all x-y boundaries * when under single-rank setup (i.e. mpi_size_world ==1). diff --git a/SRC/HYDRO_CORE/CUDA/cuda_BCsDevice.cu b/SRC/HYDRO_CORE/CUDA/cuda_BCsDevice.cu index 1749cd5c..fb1aac43 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_BCsDevice.cu +++ b/SRC/HYDRO_CORE/CUDA/cuda_BCsDevice.cu @@ -149,22 +149,6 @@ extern "C" int cuda_BCsDeviceCleanup(){ }//end cuda_moistureDeviceCleanup() -/*----->>>>> int cuda_hydroCoreDeviceSecondaryStageSetup(); --------------------------------------------------------- -* Secondary initializations at the device level for BCs -*/ -extern "C" int cuda_hydroCoreDeviceSecondaryStageSetup(float dt){ - int errorCode = CUDA_HYDRO_CORE_SUCCESS; - int BdyUpdateSteps; - - /*Compute the number of timesteps between BndyPlane Updates*/ - BdyUpdateSteps = (int) roundf(dtBdyPlaneBCs/dt); - cudaMemcpyToSymbol(BdyUpdateSteps_d, &BdyUpdateSteps, sizeof(int)); - - printf("%d/%d cuda_hydroCoreDeviceSecondaryStageSetup(): BdyUpdateSteps = %d \n",mpi_rank_world,mpi_size_world,BdyUpdateSteps); - fflush(stdout); - return(errorCode); -} - /*----->>>>> int cuda_hydroCoreDeviceBdyPlanesUpdate(); ----------------------------------------------------------------- * Utility to cycle device-sided pointers and push (copy) newest BndyPlanes from Host to Device */ diff --git a/SRC/HYDRO_CORE/CUDA/cuda_BCsDevice_cu.h b/SRC/HYDRO_CORE/CUDA/cuda_BCsDevice_cu.h index 9dca3dd8..e879d27b 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_BCsDevice_cu.h +++ b/SRC/HYDRO_CORE/CUDA/cuda_BCsDevice_cu.h @@ -59,11 +59,6 @@ extern "C" int cuda_BCsDeviceSetup(); */ extern "C" int cuda_BCsDeviceCleanup(); -/*----->>>>> int cuda_hydroCoreDeviceSecondaryStageSetup(float dt); ----------------------------------------------------------------- -* Secondary initializations at the device level for BCs -*/ -extern "C" int cuda_hydroCoreDeviceSecondaryStageSetup(float dt); - /*----->>>>> int cuda_hydroCoreDeviceBdyPlanesUpdate(); ----------------------------------------------------------------- * Utility to cycle device-sided pointers and push (copy) newest BndyPlanes from Host to Device */ diff --git a/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice.cu b/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice.cu index c869c9f6..ee48fcdf 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice.cu +++ b/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice.cu @@ -161,6 +161,87 @@ __device__ void cudaDevice_UpstreamDivAdvFlux(float* scalarField, float* scalarF } //end cudaDevice_UpstreamDivAdvFlux( +/*----->>>>> __device__ void cudaDevice_UpstreamDivAdvFluxX(); -------------------------------------------------- +* This is the cuda version of the UpstreamDivAdvFluxX routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_UpstreamDivAdvFluxX(float* scalarField, float* scalarFadv, float* u_cf, float* invD_Jac_d){ + int i,j,k; + int ijk,im1jk,ip1jk; + int iStride,jStride,kStride; + float DscalarDx; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + im1jk = (i-1)*iStride + j*jStride + k*kStride; + ip1jk = (i+1)*iStride + j*jStride + k*kStride; + DscalarDx = ( ( fmaxf(0.0,u_cf[ ip1jk ])*scalarField[ ijk ]+fminf(0.0,u_cf[ ip1jk ])*scalarField[ ip1jk ]) + -( fmaxf(0.0,u_cf[ ijk ])*scalarField[ im1jk ]+fminf(0.0,u_cf[ ijk ])*scalarField[ ijk ]) ); + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDx; + +} //end cudaDevice_UpstreamDivAdvFluxX( + +/*----->>>>> __device__ void cudaDevice_UpstreamDivAdvFluxY(); -------------------------------------------------- +* This is the cuda version of the UpstreamDivAdvFluxY routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_UpstreamDivAdvFluxY(float* scalarField, float* scalarFadv, float* v_cf, float* invD_Jac_d){ + int i,j,k; + int ijk,ijm1k,ijp1k; + int iStride,jStride,kStride; + float DscalarDy; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + ijm1k = i*iStride + (j-1)*jStride + k*kStride; + ijp1k = i*iStride + (j+1)*jStride + k*kStride; + DscalarDy = ( ( fmaxf(0.0,v_cf[ ijp1k ])*scalarField[ ijk ]+fminf(0.0,v_cf[ ijp1k ])*scalarField[ ijp1k ]) + -( fmaxf(0.0,v_cf[ ijk ])*scalarField[ ijm1k ]+fminf(0.0,v_cf[ ijk ])*scalarField[ ijk ]) ); + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDy; + +} //end cudaDevice_UpstreamDivAdvFluxY( + +/*----->>>>> __device__ void cudaDevice_UpstreamDivAdvFluxZ(); -------------------------------------------------- +* This is the cuda version of the UpstreamDivAdvFluxZ routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_UpstreamDivAdvFluxZ(float* scalarField, float* scalarFadv, float* w_cf, float* invD_Jac_d){ + int i,j,k; + int ijk,ijkm1,ijkp1; + int iStride,jStride,kStride; + float DscalarDz; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + ijkm1 = i*iStride + j*jStride + (k-1)*kStride; + ijkp1 = i*iStride + j*jStride + (k+1)*kStride; + DscalarDz = ( ( fmaxf(0.0,w_cf[ ijkp1 ])*scalarField[ ijk ]+fminf(0.0,w_cf[ ijkp1 ])*scalarField[ ijkp1 ]) + -( fmaxf(0.0,w_cf[ ijk ])*scalarField[ ijkm1 ]+fminf(0.0,w_cf[ ijk ])*scalarField[ ijk ]) ); + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDz; + +} //end cudaDevice_UpstreamDivAdvFluxZ( + /*----->>>>> __device__ void cudaDevice_SecondDivAdvFlux(); -----------------------------------------------*/ __device__ void cudaDevice_SecondDivAdvFlux(float* scalarField, float* scalarFadv, float* u_cf, float* v_cf, float* w_cf, float* invD_Jac_d){ @@ -315,6 +396,117 @@ __device__ void cudaDevice_HYB34DivAdvFlux(float* scalarField, float* scalarFadv } //end cudaDevice_HYB34DivAdvFlux() +/*----->>>>> __device__ void cudaDevice_HYB34DivAdvFluxX(); -------------------------------------------------- +* This is the cuda version of the hydro_coreHYB34DivAdvFluxX routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_HYB34DivAdvFluxX(float* scalarField, float* scalarFadv, float* u_cf, float b_hyb_p, float* invD_Jac_d){ + + int i,j,k; + int ijk,im1jk,ip1jk,im2jk,ip2jk; + int iStride,jStride,kStride; + float DscalarDx; + float one_twelfth; + float flxx_ipf,flxx_imf; + + one_twelfth = 1.0/12.0; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + im1jk = (i-1)*iStride + j*jStride + k*kStride; + ip1jk = (i+1)*iStride + j*jStride + k*kStride; + im2jk = (i-2)*iStride + j*jStride + k*kStride; + ip2jk = (i+2)*iStride + j*jStride + k*kStride; + + flxx_ipf = one_twelfth * ( 7.0*(scalarField[ ip1jk ]+scalarField[ ijk ])-(scalarField[ ip2jk ]+scalarField[ im1jk ])+ + (1.0-b_hyb_p)*copysign(1.0,u_cf[ ip1jk ])*((scalarField[ ip2jk ]-scalarField[ im1jk ])-3.0*(scalarField[ ip1jk ]-scalarField[ ijk ])) ); + flxx_imf = one_twelfth * ( 7.0*(scalarField[ ijk ]+scalarField[ im1jk ])-(scalarField[ ip1jk ]+scalarField[ im2jk ])+ + (1.0-b_hyb_p)*copysign(1.0,u_cf[ ijk ])*((scalarField[ ip1jk ]-scalarField[ im2jk ])-3.0*(scalarField[ ijk ]-scalarField[ im1jk ])) ); + DscalarDx = u_cf[ ip1jk ]*flxx_ipf - u_cf[ ijk ]*flxx_imf; + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDx; + +} //end cudaDevice_HYB34DivAdvFluxX() + +/*----->>>>> __device__ void cudaDevice_HYB34DivAdvFluxY(); -------------------------------------------------- +* This is the cuda version of the hydro_coreHYB34DivAdvFluxY routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_HYB34DivAdvFluxY(float* scalarField, float* scalarFadv, float* v_cf, float b_hyb_p, float* invD_Jac_d){ + + int i,j,k; + int ijk,ijm1k,ijp1k,ijm2k,ijp2k; + int iStride,jStride,kStride; + float DscalarDy; + float one_twelfth; + float flxy_jpf,flxy_jmf; + + one_twelfth = 1.0/12.0; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + ijm1k = i*iStride + (j-1)*jStride + k*kStride; + ijp1k = i*iStride + (j+1)*jStride + k*kStride; + ijm2k = i*iStride + (j-2)*jStride + k*kStride; + ijp2k = i*iStride + (j+2)*jStride + k*kStride; + + flxy_jpf = one_twelfth * ( 7.0*(scalarField[ ijp1k ]+scalarField[ ijk ])-(scalarField[ ijp2k ]+scalarField[ ijm1k ])+ + (1.0-b_hyb_p)*copysign(1.0,v_cf[ ijp1k ])*((scalarField[ ijp2k ]-scalarField[ ijm1k ])-3.0*(scalarField[ ijp1k ]-scalarField[ ijk ])) ); + flxy_jmf = one_twelfth * ( 7.0*(scalarField[ ijk ]+scalarField[ ijm1k ])-(scalarField[ ijp1k ]+scalarField[ ijm2k ])+ + (1.0-b_hyb_p)*copysign(1.0,v_cf[ ijk ])*((scalarField[ ijp1k ]-scalarField[ ijm2k ])-3.0*(scalarField[ ijk ]-scalarField[ ijm1k ])) ); + DscalarDy = v_cf[ ijp1k ]*flxy_jpf - v_cf[ ijk ]*flxy_jmf; + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDy; + +} //end cudaDevice_HYB34DivAdvFluxY() + +/*----->>>>> __device__ void cudaDevice_HYB34DivAdvFluxZ(); -------------------------------------------------- +* This is the cuda version of the hydro_coreHYB34DivAdvFluxZ routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_HYB34DivAdvFluxZ(float* scalarField, float* scalarFadv, float* w_cf, float b_hyb_p, float* invD_Jac_d){ + + int i,j,k; + int ijk,ijkm1,ijkp1,ijkm2,ijkp2; + int iStride,jStride,kStride; + float DscalarDz; + float one_twelfth; + float flxz_kpf,flxz_kmf; + + one_twelfth = 1.0/12.0; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + ijkm1 = i*iStride + j*jStride + (k-1)*kStride; + ijkp1 = i*iStride + j*jStride + (k+1)*kStride; + ijkm2 = i*iStride + j*jStride + (k-2)*kStride; + ijkp2 = i*iStride + j*jStride + (k+2)*kStride; + + flxz_kpf = one_twelfth * ( 7.0*(scalarField[ ijkp1 ]+scalarField[ ijk ])-(scalarField[ ijkp2 ]+scalarField[ ijkm1 ])+ + (1.0-b_hyb_p)*copysign(1.0,w_cf[ ijkp1 ])*((scalarField[ ijkp2 ]-scalarField[ ijkm1 ])-3.0*(scalarField[ ijkp1 ]-scalarField[ ijk ])) ); + flxz_kmf = one_twelfth * ( 7.0*(scalarField[ ijk ]+scalarField[ ijkm1 ])-(scalarField[ ijkp1 ]+scalarField[ ijkm2 ])+ + (1.0-b_hyb_p)*copysign(1.0,w_cf[ ijk ])*((scalarField[ ijkp1 ]-scalarField[ ijkm2 ])-3.0*(scalarField[ ijk ]-scalarField[ ijkm1 ])) ); + DscalarDz = w_cf[ ijkp1 ]*flxz_kpf - w_cf[ ijk ]*flxz_kmf; + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDz; + +} //end cudaDevice_HYB34DivAdvFluxZ() + /*----->>>>> __device__ void cudaDevice_HYB56DivAdvFlux(); -------------------------------------------------- * This is the cuda version of the hydro_coreHYB56DivAdvFlx routine from the HYDRO_CORE module */ @@ -382,6 +574,132 @@ __device__ void cudaDevice_HYB56DivAdvFlux(float* scalarField, float* scalarFadv } //end cudaDevice_HYB56DivAdvFlux( +/*----->>>>> __device__ void cudaDevice_HYB56DivAdvFluxX(); -------------------------------------------------- +* This is the cuda version of the hydro_coreHYB56DivAdvFluxX routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_HYB56DivAdvFluxX(float* scalarField, float* scalarFadv, float* u_cf, float b_hyb_p, float* invD_Jac_d){ + + int i,j,k; + int ijk,im1jk,ip1jk; + int im2jk,ip2jk; + int im3jk,ip3jk; + int iStride,jStride,kStride; + float DscalarDx; + float one_sixtieth; + float flxx_ipf,flxx_imf; + + one_sixtieth = 1.0/60.0; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + im1jk = (i-1)*iStride + j*jStride + k*kStride; + ip1jk = (i+1)*iStride + j*jStride + k*kStride; + im2jk = (i-2)*iStride + j*jStride + k*kStride; + ip2jk = (i+2)*iStride + j*jStride + k*kStride; + im3jk = (i-3)*iStride + j*jStride + k*kStride; + ip3jk = (i+3)*iStride + j*jStride + k*kStride; + + flxx_ipf = one_sixtieth * ( 37.0*(scalarField[ ip1jk ]+scalarField[ ijk ])-8.0*(scalarField[ ip2jk ]+scalarField[ im1jk ])+(scalarField[ ip3jk ]+scalarField[ im2jk ])- + (1.0-b_hyb_p)*copysign(1.0,u_cf[ ip1jk ])*((scalarField[ ip3jk ]-scalarField[ im2jk ])-5.0*(scalarField[ ip2jk ]-scalarField[ im1jk ])+10.0*(scalarField[ ip1jk ]-scalarField[ ijk ])) ); + flxx_imf = one_sixtieth * ( 37.0*(scalarField[ ijk ]+scalarField[ im1jk ])-8.0*(scalarField[ ip1jk ]+scalarField[ im2jk ])+(scalarField[ ip2jk ]+scalarField[ im3jk ])- + (1.0-b_hyb_p)*copysign(1.0,u_cf[ ijk ])*((scalarField[ ip2jk ]-scalarField[ im3jk ])-5.0*(scalarField[ ip1jk ]-scalarField[ im2jk ])+10.0*(scalarField[ ijk ]-scalarField[ im1jk ])) ); + + DscalarDx = u_cf[ ip1jk ]*flxx_ipf - u_cf[ ijk ]*flxx_imf; + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDx; + +} //end cudaDevice_HYB56DivAdvFluxX( + +/*----->>>>> __device__ void cudaDevice_HYB56DivAdvFluxY(); -------------------------------------------------- +* This is the cuda version of the hydro_coreHYB56DivAdvFluxY routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_HYB56DivAdvFluxY(float* scalarField, float* scalarFadv, float* v_cf, float b_hyb_p, float* invD_Jac_d){ + + int i,j,k; + int ijk,ijm1k,ijp1k; + int ijm2k,ijp2k; + int ijm3k,ijp3k; + int iStride,jStride,kStride; + float DscalarDy; + float one_sixtieth; + float flxy_jpf,flxy_jmf; + + one_sixtieth = 1.0/60.0; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + ijm1k = i*iStride + (j-1)*jStride + k*kStride; + ijp1k = i*iStride + (j+1)*jStride + k*kStride; + ijm2k = i*iStride + (j-2)*jStride + k*kStride; + ijp2k = i*iStride + (j+2)*jStride + k*kStride; + ijm3k = i*iStride + (j-3)*jStride + k*kStride; + ijp3k = i*iStride + (j+3)*jStride + k*kStride; + + flxy_jpf = one_sixtieth * ( 37.0*(scalarField[ ijp1k ]+scalarField[ ijk ])-8.0*(scalarField[ ijp2k ]+scalarField[ ijm1k ])+(scalarField[ ijp3k ]+scalarField[ ijm2k ])- + (1.0-b_hyb_p)*copysign(1.0,v_cf[ ijp1k ])*((scalarField[ ijp3k ]-scalarField[ ijm2k ])-5.0*(scalarField[ ijp2k ]-scalarField[ ijm1k ])+10.0*(scalarField[ ijp1k ]-scalarField[ ijk ])) ); + flxy_jmf = one_sixtieth * ( 37.0*(scalarField[ ijk ]+scalarField[ ijm1k ])-8.0*(scalarField[ ijp1k ]+scalarField[ ijm2k ])+(scalarField[ ijp2k ]+scalarField[ ijm3k ])- + (1.0-b_hyb_p)*copysign(1.0,v_cf[ ijk ])*((scalarField[ ijp2k ]-scalarField[ ijm3k ])-5.0*(scalarField[ ijp1k ]-scalarField[ ijm2k ])+10.0*(scalarField[ ijk ]-scalarField[ ijm1k ])) ); + + DscalarDy = v_cf[ ijp1k ]*flxy_jpf - v_cf[ ijk ]*flxy_jmf; + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDy; + +} //end cudaDevice_HYB56DivAdvFluxY( + +/*----->>>>> __device__ void cudaDevice_HYB56DivAdvFluxZ(); -------------------------------------------------- +* This is the cuda version of the hydro_coreHYB56DivAdvFluxZ routine from the HYDRO_CORE module +*/ +__device__ void cudaDevice_HYB56DivAdvFluxZ(float* scalarField, float* scalarFadv, float* w_cf, float b_hyb_p, float* invD_Jac_d){ + + int i,j,k; + int ijk,ijkm1,ijkp1; + int ijkm2,ijkp2; + int ijkm3,ijkp3; + int iStride,jStride,kStride; + float DscalarDz; + float one_sixtieth; + float flxz_kpf,flxz_kmf; + + one_sixtieth = 1.0/60.0; + + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + ijk = i*iStride + j*jStride + k*kStride; + ijkm1 = i*iStride + j*jStride + (k-1)*kStride; + ijkp1 = i*iStride + j*jStride + (k+1)*kStride; + ijkm2 = i*iStride + j*jStride + (k-2)*kStride; + ijkp2 = i*iStride + j*jStride + (k+2)*kStride; + ijkm3 = i*iStride + j*jStride + (k-3)*kStride; + ijkp3 = i*iStride + j*jStride + (k+3)*kStride; + + flxz_kpf = one_sixtieth * ( 37.0*(scalarField[ ijkp1 ]+scalarField[ ijk ])-8.0*(scalarField[ ijkp2 ]+scalarField[ ijkm1 ])+(scalarField[ ijkp3 ]+scalarField[ ijkm2 ])- + (1.0-b_hyb_p)*copysign(1.0,w_cf[ ijkp1 ])*((scalarField[ ijkp3 ]-scalarField[ ijkm2 ])-5.0*(scalarField[ ijkp2 ]-scalarField[ ijkm1 ])+10.0*(scalarField[ ijkp1 ]-scalarField[ ijk ])) ); + flxz_kmf = one_sixtieth * ( 37.0*(scalarField[ ijk ]+scalarField[ ijkm1 ])-8.0*(scalarField[ ijkp1 ]+scalarField[ ijkm2 ])+(scalarField[ ijkp2 ]+scalarField[ ijkm3 ])- + (1.0-b_hyb_p)*copysign(1.0,w_cf[ ijk ])*((scalarField[ ijkp2 ]-scalarField[ ijkm3 ])-5.0*(scalarField[ ijkp1 ]-scalarField[ ijkm2 ])+10.0*(scalarField[ ijk ]-scalarField[ ijkm1 ])) ); + + DscalarDz = w_cf[ ijkp1 ]*flxz_kpf - w_cf[ ijk ]*flxz_kmf; + scalarFadv[ijk] = scalarFadv[ijk] -invD_Jac_d[ijk]*DscalarDz; + +} //end cudaDevice_HYB56DivAdvFluxZ( + /*----->>>>> __device__ void cudaDevice_WENO3DivAdvFluxX(); -----------------------------------------------*/ __device__ void cudaDevice_WENO3DivAdvFluxX(float* scalarField, float* scalarFadv,float* u_cf, float* invD_Jac_d){ diff --git a/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice_cu.h b/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice_cu.h index f9234949..359ab3a9 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice_cu.h +++ b/SRC/HYDRO_CORE/CUDA/cuda_advectionDevice_cu.h @@ -50,6 +50,15 @@ __device__ void cudaDevice_calcFaceVelocities(float* hydroFlds_d, float* hydroFa __device__ void cudaDevice_UpstreamDivAdvFlux(float* scalarField, float* scalarFadv, float* u_cf, float* v_cf, float* w_cf, float* invD_Jac_d); +/*----->>>>> __device__ void cudaDevice_UpstreamDivAdvFluxX(); -------------------------------------------------- */ +__device__ void cudaDevice_UpstreamDivAdvFluxX(float* scalarField, float* scalarFadv, float* u_cf, float* invD_Jac_d); + +/*----->>>>> __device__ void cudaDevice_UpstreamDivAdvFluxY(); -------------------------------------------------- */ +__device__ void cudaDevice_UpstreamDivAdvFluxY(float* scalarField, float* scalarFadv, float* v_cf, float* invD_Jac_d); + +/*----->>>>> __device__ void cudaDevice_UpstreamDivAdvFluxZ(); -------------------------------------------------- */ +__device__ void cudaDevice_UpstreamDivAdvFluxZ(float* scalarField, float* scalarFadv, float* w_cf, float* invD_Jac_d); + /*----->>>>> __device__ void cudaDevice_SecondDivAdvFlux(); -------------------------------------------------- */ __device__ void cudaDevice_SecondDivAdvFlux(float* scalarField, float* scalarFadv, @@ -67,12 +76,30 @@ __device__ void cudaDevice_QUICKDivAdvFlux(float* scalarField, float* scalarFadv __device__ void cudaDevice_HYB34DivAdvFlux(float* scalarField, float* scalarFadv, float* u_cf, float* v_cf, float* w_cf, float b_hyb_p, float* invD_Jac_d); +/*----->>>>> __device__ void cudaDevice_HYB34DivAdvFluxX(); -------------------------------------------------- */ +__device__ void cudaDevice_HYB34DivAdvFluxX(float* scalarField, float* scalarFadv, float* u_cf, float b_hyb_p, float* invD_Jac_d); + +/*----->>>>> __device__ void cudaDevice_HYB34DivAdvFluxY(); -------------------------------------------------- */ +__device__ void cudaDevice_HYB34DivAdvFluxY(float* scalarField, float* scalarFadv, float* v_cf, float b_hyb_p, float* invD_Jac_d); + +/*----->>>>> __device__ void cudaDevice_HYB34DivAdvFluxZ(); -------------------------------------------------- */ +__device__ void cudaDevice_HYB34DivAdvFluxZ(float* scalarField, float* scalarFadv, float* w_cf, float b_hyb_p, float* invD_Jac_d); + /*----->>>>> __device__ void cudaDevice_HYB56DivAdvFlux(); -------------------------------------------------- * This is the cuda version of the HYB56DivAdvFlux routine from the HYDRO_CORE module */ __device__ void cudaDevice_HYB56DivAdvFlux(float* scalarField, float* scalarFadv, float* u_cf, float* v_cf, float* w_cf, float b_hyb_p, float* invD_Jac_d); +/*----->>>>> __device__ void cudaDevice_HYB56DivAdvFluxX(); -------------------------------------------------- */ +__device__ void cudaDevice_HYB56DivAdvFluxX(float* scalarField, float* scalarFadv, float* u_cf, float b_hyb_p, float* invD_Jac_d); + +/*----->>>>> __device__ void cudaDevice_HYB56DivAdvFluxY(); -------------------------------------------------- */ +__device__ void cudaDevice_HYB56DivAdvFluxY(float* scalarField, float* scalarFadv, float* v_cf, float b_hyb_p, float* invD_Jac_d); + +/*----->>>>> __device__ void cudaDevice_HYB56DivAdvFluxZ(); -------------------------------------------------- */ +__device__ void cudaDevice_HYB56DivAdvFluxZ(float* scalarField, float* scalarFadv, float* w_cf, float b_hyb_p, float* invD_Jac_d); + /*----->>>>> __device__ void cudaDevice_WENO3DivAdvFluxX(); -------------------------------------------------- */ __device__ void cudaDevice_WENO3DivAdvFluxX(float* scalarField, float* scalarFadv,float* u_cf, float* invD_Jac_d); diff --git a/SRC/HYDRO_CORE/CUDA/cuda_cellpertDevice.cu b/SRC/HYDRO_CORE/CUDA/cuda_cellpertDevice.cu index c2133cb8..a0c47fdb 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_cellpertDevice.cu +++ b/SRC/HYDRO_CORE/CUDA/cuda_cellpertDevice.cu @@ -46,7 +46,7 @@ extern "C" int cuda_cellpertDeviceSetup(){ cudaMemcpyToSymbol(cellpert_ktop_d, &cellpert_ktop, sizeof(int)); Nelems1d_xy = (Nx/cellpert_gppc+min(Nx%cellpert_gppc,1))*(2*cellpert_ndbc+min(Ny%cellpert_gppc,1)) + (Ny/cellpert_gppc-2*cellpert_ndbc)*(2*cellpert_ndbc+min(Nx%cellpert_gppc,1)); - Nelems1d = (size_t)(Nelems1d_xy*(cellpert_ktop-cellpert_kbottom+1)); + Nelems1d = (size_t)(Nelems1d_xy*(Nz-cellpert_kbottom+1)); fecuda_DeviceMalloc(Nelems1d, &randcp_d); return(errorCode); @@ -92,7 +92,7 @@ extern "C" int cuda_hydroCoreDeviceBuildCPmethod(int simTime_it){ // uniform distribution of pseudo-random numbers on randcp_d (1d-array) n_xy = (Nx/cellpert_gppc+min(Nx%cellpert_gppc,1))*(2*cellpert_ndbc+min(Ny%cellpert_gppc,1)) + (Ny/cellpert_gppc-2*cellpert_ndbc)*(2*cellpert_ndbc+min(Nx%cellpert_gppc,1)); - n_tot = n_xy*(cellpert_ktop-cellpert_kbottom+1); + n_tot = n_xy*(Nz-cellpert_kbottom+1); curandCreateGenerator(&gen,CURAND_RNG_PSEUDO_DEFAULT); curandSetPseudoRandomGeneratorSeed(gen,(unsigned long long)simTime_it); @@ -218,7 +218,7 @@ __device__ void cudaDevice_CellPerturbation(int i_ind, int j_ind, int k_ind, int ncx = (Nx_tot/gppc)+ ncx_p; ncy_p = min(Ny_tot%gppc,1); ncy = (Ny_tot/gppc); - ncz = cellpert_ktop_d-cellpert_kbottom_d+1; + ncz = Nz-cellpert_kbottom_d+1; // only domain boundary ring nc_xy = ncx*(2*ndbc+ncy_p) + (ncy-2*ndbc)*(2*ndbc+ncx_p); @@ -313,7 +313,7 @@ __device__ void cudaDevice_CellPerturbationMasked(int i_ind, int j_ind, int k_in ncx = (Nx_tot/gppc)+ ncx_p; ncy_p = min(Ny_tot%gppc,1); ncy = (Ny_tot/gppc); - ncz = cellpert_ktop_d-cellpert_kbottom_d+1; + ncz = Nz-cellpert_kbottom_d+1; // only domain boundary ring nc_xy = ncx*(2*ndbc+ncy_p) + (ncy-2*ndbc)*(2*ndbc+ncx_p); diff --git a/SRC/HYDRO_CORE/CUDA/cuda_coriolisDevice.cu b/SRC/HYDRO_CORE/CUDA/cuda_coriolisDevice.cu index c9b2327d..6701a269 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_coriolisDevice.cu +++ b/SRC/HYDRO_CORE/CUDA/cuda_coriolisDevice.cu @@ -26,6 +26,7 @@ __constant__ float corioLS_fact_d; /*large-scale forcing factor on Co */ extern "C" int cuda_coriolisDeviceSetup(){ int errorCode = CUDA_CORIOLIS_SUCCESS; + cudaMemcpyToSymbol(coriolisSelector_d, &coriolisSelector, sizeof(int)); cudaMemcpyToSymbol(corioConstHorz_d, &corioConstHorz, sizeof(float)); cudaMemcpyToSymbol(corioConstVert_d, &corioConstVert, sizeof(float)); @@ -52,10 +53,17 @@ extern "C" int cuda_coriolisDeviceCleanup(){ */ __device__ void cudaDevice_calcCoriolis(float* Frhs_u, float* Frhs_v, float* Frhs_w, float* rho, float* uMom, float* vMom, float* wMom, - float* rhoBS, float* uBS, float* vBS, float* wBS){ + float* rhoBS, float* uBS, float* vBS, float* wBS, + float* lat){ + float pi = acosf(-1.0); + float lat_factH; + float lat_factV; + + lat_factH = sinf(pi/180.0*(*lat)); + lat_factV = cosf(pi/180.0*(*lat)); - *Frhs_u = *Frhs_u + ( corioConstHorz_d*((*vMom)/(*rho)-corioLS_fact_d*(*vBS)/(*rhoBS)) - -corioConstVert_d*((*wMom)/(*rho)-corioLS_fact_d*(*wBS)/(*rhoBS)) ); - *Frhs_v = *Frhs_v - ( corioConstHorz_d*((*uMom)/(*rho)-corioLS_fact_d*(*uBS)/(*rhoBS)) ); - *Frhs_w = *Frhs_w + ( corioConstVert_d*((*uMom)/(*rho)-corioLS_fact_d*(*uBS)/(*rhoBS)) ); + *Frhs_u = *Frhs_u + ( corioConstHorz_d*lat_factH*((*vMom)/(*rho)-corioLS_fact_d*(*vBS)/(*rhoBS)) + -corioConstVert_d*lat_factV*((*wMom)/(*rho)-corioLS_fact_d*(*wBS)/(*rhoBS)) ); + *Frhs_v = *Frhs_v - ( corioConstHorz_d*lat_factH*((*uMom)/(*rho)-corioLS_fact_d*(*uBS)/(*rhoBS)) ); + *Frhs_w = *Frhs_w + ( corioConstVert_d*lat_factV*((*uMom)/(*rho)-corioLS_fact_d*(*uBS)/(*rhoBS)) ); } // end cudaDevice_calcCoriolis() diff --git a/SRC/HYDRO_CORE/CUDA/cuda_coriolisDevice_cu.h b/SRC/HYDRO_CORE/CUDA/cuda_coriolisDevice_cu.h index 96779c7f..a36a0f9a 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_coriolisDevice_cu.h +++ b/SRC/HYDRO_CORE/CUDA/cuda_coriolisDevice_cu.h @@ -43,6 +43,7 @@ extern "C" int cuda_coriolisDeviceCleanup(); */ __device__ void cudaDevice_calcCoriolis(float* Frhs_u, float* Frhs_v, float* Frhs_w, float* rho, float* uMom, float* vMom, float* wMom, - float* rhoBS, float* uBS, float* vBS, float* wBS); + float* rhoBS, float* uBS, float* vBS, float* wBS, + float* lat); #endif // _CORIOLIS_CUDADEV_CU_H diff --git a/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu b/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu index 5c691ad1..3cf070ca 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu +++ b/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -48,6 +49,7 @@ #include "cuda_moistureDevice.cu" #include "cuda_filtersDevice.cu" #include "cuda_cellpertDevice.cu" +#include "cuda_towersDevice.cu" #ifdef URBAN_EXT #include "cuda_urbanDevice.cu" @@ -311,6 +313,7 @@ extern "C" int cuda_hydroCoreDeviceCleanup(){ errorCode = cuda_filtersDeviceCleanup(); } + errorCode = cuda_towersDeviceCleanup(); #ifdef URBAN_EXT /* URBAN */ if (urbanSelector > 0){ @@ -327,6 +330,25 @@ extern "C" int cuda_hydroCoreDeviceCleanup(){ }//end cuda_hydroCoreDeviceCleanup() +/*----->>>>> int cuda_hydroCoreDeviceSecondaryStageSetup(); --------------------------------------------------------- +* Secondary initializations at the device level for BCs and TOWERS submodules +*/ +extern "C" int cuda_hydroCoreDeviceSecondaryStageSetup(float dt, int batchSize){ + int errorCode = CUDA_HYDRO_CORE_SUCCESS; + int BdyUpdateSteps; + + /*Initialize device-level TOWER submodule */ + errorCode = cuda_towersDeviceSetup(batchSize, rank_nTowers, towerInstanceSize, towerSurfInstanceSize); + + /*Compute the number of timesteps between BndyPlane Updates*/ + BdyUpdateSteps = (int) roundf(dtBdyPlaneBCs/dt); + cudaMemcpyToSymbol(BdyUpdateSteps_d, &BdyUpdateSteps, sizeof(int)); + + printf("%d/%d cuda_hydroCoreDeviceSecondaryStageSetup(): BdyUpdateSteps = %d \n",mpi_rank_world,mpi_size_world,BdyUpdateSteps); + fflush(stdout); + return(errorCode); +} + /*----->>>>> extern "C" int cuda_hydroCoreDeviceBuildFrhs(); -------------------------------------------------- * This routine provides the externally callable cuda-kernel call to perform a complete hydroCore build_Frhs */ @@ -508,7 +530,7 @@ extern "C" int cuda_hydroCoreDeviceBuildFrhs(float simTime, int simTime_it, int cudaDevice_hydroCoreComplete<<>>(simTime, simTime_it, dt, timeStage, numRKstages, hydroFlds_d, hydroFldsFrhs_d, hydroFaceVels_d, hydroBaseStateFlds_d, hydroTauFlds_d, sgstkeScalars_d, sgstkeScalarsFrhs_d, moistScalars_d, moistScalarsFrhs_d, moistTauFlds_d, - J13_d, J23_d, J31_d, J32_d, J33_d, invD_Jac_d, zPos_d); + J13_d, J23_d, J31_d, J32_d, J33_d, invD_Jac_d, zPos_d, lat_d); gpuErrchk( cudaGetLastError() ); gpuErrchk( cudaDeviceSynchronize() ); @@ -760,21 +782,37 @@ __global__ void cudaDevice_hydroCoreCommence(int simTime_it, float* hydroFlds_d, fld = &sgstkeScalars_d[fldStride*iFld]; fldBS = &sgstkeScalarsFrhs_d[fldStride*iFld]; // Frhs forcing iwas set to zero, so it can be used here as zero-valued base state if(hydroBCs_d == 1){ //Using LAD BCs - cudaDevice_VerticalAblBCs(iFld, fld, fldBS); - if(rankXid_d == 0){ + cudaDevice_VerticalAblZeroGradBCs(fld); + if (iFld == 0){ // TKE_0 + if(rankXid_d == 0){ + cudaDevice_westBdyBCs(iFld+Nhydro_d, timeWeight, fld, YZBdyPlanes_d, YZBdyPlanesNext_d); + } + if(rankXid_d == numProcsX_d-1){ + cudaDevice_eastBdyBCs(iFld+Nhydro_d, timeWeight, fld, YZBdyPlanes_d, YZBdyPlanesNext_d); + } + if(rankYid_d == 0){ + cudaDevice_southBdyBCs(iFld+Nhydro_d, timeWeight, fld, XZBdyPlanes_d, XZBdyPlanesNext_d); + } + if(rankYid_d == numProcsY_d-1){ + cudaDevice_northBdyBCs(iFld+Nhydro_d, timeWeight, fld, XZBdyPlanes_d, XZBdyPlanesNext_d); + } + cudaDevice_ceilingBdyBCs(iFld+Nhydro_d, timeWeight, fld, XYBdyPlanes_d, XYBdyPlanesNext_d); + }else{ // all other TKE scales + if(rankXid_d == 0){ cudaDevice_lateralTKEBdyBCs(iFld, fld, fldBS, 0); - } - if(rankXid_d == numProcsX_d-1){ + } + if(rankXid_d == numProcsX_d-1){ cudaDevice_lateralTKEBdyBCs(iFld, fld, fldBS, 1); - } - if(rankYid_d == 0){ + } + if(rankYid_d == 0){ cudaDevice_lateralTKEBdyBCs(iFld, fld, fldBS, 2); - } - if(rankYid_d == numProcsY_d-1){ + } + if(rankYid_d == numProcsY_d-1){ cudaDevice_lateralTKEBdyBCs(iFld, fld, fldBS, 3); - } + } + } // end if (iFld == 0) }else if (hydroBCs_d == 2){ - cudaDevice_VerticalAblBCs(iFld, fld, fldBS); // to apply zero-gradient lower boundary BCs + cudaDevice_VerticalAblZeroGradBCs(fld); if(numProcsX_d==1){ cudaDevice_HorizontalPeriodicXdirBCs(iFld, fld); }//periodic and single rank in X-dir --> implies no MPI exchanges made so perform on-device exchange @@ -795,18 +833,18 @@ __global__ void cudaDevice_hydroCoreCommence(int simTime_it, float* hydroFlds_d, if(hydroBCs_d == 1){ //Using LAD BCs cudaDevice_VerticalAblBCs(iFld, fld, fldBS); if(rankXid_d == 0){ - cudaDevice_westBdyBCs(iFld+Nhydro_d, timeWeight, fld, YZBdyPlanes_d, YZBdyPlanesNext_d); + cudaDevice_westBdyBCs(iFld+Nhydro_d+1, timeWeight, fld, YZBdyPlanes_d, YZBdyPlanesNext_d); } if(rankXid_d == numProcsX_d-1){ - cudaDevice_eastBdyBCs(iFld+Nhydro_d, timeWeight, fld, YZBdyPlanes_d, YZBdyPlanesNext_d); + cudaDevice_eastBdyBCs(iFld+Nhydro_d+1, timeWeight, fld, YZBdyPlanes_d, YZBdyPlanesNext_d); } if(rankYid_d == 0){ - cudaDevice_southBdyBCs(iFld+Nhydro_d, timeWeight, fld, XZBdyPlanes_d, XZBdyPlanesNext_d); + cudaDevice_southBdyBCs(iFld+Nhydro_d+1, timeWeight, fld, XZBdyPlanes_d, XZBdyPlanesNext_d); } if(rankYid_d == numProcsY_d-1){ - cudaDevice_northBdyBCs(iFld+Nhydro_d, timeWeight, fld, XZBdyPlanes_d, XZBdyPlanesNext_d); + cudaDevice_northBdyBCs(iFld+Nhydro_d+1, timeWeight, fld, XZBdyPlanes_d, XZBdyPlanesNext_d); } - cudaDevice_ceilingBdyBCs(iFld+Nhydro_d, timeWeight, fld, XYBdyPlanes_d, XYBdyPlanesNext_d); + cudaDevice_ceilingBdyBCs(iFld+Nhydro_d+1, timeWeight, fld, XYBdyPlanes_d, XYBdyPlanesNext_d); }else if (hydroBCs_d == 2){ cudaDevice_VerticalAblZeroGradBCs(fld); // to apply zero-gradient bottom/top BCs if(numProcsX_d==1){ @@ -864,16 +902,18 @@ __global__ void cudaDevice_hydroCoreCommenceRhoInvPresPert(float* hydroFlds_d, f } // end cudaDevice_hydroCoreCommenceRhoInvPresPert() __global__ void cudaDevice_hydroCoreComplete(float simTime, int simTime_it, float dt, int timeStage, int numRKstages, - float* hydroFlds, float* hydroFldsFrhs, - float* hydroFaceVels, float* hydroBaseStateFlds, - float* hydroTauFlds, - float* sgstkeScalars, float* sgstkeScalarsFrhs, - float* moistScalars, float* moistScalarsFrhs, float* moistTauFlds, - float* J13_d, float* J23_d, float* J31_d, float* J32_d, float* J33_d, float* invD_Jac_d, float* zPos_d){ + float* hydroFlds, float* hydroFldsFrhs, + float* hydroFaceVels, float* hydroBaseStateFlds, + float* hydroTauFlds, + float* sgstkeScalars, float* sgstkeScalarsFrhs, + float* moistScalars, float* moistScalarsFrhs, float* moistTauFlds, + float* J13_d, float* J23_d, float* J31_d, float* J32_d, float* J33_d, + float* invD_Jac_d, float* zPos_d, float* lat_d){ - int i,j,k,ijk; + int i,j,k,ijk,ij; int iFld,fldStride; int iStride,jStride,kStride; + int iStride2d,jStride2d; float* rho; float* rho_BS; float* u_cf; @@ -894,6 +934,8 @@ __global__ void cudaDevice_hydroCoreComplete(float simTime, int simTime_it, floa iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); jStride = (Nz_d+2*Nh_d); kStride = 1; + iStride2d = (Ny_d+2*Nh_d); + jStride2d = 1; rho = &hydroFlds[fldStride*RHO_INDX]; rho_BS = &hydroBaseStateFlds[fldStride*RHO_INDX_BS]; u_cf = &hydroFaceVels[fldStride*0]; @@ -902,7 +944,7 @@ __global__ void cudaDevice_hydroCoreComplete(float simTime, int simTime_it, floa if((i >= iMin_d)&&(i < iMax_d) && (j >= jMin_d)&&(j < jMax_d) && - (k >= kMin_d)&&(k < kMax_d) ){ + (k >= kMin_d+3)&&(k < kMax_d) ){ // skipping the first 3 vertical levels for(iFld=0; iFld < Nhydro_d; iFld++){ fld = &hydroFlds[fldStride*iFld]; fldFrhs = &hydroFldsFrhs[fldStride*iFld]; @@ -926,24 +968,6 @@ __global__ void cudaDevice_hydroCoreComplete(float simTime, int simTime_it, floa } else { // defaults to 1st-order upwinding cudaDevice_UpstreamDivAdvFlux(fld, fldFrhs, u_cf, v_cf, w_cf, invD_Jac_d); } - if(iFld==W_INDX){ - if(dampingLayerSelector_d > 0){ // RAYLEIGH DAMPING ON W ******!!!!!!!! - cudaDevice_topRayleighDampingLayerForcing(fld, fldFrhs, - &rho[0], &rho_BS[0], zPos_d); - } //end if dampingLayerSelector > 0 - if(buoyancySelector_d > 0){ // BUOYANCY SOURCE?SINK OF W ******!!!!!!!! - ijk = i*iStride + j*jStride + k*kStride; - if (moistureSelector_d>0){ - if(moistureNvars_d==1){ - cudaDevice_calcBuoyancyMoistNvar1(&fldFrhs[ijk], &rho[ijk], &rho_BS[ijk],&moistScalars[ijk]); - }else if(moistureNvars_d==2){ - cudaDevice_calcBuoyancyMoistNvar2(&fldFrhs[ijk], &rho[ijk], &rho_BS[ijk],&moistScalars[ijk],&moistScalars[fldStride+ijk]); - } - }else{ - cudaDevice_calcBuoyancy(&fldFrhs[ijk], &rho[ijk], &rho_BS[ijk]); - } - } //end if buoyancySelector > 0 - }//end if iFld==W_INDX }//for iFld if ((turbulenceSelector_d>0) && (TKESelector_d>0)){ // : advection of SGSTKE fields for(iFld=0; iFld < TKESelector_d; iFld++){ @@ -1012,9 +1036,96 @@ __global__ void cudaDevice_hydroCoreComplete(float simTime, int simTime_it, floa } } } + }//end if in the range of non-halo cells (skipping the first 3 vertical levels) + // lower the order of adection as the surface is approached + if((i >= iMin_d)&&(i < iMax_d) && + (j >= jMin_d)&&(j < jMax_d) && + (k >= kMin_d)&&(k < kMin_d+3) ){ // only first 3 vertical grid levels + for(iFld=0; iFld < Nhydro_d; iFld++){ + fld = &hydroFlds[fldStride*iFld]; + fldFrhs = &hydroFldsFrhs[fldStride*iFld]; + /* Calculate scalar, cell-valued divergence of the advective flux */ + cudaDevice_HYB34DivAdvFluxX(fld, fldFrhs, u_cf, b_hyb_d, invD_Jac_d); + cudaDevice_HYB34DivAdvFluxY(fld, fldFrhs, v_cf, b_hyb_d, invD_Jac_d); + if (k == kMin_d+2) { // hybrid 3rd-4th order + cudaDevice_HYB34DivAdvFluxZ(fld, fldFrhs, w_cf, b_hyb_d, invD_Jac_d); + } else { // 1st-order upwinding + cudaDevice_UpstreamDivAdvFluxZ(fld, fldFrhs, w_cf, invD_Jac_d); + } + }//for iFld + if ((turbulenceSelector_d>0) && (TKESelector_d>0)){ // : advection of SGSTKE fields + for(iFld=0; iFld < TKESelector_d; iFld++){ + fld = &sgstkeScalars[fldStride*iFld]; + fldFrhs = &sgstkeScalarsFrhs[fldStride*iFld]; + TKEAdvSelector_flag = TKEAdvSelector_d; + TKEAdvSelector_b_hyb_flag = TKEAdvSelector_b_hyb_d; + /* Calculate scalar, cell-valued divergence of the advective flux */ + cudaDevice_HYB34DivAdvFluxX(fld, fldFrhs, u_cf, TKEAdvSelector_b_hyb_flag, invD_Jac_d); + cudaDevice_HYB34DivAdvFluxY(fld, fldFrhs, v_cf, TKEAdvSelector_b_hyb_flag, invD_Jac_d); + if (k == kMin_d+2) { // hybrid 3rd-4th order + cudaDevice_HYB34DivAdvFluxZ(fld, fldFrhs, w_cf, TKEAdvSelector_b_hyb_flag, invD_Jac_d); + } else { // 1st-order upwinding + cudaDevice_UpstreamDivAdvFluxZ(fld, fldFrhs, w_cf, invD_Jac_d); + } + } + } + + if ((moistureSelector_d>0) && (moistureNvars_d>0)){ // : advection of moisture fields + for(iFld=0; iFld < moistureNvars_d; iFld++){ + fld = &moistScalars[fldStride*iFld]; + fldFrhs = &moistScalarsFrhs[fldStride*iFld]; + if (iFld==0){ // water vapor + cudaDevice_HYB34DivAdvFluxX(fld, fldFrhs, u_cf, moistureAdvSelectorQv_b_d, invD_Jac_d); + cudaDevice_HYB34DivAdvFluxY(fld, fldFrhs, v_cf, moistureAdvSelectorQv_b_d, invD_Jac_d); + if (k == kMin_d+2) { // hybrid 3rd-4th order + cudaDevice_HYB34DivAdvFluxZ(fld, fldFrhs, w_cf, moistureAdvSelectorQv_b_d, invD_Jac_d); + } else { // 1st-order upwinding + cudaDevice_UpstreamDivAdvFluxZ(fld, fldFrhs, w_cf, invD_Jac_d); + } + } else { // non-qv moisture species (non-oscillatory schemes) + if (moistureAdvSelectorQi_d == 0) { // 1st-order upstream + cudaDevice_UpstreamDivAdvFlux(fld, fldFrhs, u_cf, v_cf, w_cf, invD_Jac_d); + } else { + cudaDevice_WENO3DivAdvFluxX(fld, fldFrhs, u_cf, invD_Jac_d); + cudaDevice_WENO3DivAdvFluxY(fld, fldFrhs, v_cf, invD_Jac_d); + if (k == kMin_d+2) { // 3rd-order WENO + cudaDevice_WENO3DivAdvFluxZ(fld, fldFrhs, w_cf, invD_Jac_d); + } else { // 1st-order upstream + cudaDevice_UpstreamDivAdvFluxZ(fld, fldFrhs, w_cf, invD_Jac_d); + } + } + } + } + } + }//end if in the range of non-halo cells (only first 3 vertical grid levels) + + if((i >= iMin_d)&&(i < iMax_d) && + (j >= jMin_d)&&(j < jMax_d) && + (k >= kMin_d)&&(k < kMax_d) ){ + // W terms + iFld=W_INDX; + fld = &hydroFlds[fldStride*iFld]; + fldFrhs = &hydroFldsFrhs[fldStride*iFld]; + if(dampingLayerSelector_d > 0){ // RAYLEIGH DAMPING ON W ******!!!!!!!! + cudaDevice_topRayleighDampingLayerForcing(fld, fldFrhs, + &rho[0], &rho_BS[0], zPos_d); + } //end if dampingLayerSelector > 0 + if(buoyancySelector_d > 0){ // BUOYANCY SOURCE?SINK OF W ******!!!!!!!! + ijk = i*iStride + j*jStride + k*kStride; + if (moistureSelector_d>0){ + if(moistureNvars_d==1){ + cudaDevice_calcBuoyancyMoistNvar1(&fldFrhs[ijk], &rho[ijk], &rho_BS[ijk],&moistScalars[ijk]); + }else if(moistureNvars_d==2){ + cudaDevice_calcBuoyancyMoistNvar2(&fldFrhs[ijk], &rho[ijk], &rho_BS[ijk],&moistScalars[ijk],&moistScalars[fldStride+ijk]); + } + }else{ + cudaDevice_calcBuoyancy(&fldFrhs[ijk], &rho[ijk], &rho_BS[ijk]); + } + } //end if buoyancySelector > 0 if(coriolisSelector_d > 0){ ijk = i*iStride + j*jStride + k*kStride; + ij = i*iStride2d + j*jStride2d; cudaDevice_MomentumBS(U_INDX, zPos_d[ijk], &hydroBaseStateFlds[RHO_INDX_BS*fldStride+ijk], &MomBSval[0]); cudaDevice_MomentumBS(V_INDX, zPos_d[ijk], &hydroBaseStateFlds[RHO_INDX_BS*fldStride+ijk], &MomBSval[1]); cudaDevice_MomentumBS(W_INDX, zPos_d[ijk], &hydroBaseStateFlds[RHO_INDX_BS*fldStride+ijk], &MomBSval[2]); @@ -1028,9 +1139,11 @@ __global__ void cudaDevice_hydroCoreComplete(float simTime, int simTime_it, floa &hydroBaseStateFlds[RHO_INDX_BS*fldStride+ijk], &MomBSval[0], &MomBSval[1], - &MomBSval[2]); + &MomBSval[2], + &lat_d[ij]); } //end if coriolisSelector_d > 0 }//end if in the range of non-halo cells + if((turbulenceSelector_d > 0) && ((physics_oneRKonly_d==0) || (timeStage==numRKstages))){ cudaDevice_hydroCoreCalcTurbMixing( &hydroFldsFrhs[fldStride*U_INDX], @@ -1293,7 +1406,7 @@ extern "C" int cuda_hydroCoreInitFieldsDevice(){ * This function handles the synchronization to host of on-device (GPU) fields by executing the appropriate sequence * of cudaMemcpyDeviceiToHost data transfers. */ -extern "C" int cuda_hydroCoreSynchFieldsFromDevice(){ +extern "C" int cuda_hydroCoreSynchFieldsFromDevice(int batchSize){ int errorCode = CUDA_HYDRO_CORE_SUCCESS; int Nelems; int Nelems2d; @@ -1364,6 +1477,12 @@ extern "C" int cuda_hydroCoreSynchFieldsFromDevice(){ } } #endif + /* TOWERS */ + if(rank_nTowers > 0){ + gpuErrchk( cudaMemcpy(towersData, towersData_d, batchSize*rank_nTowers*towerInstanceSize*sizeof(float), cudaMemcpyDeviceToHost) ); + gpuErrchk( cudaMemcpy(towersSurfData, towersSurfData_d, batchSize*rank_nTowers*towerSurfInstanceSize*sizeof(float), cudaMemcpyDeviceToHost) ); + } + gpuErrchk( cudaPeekAtLastError() ); /*Check for errors in the cudaMemCpy calls*/ //#ifdef DEBUG #if 1 @@ -1372,7 +1491,7 @@ extern "C" int cuda_hydroCoreSynchFieldsFromDevice(){ fflush(stdout); MPI_Barrier(MPI_COMM_WORLD); #endif - + return(errorCode); }//end cuda_hydroCoreSynchFieldsFromDevice() diff --git a/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice_cu.h b/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice_cu.h index f6f67c12..df3d9cbe 100644 --- a/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice_cu.h +++ b/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice_cu.h @@ -80,6 +80,9 @@ extern float *hydroRhoInv_d; //storage for 1.0/rho /*---CELL PERTURBATION METHOD*/ #include +/*---TOWERS*/ +#include + #ifdef URBAN_EXT /*URBAN */ #include @@ -118,6 +121,11 @@ extern "C" int cuda_hydroCoreDeviceSetup(); */ extern "C" int cuda_hydroCoreDeviceCleanup(); +/*----->>>>> int cuda_hydroCoreDeviceSecondaryStageSetup(float dt); ----------------------------------------------------------------- +* Secondary initializations at the device level for BCs +*/ +extern "C" int cuda_hydroCoreDeviceSecondaryStageSetup(float dt, int batchSize); + /*----->>>>> extern "C" int cuda_hydroCoreDeviceBuildFrhs(); -------------------------------------------------- * This routine provides the externally callable cuda-kernel call to perform a complete hydroCore build_Frhs */ @@ -162,7 +170,7 @@ __global__ void cudaDevice_hydroCoreComplete(float simTime, int simTime_it, floa float* hydroFaceVels, float* hydroBaseStateFlds, float* hydroTauFlds, float* sgstkeScalars, float* sgstkeScalarsFrhs, float* moistScalars, float* moistScalarsFrhs, float* moistTauFlds, - float* J13_d, float* J23_d, float* J31_d, float* J32_d, float* J33_d, float* invD_Jac_d, float* zPos_d); + float* J13_d, float* J23_d, float* J31_d, float* J32_d, float* J33_d, float* invD_Jac_d, float* zPos_d, float* lat_d); /*----->>>>> __device__ void cudaDevice_SetRhoInv(); -------------------------------------------------- * This is the cuda version of the SetRhoInv routine from the HYDRO_CORE module */ @@ -181,8 +189,8 @@ extern "C" int cuda_hydroCoreInitFieldsDevice(); /*----->>>>> extern "C" int cuda_hydroCoreSynchFieldsFromDevice(); -------------------------------------------------- * This function handles the synchronization to host of on-device (GPU) fields by executing the appropriate sequence -* of cudaMemcpyDeviceiToHost data transfers. +* of cudaMemcpyDeviceiToHost data transfers. This now includes virtual towers. */ -extern "C" int cuda_hydroCoreSynchFieldsFromDevice(); +extern "C" int cuda_hydroCoreSynchFieldsFromDevice(int batchSize); #endif // _HYDRO_CORE_CUDADEV_CU_H diff --git a/SRC/HYDRO_CORE/CUDA/cuda_towersDevice.cu b/SRC/HYDRO_CORE/CUDA/cuda_towersDevice.cu new file mode 100644 index 00000000..0ab84618 --- /dev/null +++ b/SRC/HYDRO_CORE/CUDA/cuda_towersDevice.cu @@ -0,0 +1,172 @@ +/* FastEddy®: SRC/HYDRO_CORE/CUDA/cuda_towersDevice.cu +* ©2016 University Corporation for Atmospheric Research +* +* This file is licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + + +/*---TOWERS*/ +__constant__ int rank_nTowers_d; /*Number of towers within an mpi_rank_world subdomain*/ +__constant__ int towerInstanceSize_d; /*size of a timestep instance of all tower data*/ +__constant__ int towerSurfInstanceSize_d; /*size of a timestep instance of all tower surface data*/ +float* towersData_d; /* Data Structure to store virtual tower data on device*/ +float* towersSurfData_d; /* Data Structure to store virtual tower surface data on device*/ +int* tower_iInds_d; /*rank-centric i-indices */ +int* tower_jInds_d; /*rank-centric j-indices */ + +/*#################------------ TOWERS submodule function definitions ------------------#############*/ +/*----->>>>> int cuda_towersDeviceSetup(); --------------------------------------------------------- +* Used to cudaMemcpy parameters and allocate arrays for the TOWERS_CUDA submodule. +*/ +extern "C" int cuda_towersDeviceSetup(int NtBatch, int rank_nTowers, int towerInstanceSize, int towerSurfInstanceSize){ + int errorCode = CUDA_TOWERS_SUCCESS; + size_t NelemsTowers; + size_t NelemsTowersSurf; + + cudaMemcpyToSymbol(rank_nTowers_d, &rank_nTowers, sizeof(int)); + cudaMemcpyToSymbol(towerInstanceSize_d, &towerInstanceSize, sizeof(int)); + cudaMemcpyToSymbol(towerSurfInstanceSize_d, &towerSurfInstanceSize, sizeof(int)); + + if(rank_nTowers > 0){ + NelemsTowers = (size_t)(NtBatch*rank_nTowers*towerInstanceSize); + fecuda_DeviceMalloc(NelemsTowers, &towersData_d); + NelemsTowersSurf = (size_t)(NtBatch*rank_nTowers*towerSurfInstanceSize); + fecuda_DeviceMalloc(NelemsTowersSurf, &towersSurfData_d); + fecuda_DeviceMallocInt(rank_nTowers, &tower_iInds_d); + cudaMemcpy(tower_iInds_d, tower_iInds, rank_nTowers*sizeof(int), cudaMemcpyHostToDevice); + fecuda_DeviceMallocInt(rank_nTowers, &tower_jInds_d); + cudaMemcpy(tower_jInds_d, tower_jInds, rank_nTowers*sizeof(int), cudaMemcpyHostToDevice); + } + printf("cuda_towersDeviceSetup: NelemsTowers = %d, NelemsTowersSurf = %d\n",NelemsTowers,NelemsTowersSurf); + fflush(stdout); + return(errorCode); +} //end cuda_towersDeviceSetup() + +/*----->>>>> extern "C" int cuda_towersDeviceCleanup(); ----------------------------------------------------------- +Used to free all malloced memory by the TOWERS submodule. +*/ + +extern "C" int cuda_towersDeviceCleanup(){ + int errorCode = CUDA_TOWERS_SUCCESS; + + /* Free any TOWERS submodule arrays */ + if(rank_nTowers > 0){ + cudaFree(towersData_d); + cudaFree(towersSurfData_d); + cudaFree(tower_iInds_d); + cudaFree(tower_jInds_d); + + } + return(errorCode); + +}//end cuda_towersDeviceCleanup() + +/*----->>>>> __global__ void cudaDevice_towerAppendBuffers(); ------------------------------------------ +* This is the gloabl-entry kernel routine for updating tower device-sided buffers +*/ +__global__ void cudaDevice_towerAppendBuffers(int itBatch, int batchSize, + float *towersData_d, float *towersSurfData_d, int *tower_iInds_d, int *tower_jInds_d, + int Nhydro, float *hydroFlds_d, + int Nsgstke, float *sgstkeScalars_d, + int Nmoist, float *moistScalars_d, + int NauxSc, float *hydroAuxScalars_d, + int Ntaus, float *hydroTauFlds_d, + int NtausMoist, float *moistTauFlds_d, + float *z0m_d, float *z0t_d, float *tskin_d, float *qskin_d, + float *fricVel_d, float *invOblen_d, float *htFlux_d, float *qFlux_d){ + + int towerCount; + int i,j,k; + int ijk; + int ij; + int iStride,jStride,kStride; + int iFld,fldStride; + int towerBaseAddress; + int towerSurfBaseAddress; + int towerFld_size; + int towerFld_cnt; + int towIndx; + /*Establish necessary indices for spatial locality*/ + i = (blockIdx.x)*blockDim.x + threadIdx.x; + j = (blockIdx.y)*blockDim.y + threadIdx.y; + k = (blockIdx.z)*blockDim.z + threadIdx.z; + iStride = (Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + jStride = (Nz_d+2*Nh_d); + kStride = 1; + fldStride = (Nx_d+2*Nh_d)*(Ny_d+2*Nh_d)*(Nz_d+2*Nh_d); + towerFld_size = Nz_d; + + for(towerCount = 0; towerCount < rank_nTowers_d; towerCount++){ + if((i == tower_iInds_d[towerCount])&& + (j == tower_jInds_d[towerCount]) && + (k >= kMin_d)&&(k < kMax_d) ){ + + towerBaseAddress = towerCount*(batchSize*towerInstanceSize_d); + towerFld_cnt = 0; + ijk = i*iStride + j*jStride + k*kStride; + ij = i*(Ny_d+2*Nh_d) + j; + for(iFld=0; iFld < Nhydro; iFld++){ + towIndx = towerBaseAddress + itBatch*towerInstanceSize_d + towerFld_cnt*towerFld_size + k-Nh_d; + towersData_d[towIndx] = hydroFlds_d[iFld*fldStride+ijk]; + towerFld_cnt += 1; + } + for(iFld=0; iFld < Nsgstke; iFld++){ + towIndx = towerBaseAddress + itBatch*towerInstanceSize_d + towerFld_cnt*towerFld_size + k-Nh_d; + towersData_d[towIndx] = sgstkeScalars_d[iFld*fldStride+ijk]; + towerFld_cnt += 1; + } + for(iFld=0; iFld < Nmoist; iFld++){ + towIndx = towerBaseAddress + itBatch*towerInstanceSize_d + towerFld_cnt*towerFld_size + k-Nh_d; + towersData_d[towIndx] = moistScalars_d[iFld*fldStride+ijk]; + towerFld_cnt += 1; + } + for(iFld=0; iFld < NauxSc; iFld++){ + towIndx = towerBaseAddress + itBatch*towerInstanceSize_d + towerFld_cnt*towerFld_size + k-Nh_d; + towersData_d[towIndx] = hydroAuxScalars_d[iFld*fldStride+ijk]; + towerFld_cnt += 1; + } + for(iFld=0; iFld < Ntaus; iFld++){ + towIndx = towerBaseAddress + itBatch*towerInstanceSize_d + towerFld_cnt*towerFld_size + k-Nh_d; + towersData_d[towIndx] = hydroTauFlds_d[iFld*fldStride+ijk]; + towerFld_cnt += 1; + } + for(iFld=0; iFld < NtausMoist; iFld++){ + towIndx = towerBaseAddress + itBatch*towerInstanceSize_d + towerFld_cnt*towerFld_size + k-Nh_d; + towersData_d[towIndx] = moistTauFlds_d[iFld*fldStride+ijk]; + towerFld_cnt += 1; + } + if(k == kMin_d){ + towerSurfBaseAddress = towerCount*(batchSize*towerSurfInstanceSize_d); + towIndx = towerSurfBaseAddress + itBatch*towerSurfInstanceSize_d; + towersSurfData_d[towIndx] = z0m_d[ij]; + towIndx += 1; //only a single surface value so increment by 1 + towersSurfData_d[towIndx] = z0t_d[ij]; + towIndx += 1; //only a single surface value so increment by 1 + towersSurfData_d[towIndx] = tskin_d[ij]; + towIndx += 1; //only a single surface value so increment by 1 + towersSurfData_d[towIndx] = fricVel_d[ij]; + towIndx += 1; //only a single surface value so increment by 1 + towersSurfData_d[towIndx] = invOblen_d[ij]; + towIndx += 1; //only a single surface value so increment by 1 + towersSurfData_d[towIndx] = htFlux_d[ij]; + towIndx += 1; //only a single surface value so increment by 1 + if(Nmoist > 0){ + towersSurfData_d[towIndx] = qskin_d[ij]; + towIndx += 1; //only a single surface value so increment by 1 + towersSurfData_d[towIndx] = qFlux_d[ij]; + towIndx += 1; //only a single surface value so increment by 1 + }//end if Nmoist > 0 + }//end if k == kMin_d + }//end if i= tower_iInds_d[towerCount] && ... + }//end for towerCount +}//end cudaDevice_towerAppendBuffers diff --git a/SRC/HYDRO_CORE/CUDA/cuda_towersDevice_cu.h b/SRC/HYDRO_CORE/CUDA/cuda_towersDevice_cu.h new file mode 100644 index 00000000..cb1d96bb --- /dev/null +++ b/SRC/HYDRO_CORE/CUDA/cuda_towersDevice_cu.h @@ -0,0 +1,59 @@ +/* FastEddy®: SRC/HYDRO_CORE/CUDA/cuda_towersDevice_cu.h +* ©2016 University Corporation for Atmospheric Research +* +* This file is licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ +#ifndef _TOWERS_CUDADEV_CU_H +#define _TOWERS_CUDADEV_CU_H + +/*coriolis return codes */ +#define CUDA_TOWERS_SUCCESS 0 + +/*##############------------------- TOWERS submodule variable declarations ---------------------#################*/ +/*---TOWERS*/ +extern __constant__ int rank_nTowers_d; /*Number of towers within an mpi_rank_world subdomain*/ +extern __constant__ int towerInstanceSize_d; /*size of a timestep instance of all tower data*/ +extern __constant__ int towerSurfInstanceSize_d; /*size of a timestep instance of all tower surface data*/ +extern float* towersData_d; /* Data Structure to store virtual tower data on device*/ +extern float* towersSurfData_d; /* Data Structure to store virtual tower surface ata on device*/ +extern int* tower_iInds_d; /*rank-centric i-indices */ +extern int* tower_jInds_d; /*rank-centric j-indices */ + +/*##############-------------- TOWERS_CUDADEV submodule function declarations ------------------############*/ + +/*----->>>>> int cuda_towersDeviceSetup(); --------------------------------------------------------- +* Used to cudaMemcpy parameters and allocate arrays for the TOWERS_CUDA submodule. +*/ +extern "C" int cuda_towersDeviceSetup(int NtBatch, int rank_nTowers, int towerInstanceSize, int towerSurfInstanceSize); + +/*----->>>>> extern "C" int cuda_towersDeviceCleanup(); ----------------------------------------------------------- +* Used to free all malloced memory by the TOWERS submodule. +*/ +extern "C" int cuda_towersDeviceCleanup(); + +/*----->>>>> __global__ void cudaDevice_towerAppendBuffers(); ------------------------------------------ +* This is the gloabl-entry kernel routine for updating tower device-sided buffers +*/ +//__global__ void cudaDevice_towerAppendBuffers(int itBatch, int NtBatch, +__global__ void cudaDevice_towerAppendBuffers(int itBatch, int batchSize, + float *towersData_d, float *towersSurfData_d, int *tower_iInds_d, int *tower_jInds_d, + int Nhydro, float *hydroFlds_d, + int Nsgstke, float *sgstkeScalars_d, + int Nmoist, float *moistScalars_d, + int NauxSc, float *hydroAuxScalars_d, + int Ntaus, float *hydroTauFlds_d, + int NtausMoist, float *moistTauFlds_d, + float *z0m_d, float *z0t_d, float *tskin_d, float *qskin_d, + float *fricVel_d, float *invOblen_d, float *htFlux_d, float *qFlux_d); + +#endif // _TOWERS_CUDADEV_CU_H diff --git a/SRC/HYDRO_CORE/hydro_core.c b/SRC/HYDRO_CORE/hydro_core.c index 33084d36..f398577b 100644 --- a/SRC/HYDRO_CORE/hydro_core.c +++ b/SRC/HYDRO_CORE/hydro_core.c @@ -283,6 +283,20 @@ float thetaAmplitude; /* Initial theta perturbation (maximum amplitude in K)*/ int physics_oneRKonly; /* selector to apply physics RHS forcing only at the latest RK stage */ +/*---VIRTUAL TOWERS*/ +int *towerIDs; +int *tower_iInds; +int *tower_jInds; +int rank_nTowers; +float *tower_xOffsets; +float *tower_yOffsets; +double *tower_LonOffsets; +double *tower_LatOffsets; +int towerInstanceSize; +int towerSurfInstanceSize; +float *towersData; +float *towersSurfData; + /*###################------------------- HYDRO_CORE module function definitions ---------------------#################*/ /*----->>>>> int hydro_coreGetParams(); ---------------------------------------------------------------------- @@ -325,6 +339,7 @@ int hydro_coreGetParams(){ errorCode = queryIntegerParameter("TKEAdvSelector", &TKEAdvSelector, 0, 6, PARAM_OPTIONAL); TKEAdvSelector_b_hyb = 0.0; //Default to 0.0 errorCode = queryFloatParameter("TKEAdvSelector_b_hyb", &TKEAdvSelector_b_hyb, 0.0, 1.0, PARAM_OPTIONAL); + if (turbulenceSelector == 1){ if (TKESelector == 0){ c_s = 0.18; //Default to 0.18 @@ -430,14 +445,13 @@ int hydro_coreGetParams(){ // cellpertSelector = 0; // Default to off errorCode = queryIntegerParameter("cellpertSelector", &cellpertSelector, 0, 1, PARAM_OPTIONAL); - cellpert_nts = 500; // Default to 500 time steps - errorCode = queryIntegerParameter("cellpert_nts", &cellpert_nts, 0, 1e+6, PARAM_OPTIONAL); if (cellpertSelector > 0){ - errorCode = queryIntegerParameter("cellpertSelector", &cellpertSelector, 0, 1, PARAM_OPTIONAL); cellpert_sw2b = 0; // Default to 0 errorCode = queryIntegerParameter("cellpert_sw2b", &cellpert_sw2b, 0, 3, PARAM_OPTIONAL); cellpert_amp = 0.5; // Default to 0.5 K errorCode = queryFloatParameter("cellpert_amp", &cellpert_amp, 0.0, 20.0, PARAM_OPTIONAL); + cellpert_nts = 500; // Default to 500 time steps + errorCode = queryIntegerParameter("cellpert_nts", &cellpert_nts, 0, 1e+6, PARAM_OPTIONAL); cellpert_gppc = 8; // Default to 8 grid points per cell errorCode = queryIntegerParameter("cellpert_gppc", &cellpert_gppc, 0, 50, PARAM_OPTIONAL); cellpert_ndbc = 3; // Default to 3 cells @@ -446,15 +460,15 @@ int hydro_coreGetParams(){ errorCode = queryIntegerParameter("cellpert_kbottom", &cellpert_kbottom, 1, 10, PARAM_OPTIONAL); cellpert_ktop = 20; // Default to 20th grid point above surface errorCode = queryIntegerParameter("cellpert_ktop", &cellpert_ktop, 0, 200, PARAM_OPTIONAL); - if (cellpert_ktop > Nz){ - cellpert_ktop = Nz-10; - } cellpert_tvcp = 0; // Default to off errorCode = queryIntegerParameter("cellpert_tvcp", &cellpert_tvcp, 0, 1, PARAM_OPTIONAL); cellpert_eckert = 0.2; // Default to Ec = 0.2 errorCode = queryFloatParameter("cellpert_eckert", &cellpert_eckert, 0.0, 10.0, PARAM_OPTIONAL); cellpert_tsfact = 1.0; // Default to cellpert_tsfact = 1.0 errorCode = queryFloatParameter("cellpert_tsfact", &cellpert_tsfact, 0.0, 10.0, PARAM_OPTIONAL); + if (cellpert_ktop > Nz){ + cellpert_ktop = Nz; + } } // lsfSelector = 0; // Default to off @@ -527,7 +541,9 @@ int hydro_coreGetParams(){ errorCode = queryFloatParameter("moistureCondTscale", &moistureCondTscale, 1e-4, 1000.0, PARAM_MANDATORY); errorCode = queryIntegerParameter("moistureCondBasePres", &moistureCondBasePres, 0, 1, PARAM_MANDATORY); errorCode = queryFloatParameter("moistureMPcallTscale", &moistureMPcallTscale, 1e-4, 1000.0, PARAM_MANDATORY); - errorCode = queryFloatParameter("surflayer_wq", &surflayer_wq, -5e+0, 5e+0, PARAM_MANDATORY); + if (surflayerSelector == 1){ + errorCode = queryFloatParameter("surflayer_wq", &surflayer_wq, -5e+0, 5e+0, PARAM_MANDATORY); + } if (surflayerSelector == 2){ errorCode = queryFloatParameter("surflayer_qr", &surflayer_qr, -1e+1, 1e+1, PARAM_MANDATORY); errorCode = queryIntegerParameter("surflayer_qskin_input", &surflayer_qskin_input, 0, 1, PARAM_OPTIONAL); @@ -615,7 +631,7 @@ int hydro_coreGetParams(){ } //endif srcAuxScFile == NULL... }// endif NhydroAuxScalars > 0 stabilityScheme = 2; //Default to 2 - errorCode = queryIntegerParameter("stabilityScheme", &stabilityScheme, 2, 2, PARAM_MANDATORY); + errorCode = queryIntegerParameter("stabilityScheme", &stabilityScheme, 1, 2, PARAM_MANDATORY); temp_grnd = 300.0; //Default to 300.0-(Kelvin) = 80.33-(Fahrenheit) = 26.85-(Celsius) errorCode = queryFloatParameter("temp_grnd", &temp_grnd, FLT_MIN, FLT_MAX, PARAM_MANDATORY); pres_grnd = 1.0e5; //Default to refPressure 100,000-(pascals) = 1000-(millibars) @@ -678,7 +694,6 @@ int hydro_coreInit(){ char moistName[MAX_HC_FLDNAME_LENGTH]; char moistName_base[MAX_HC_FLDNAME_LENGTH]; char moistName_tmp[MAX_HC_FLDNAME_LENGTH]; - float pi; int fldStride; float z1oz0,z1,z1ozt0; int strLength; @@ -840,9 +855,9 @@ int hydro_coreInit(){ printParameter("stabilityScheme", "Scheme used to set hydrostatic, stability-dependent Base-State EOS fields"); printParameter("temp_grnd", "Air Temperature (K) at the ground used to set hydrostatic Base-State EOS fields"); printParameter("pres_grnd", "Pressure (Pa) at the ground used to set hydrostatic Base-State EOS fields"); - printParameter("zStableBottom", "Height (m) of the first stable upper-layer when stabilityScheme = 1 or 2"); + printParameter("zStableBottom", "Height (m) of the first stable upper-layer when stabilityScheme = 2"); printParameter("stableGradient", - "Vertical gradient (K/m) of the first stable upper-layer when stabilityScheme = 1 or 2"); + "Vertical gradient (K/m) of the first stable upper-layer when stabilityScheme = 2"); printParameter("zStableBottom2", "Height (m) of the second stable upper-layer when stabilityScheme = 2"); printParameter("zStableBottom3", "Height (m) of the third stable upper-layer when stabilityScheme = 2"); printParameter("stableGradient2", @@ -1343,7 +1358,18 @@ int hydro_coreInit(){ errorCode=GADInit(); #endif MPI_Barrier(MPI_COMM_WORLD); - + + /*Initialize lat & lon arrays if no initial condition file was provided (cold-start) */ + if(inFile == NULL){ + for(i=iMin-Nh; i < iMax+Nh; i++){ + for(j=jMin-Nh; j < jMax+Nh; j++){ + ij = i*(Nyp+2*Nh)+j; + lat[ij] = coriolisLatitude; + lon[ij] = 0.0; // longitude is zero in idealized fresh start runs + } + } + } + /* Provide intial approximation for the momentum and heat exchange coefficient at all surface locations*/ k = kMin; for(i=iMin-Nh; i < iMax+Nh; i++){ @@ -1521,12 +1547,12 @@ int hydro_coreInit(){ fflush(stdout); } if( moistureSelector > 0){ - nBndyVars = Nhydro+moistureNvars; + nBndyVars = Nhydro+1+moistureNvars; // +1 is for TKE_0 nSurfBndyVars = 2; //Only allows tskin and qskin }else{ - nBndyVars = Nhydro; + nBndyVars = Nhydro+1; // +1 is for TKE_0 nSurfBndyVars = 1; //Only allows tskin - } //end if moisture is on else not //NOTE: Doesn't handle any AuxScalars or TKE-related Prog. variables. + } //end if moisture is on else not //NOTE: Doesn't handle any AuxScalars Prog. variables. XZBdyPlanesGlobal = (float *) malloc( 2*(nBndyVars)*Nx*Nz*sizeof(float) ); YZBdyPlanesGlobal = (float *) malloc( 2*(nBndyVars)*Ny*Nz*sizeof(float) ); XYBdyPlanesGlobal = (float *) malloc( 2*(nBndyVars)*Nx*Ny*sizeof(float) ); @@ -1558,11 +1584,10 @@ int hydro_coreInit(){ Rv_Rg = R_vapor/R_gas; /* Ratio R_vapor/R_gas*/ /* Coriolis-term constants */ - pi = acos(-1); if(coriolisSelector > 0){ - corioConstHorz = 1.45842e-4*sin(pi/180.0*coriolisLatitude); //1.45842e-4 = 2*Earth-Omega + corioConstHorz = 1.45842e-4; //1.45842e-4 = 2*Earth-Omega if(coriolisSelector > 1){ - corioConstVert = 1.45842e-4*cos(pi/180.0*coriolisLatitude); + corioConstVert = 1.45842e-4; }else{ corioConstVert = 0.0; } //end if vert @@ -1736,7 +1761,7 @@ int hydro_coreSetBaseState(){ rhoBase = &hydroBaseStateFlds[RHO_INDX_BS*fldStride]; thetaBase = &hydroBaseStateFlds[THETA_INDX_BS*fldStride]; /* ----Based on stabilityScheme setup Base-State rho,theta, and pressure profiles */ - if(stabilityScheme == 0){ /* None, constant density, theta (potential temperature), and pressure fields */ + if(stabilityScheme == 1){ /* None, constant density, theta (potential temperature), and pressure fields -> laboratory scale simulations */ for(i=iMin-Nh; i < iMax+Nh; i++){ // Cover the halos in X for(j=jMin-Nh; j < jMax+Nh; j++){ // Cover the halos in Y for(k=kMin-Nh; k < kMax+Nh; k++){ // Cover the halos in Z @@ -1747,38 +1772,6 @@ int hydro_coreSetBaseState(){ } //end for(k... } // end for(j... } // end for(i... - printf("stabilityScheme == 0: Base State setup complete.\n"); - }else if(stabilityScheme == 1){ /* stable linear potential temperature profile above some height zStableBottom, - neutral below zStableBottom*/ - for(i=iMin-Nh; i < iMax+Nh; i++){ // Cover the halos in X - for(j=jMin-Nh; j < jMax+Nh; j++){ // Cover the halos in Y - for(k=kMin-Nh; k < kMax+Nh; k++){ // Cover the halos in Z - ijk = i*(Nyp+2*Nh)*(Nzp+2*Nh)+j*(Nzp+2*Nh)+k; - if(zPos[ijk] <= zStableBottom){ //This point is within the neutral lower-layer - thetaBase[ijk] = theta_grnd; - hydroBaseStatePres[ijk] = refPressure*pow( (-accel_g/cp_gas)*( zPos[ijk]/theta_grnd ) - +pow(pres_grnd/refPressure,R_cp) //base of the first pow (...) - ,cp_R); //exponent of the first pow(...) - }else{ //This point is within the stable upper-layer - //Set theta - thetaBase[ijk] = theta_grnd + stableGradient*(zPos[ijk]-zStableBottom); - //set base state pressure - hydroBaseStatePres[ijk] = refPressure*pow( (-accel_g/cp_gas)*( zStableBottom/theta_grnd - +(1.0/stableGradient)*log(1.0+stableGradient*(zPos[ijk]-zStableBottom)/theta_grnd)) - +pow(pres_grnd/refPressure,R_cp) //base of the first pow (...) - ,cp_R); //exponent of the first pow(...) - } //end zPos[ijk >= zStableBottom - //back out base state air temperature - BS_Temp = thetaBase[ijk]*pow( hydroBaseStatePres[ijk]/refPressure,R_cp); - //back out base state density - rhoBase[ijk] = hydroBaseStatePres[ijk]/(BS_Temp*R_gas); - //Given this density set the flux form of the potential temperature prognostic field (rho*theta) - thetaBase[ijk] = thetaBase[ijk]*rhoBase[ijk]; - //Finally recast the base state pressure in a "discretisation-consistent" manner - hydroBaseStatePres[ijk] = pow(thetaBase[ijk]*constant_1, cp_cv); //This minimizes round off under the pressure formulation in calcPerturbationPRessure() - } //end for(k... - } // end for(j... - } // end for(i... printf("stabilityScheme == 1: Base State setup complete.\n"); }else if(stabilityScheme == 2){ for(i=iMin-Nh; i < iMax+Nh; i++){ // Cover the halos in X @@ -1830,25 +1823,6 @@ int hydro_coreSetBaseState(){ } // end for(j... } // end for(i... printf("stabilityScheme == 2: Base State setup complete.\n"); - }else if(stabilityScheme == 3){ - printf("stabilityScheme == 3: ERROR: No scheme implemented for stabilityScheme == 3, use instead 1, 2, or 4!! \n"); - }else if(stabilityScheme == 4){ /*Experimental setup for constant rho and constant theta profiles. - Use only for total domain vertical extent < 10m. */ - rho_grnd = 1.1; - for(i=iMin-Nh; i < iMax+Nh; i++){ // Cover the halos in X - for(j=jMin-Nh; j < jMax+Nh; j++){ // Cover the halos in Y - for(k=kMin-Nh; k < kMax+Nh; k++){ // Cover the halos in Z - ijk = i*(Nyp+2*Nh)*(Nzp+2*Nh)+j*(Nzp+2*Nh)+k; - if(zPos[ijk] <= 0.5){ - rhoBase[ijk] = rho_grnd; - }else{ - rhoBase[ijk] = rho_grnd; - } - thetaBase[ijk] = rho_grnd*theta_grnd; - hydroBaseStatePres[ijk] = pow(thetaBase[ijk]*constant_1,cp_cv); - } //end for(k... - } // end for(j... - } // end for(i... } //end if-else... stabilityScheme... if(inFile == NULL){ @@ -1910,34 +1884,168 @@ int hydro_coreSetBaseState(){ } // end for(j... } // end for(i... } //endif thetaPerturbationSwitch==1 - - }else{ //Initial conditions were provided... - if(stabilityScheme==3){ - for(iFld=0; iFld < 2; iFld++){ - switch (iFld){ - case 0: - fldBase = &hydroFlds[RHO_INDX*fldStride]; - fldBaseBS = &hydroBaseStateFlds[RHO_INDX_BS*fldStride]; - break; - case 1: - fldBase = &hydroFlds[THETA_INDX*fldStride]; - fldBaseBS = &hydroBaseStateFlds[THETA_INDX_BS*fldStride]; - break; - } - for(i=iMin-Nh; i < iMax+Nh; i++){ // Cover the halos in X - for(j=jMin-Nh; j < jMax+Nh; j++){ // Cover the halos in Y - for(k=kMin-Nh; k < kMax+Nh; k++){ // Cover the halos in Z - ijk = i*(Nyp+2*Nh)*(Nzp+2*Nh)+j*(Nzp+2*Nh)+k; - fldBaseBS[ijk] = fldBase[ijk]; - } //end for(k... - } // end for(j... - } // end for(i... - }//end if-else iFld==0 - }//end if stabilityScheme==3 }//If no initial conditions were specified return(errorCode); }// end coreSetBaseState +/*----->>>>> int hydro_coreAllocateTowersDataStructure(); --------------------------------------------------- +* Utility to allocate virtual tower data structures on appropriate ranks +*/ +int hydro_coreAllocateTowersDataStructure(int nProfs, ioProfiles_t towProfs, int NtBatch){ + int errorCode = HYDRO_CORE_SUCCESS; + int itower; + int nElems; + int nSurfElems; + int towerCount; + int i,j,k,ij,ijk; + int iStride,jStride,kStride,fldStride; + int towerBaseAddress; + int towerSurfBaseAddress; + int towerFld_size; + int towerFld_cnt; + int towIndx; + int iFld; + + rank_nTowers = 0; + towerInstanceSize = Nz*(registered3dVars-4); // r3dV-4 since no x,y,zPos, or pressure + towerSurfInstanceSize = (registered2dVars-(3+surflayer_offshore)); // r2dV-(3+surflayer_offshore) since no (topoPos, lat, lon + sea_mask) + //Count the number of towers in a given mpi_rank's subdomain + for(itower = 0; itower < nProfs; itower++){ + if(towProfs.mpi_ranks[itower]==mpi_rank_world){ + rank_nTowers = rank_nTowers + 1; + } + } + if(rank_nTowers > 0){ + //Allocate and set the per-rank towerIDs + towerIDs = (int *) malloc(rank_nTowers*sizeof(int)); + towerCount=0; + for(itower = 0; itower < nProfs; itower++){ + if(towProfs.mpi_ranks[itower]==mpi_rank_world){ + towerIDs[towerCount]=towProfs.profIDs[itower]; + towerCount=towerCount+1; + } + } + + //Calculate the number of float data elements + nElems = rank_nTowers*NtBatch*towerInstanceSize; + nSurfElems = rank_nTowers*NtBatch*towerSurfInstanceSize; + //Allocate the tower data structure + towersData = (float *) malloc(nElems*sizeof(float)); + towersSurfData = (float *) malloc(nSurfElems*sizeof(float)); + printf("%d/%d: NtBatch = %d, rank_nTowers = %d, towerInstanceSize = %d, nElems = %d, towerSurfInstanceSize = %d, nSurfElems = %d\n", + mpi_rank_world,mpi_size_world,NtBatch,rank_nTowers,towerInstanceSize,nElems,towerSurfInstanceSize,nSurfElems); + + //Now identify the mpi_rank-specific i,j indices for each tower in the mpi_rank's subdomain + tower_iInds = (int *) malloc(rank_nTowers*sizeof(int)); + tower_jInds = (int *) malloc(rank_nTowers*sizeof(int)); + if(towerProfiles.coordType == 0){ + tower_LonOffsets = (double *) malloc(rank_nTowers*sizeof(double)); + tower_LatOffsets = (double *) malloc(rank_nTowers*sizeof(double)); + }else{ + tower_xOffsets = (float *) malloc(rank_nTowers*sizeof(float)); + tower_yOffsets = (float *) malloc(rank_nTowers*sizeof(float)); + } + for(towerCount = 0; towerCount < rank_nTowers; towerCount++){ + //Call an index finding function from the grid module. + if(towerProfiles.coordType == 0){ + errorCode = gridGetIJindsFromLatLonPosition(towerProfiles.coordsLon[towerIDs[towerCount]], + towerProfiles.coordsLat[towerIDs[towerCount]], + &tower_iInds[towerCount], &tower_jInds[towerCount]); + errorCode = gridGetLatLonOffsetsFromCellIndices(towerProfiles.coordsLon[towerIDs[towerCount]], + towerProfiles.coordsLat[towerIDs[towerCount]], + tower_iInds[towerCount], tower_jInds[towerCount], + &tower_LonOffsets[towerCount], &tower_LatOffsets[towerCount]); + printf("%d/%d: towerCount = %d, towerID = %d, (Lat,Lon) = (%f,%f), (iInd,jInd) = (%d,%d), (LatOffset,LonOffset) = (%f,%f))\n", + mpi_rank_world,mpi_size_world,towerCount,towerIDs[towerCount], + towerProfiles.coordsLat[towerIDs[towerCount]],towerProfiles.coordsLon[towerIDs[towerCount]], + tower_iInds[towerCount],tower_jInds[towerCount],tower_LatOffsets[towerCount], tower_LonOffsets[towerCount]); + }else{ + errorCode = gridGetIJindsFromXYPosition(towerProfiles.coordsWE[towerIDs[towerCount]], + towerProfiles.coordsSN[towerIDs[towerCount]], + &tower_iInds[towerCount], &tower_jInds[towerCount]); + errorCode = gridGetXYOffsetsFromCellIndices(towerProfiles.coordsWE[towerIDs[towerCount]], + towerProfiles.coordsSN[towerIDs[towerCount]], + tower_iInds[towerCount], tower_jInds[towerCount], + &tower_xOffsets[towerCount], &tower_yOffsets[towerCount]); + printf("%d/%d: towerCount = %d, towerID = %d, (x,y) = (%f,%f), (iInd,jInd) = (%d,%d), (xOffset,yOffset) = (%f,%f))\n", + mpi_rank_world,mpi_size_world,towerCount,towerIDs[towerCount], + towerProfiles.coordsWE[towerIDs[towerCount]],towerProfiles.coordsSN[towerIDs[towerCount]], + tower_iInds[towerCount],tower_jInds[towerCount],tower_xOffsets[towerCount], tower_yOffsets[towerCount]); + }//end if coordType == 0, else + fflush(stdout); + } + + //Initialize the towerData and towerSurfData values + iStride = (Nyp+2*Nh)*(Nzp+2*Nh); + jStride = (Nzp+2*Nh); + kStride = 1; + + fldStride = (Nxp+2*Nh)*(Nyp+2*Nh)*(Nzp+2*Nh); + towerFld_size = Nzp; + for(towerCount = 0; towerCount < rank_nTowers; towerCount++){ + towerBaseAddress = towerCount*(NtBatch*towerInstanceSize); + + i = tower_iInds[towerCount]; + j = tower_jInds[towerCount]; + + for(k=kMin; k < kMax; k++){ + towerFld_cnt = 0; + ijk = i*iStride + j*jStride + k*kStride; + for(iFld=0; iFld < Nhydro; iFld++){ + towIndx = towerBaseAddress + towerFld_cnt*towerFld_size + k-Nh; + towersData[towIndx] = hydroFlds[iFld*fldStride+ijk]; + towerFld_cnt += 1; + }//end for iFld + for(iFld=0; iFld < TKESelector*turbulenceSelector; iFld++){ + towIndx = towerBaseAddress + towerFld_cnt*towerFld_size + k-Nh; + towersData[towIndx] = sgstkeScalars[iFld*fldStride+ijk]; + towerFld_cnt += 1; + } + for(iFld=0; iFld < moistureNvars*moistureSelector; iFld++){ + towIndx = towerBaseAddress + towerFld_cnt*towerFld_size + k-Nh; + towersData[towIndx] = moistScalars[iFld*fldStride+ijk]; + towerFld_cnt += 1; + } + for(iFld=0; iFld < NhydroAuxScalars; iFld++){ + towIndx = towerBaseAddress + towerFld_cnt*towerFld_size + k-Nh; + towersData[towIndx] = hydroAuxScalars[iFld*fldStride+ijk]; + towerFld_cnt += 1; + } + for(iFld=0; iFld < hydroSubGridWrite*9; iFld++){ //There are 6 Tau^i-j and 3 tau^Theta-j + towIndx = towerBaseAddress + towerFld_cnt*towerFld_size + k-Nh; + towersData[towIndx] = hydroTauFlds[iFld*fldStride+ijk]; + towerFld_cnt += 1; + } + if(k == kMin){ + ij = i*(Nyp+2*Nh) + j; + towerSurfBaseAddress = towerCount*(NtBatch*towerSurfInstanceSize); + towIndx = towerSurfBaseAddress; + towersSurfData[towIndx] = z0m[ij]; + towIndx += 1; //only a single surface value so increment by 1 + towersSurfData[towIndx] = z0t[ij]; + towIndx += 1; //only a single surface value so increment by 1 + towersSurfData[towIndx] = tskin[ij]; + towIndx += 1; //only a single surface value so increment by 1 + towersSurfData[towIndx] = fricVel[ij]; + towIndx += 1; //only a single surface value so increment by 1 + towersSurfData[towIndx] = invOblen[ij]; + towIndx += 1; //only a single surface value so increment by 1 + towersSurfData[towIndx] = htFlux[ij]; + towIndx += 1; //only a single surface value so increment by 1 + if(moistureNvars*moistureSelector > 0){ + towersSurfData[towIndx] = qskin[ij]; + towIndx += 1; //only a single surface value so increment by 1 + towersSurfData[towIndx] = qFlux[ij]; + towIndx += 1; //only a single surface value so increment by 1 + }//end if Nmoist > 0 + } + }// end for k + } + }// end if rank_nTowers > 0 + + return(errorCode); +} //end hydro_coreAllocateProfilesDataStructure() + /*----->>>>> int hydro_coreSetupBndyPlanesAllRanks(); --------------------------------------------------- * Utility to read/scatter (across ranks as appropriate) the next set of BdyPlanes in the series */ @@ -1971,19 +2079,17 @@ int hydro_coreSetupBndyPlanesAllRanks(){ sprintf(fieldName,"theta"); fieldIndex = 4; errorCode = hydro_coreReadFieldBndyPlanes(ncid, fieldName, fieldIndex); + sprintf(fieldName,"TKE_0"); + fieldIndex = 5; + errorCode = hydro_coreReadFieldBndyPlanes(ncid, fieldName, fieldIndex); if(moistureSelector > 0){ if(moistureNvars > 0){ sprintf(fieldName,"qv"); - fieldIndex = 5; + fieldIndex = 6; errorCode = hydro_coreReadFieldBndyPlanes(ncid, fieldName, fieldIndex); } if(moistureNvars > 1){ sprintf(fieldName,"ql"); - fieldIndex = 6; - errorCode = hydro_coreReadFieldBndyPlanes(ncid, fieldName, fieldIndex); - } - if(moistureNvars > 2){ - sprintf(fieldName,"qr"); fieldIndex = 7; errorCode = hydro_coreReadFieldBndyPlanes(ncid, fieldName, fieldIndex); } @@ -3289,6 +3395,21 @@ int hydro_coreCleanup(){ memReleaseFloat(hydroAuxScalarsFrhs); } //end if NhydroAuxScalars + if(rank_nTowers > 0){ + free(towerIDs); + free(tower_iInds); + free(tower_jInds); + if(towerProfiles.coordType == 0){ + free(tower_LonOffsets); + free(tower_LatOffsets); + }else{ + free(tower_xOffsets); + free(tower_yOffsets); + }//end if coordType == 0, else... + free(towersData); + free(towersSurfData); + } + #ifdef GAD_EXT if(GADSelector > 0){ errorCode = GADCleanup(); @@ -3429,14 +3550,16 @@ int hydro_coreAddFieldAttributes(char *fieldName, int isForcing) { {"ql", "g kg-1", "Cloud liquid water mixing ratio", "cloud_liquid_water_mixing_ratio"}, {"fricVel", "m s-1", "Surface friction velocity", "surface_friction_velocity"}, {"htFlux", "K m s-1", "Surface sensible heat flux", "surface_upward_sensible_heat_flux"}, - {"qFlux", "kg kg-1 m s-1", "Surface latent heat flux", "surface_upward_latent_heat_flux"}, + {"qFlux", "g kg-1 m s-1", "Surface latent heat flux", "surface_upward_latent_heat_flux"}, {"tskin", "K", "Surface skin temperature", "surface_temperature"}, - {"qskin", "kg kg-1", "Surface skin water vapor mixing ratio", NULL}, + {"qskin", "g kg-1", "Surface skin water vapor mixing ratio", NULL}, {"z0m", "m", "Roughness length for momentum", "surface_roughness_length_for_momentum_in_air"}, {"z0t", "m", "Roughness length for heat", "surface_roughness_length_for_heat_in_air"}, - {"invOblen", "m-1", "Inverse Obukhov length", "atmosphere_boundary_layer_thickness"}, + {"invOblen", "m-1", "Inverse Obukhov length", NULL}, {"CanopyLAD", "m-1", "Leaf area density", "leaf_area_density"}, {"SeaMask", "-", "Sea mask", "sea_area_fraction"}, + {"lat", "degree_north", "Latitude", "latitude"}, + {"lon", "degree_east", "Longitude", "longitude"}, {NULL, NULL, NULL, NULL} // End marker }; diff --git a/SRC/HYDRO_CORE/hydro_core.h b/SRC/HYDRO_CORE/hydro_core.h index 38214e37..bc765adf 100644 --- a/SRC/HYDRO_CORE/hydro_core.h +++ b/SRC/HYDRO_CORE/hydro_core.h @@ -285,6 +285,20 @@ extern float thetaPerturbationAmplitude; /* Initial theta perturbations maximum extern int physics_oneRKonly; /* selector to apply physics RHS forcing only at the latest RK stage */ +/*---HYDRO_CORE_TOWERS*/ +extern int *towerIDs; +extern int *tower_iInds; +extern int *tower_jInds; +extern int rank_nTowers; +extern float *tower_xOffsets; +extern float *tower_yOffsets; +extern double *tower_LonOffsets; +extern double *tower_LatOffsets; +extern int towerInstanceSize; +extern int towerSurfInstanceSize; +extern float *towersData; +extern float *towersSurfData; + /*###################------------- HYDRO_CORE module function declarations ---------------------#################*/ /*----->>>>> int hydro_coreGetParams(); ---------------------------------------------------------------------- @@ -319,6 +333,11 @@ int hydro_coreGetFieldName(char * fldName, int iFld); */ int hydro_coreSetBaseState(); +/*----->>>>> int hydro_coreAllocateProfilesDataStructure(); --------------------------------------------------- +* Utility to allocate virtual tower data structures on appropriate ranks +*/ +int hydro_coreAllocateTowersDataStructure(int nProfs, ioProfiles_t towProfs, int NtBatch); + /*----->>>>> int hydro_coreSetupBndyPlanesAllRanks(); --------------------------------------------------- * Utility to read/scatter (across ranks as appropriate) the next set of BdyPlanes in the series */ @@ -340,7 +359,7 @@ int hydro_coreScatterFieldBndyPlanes(int Nfields); */ int hydro_coreReadNextBndyPlanesFile(); -/*----->>>>> int hydroi_coreTVCP(); ----------------------------------------------------------- +/*----->>>>> int hydro_coreTVCP(); ----------------------------------------------------------- * Updates model parameters used by the CELLPERT submodule from dynamic lateral BNDY conditions. */ int hydro_coreTVCP(float dt); diff --git a/SRC/IO/io.c b/SRC/IO/io.c index 76b3ebf3..98ebc061 100644 --- a/SRC/IO/io.c +++ b/SRC/IO/io.c @@ -32,6 +32,13 @@ size_t start2d[MAXDIMS]; size_t count2dTD[MAXDIMS]; size_t start2dTD[MAXDIMS]; + +// Include the netCDF-centric source code +#include + +// Include the unformatted N-to-N binary-centric source code +#include + /*######################------------------- IO module variable definitions ---------------------#################*/ /* Parameters */ int ioOutputMode; /*0: N-to-1 gather and write to a netCDF file, 1: N-to-N writes of FastEddy binary files*/ @@ -45,6 +52,9 @@ int frqOutput; /*frequency (in timesteps) at which to produce output; should char *outSubString; /*subString portion of outFile holding element-in-series as in path/base.substring */ char *outFileName; /*full name instance of outFileName = path/base.substring */ char *inFileName; /*full name instance of inFileName = path/infile */ +int registeredVars; /* Total number of variables registered with the primary ioVarsList */ +int registered3dVars; /* Number of 3-dimensional variables registered with the primary ioVarsList */ +int registered2dVars; /* Number of 2-dimensional variables registered with the primary ioVarsList */ /*IO-Buffers*/ float *ioBuffField; @@ -56,6 +66,14 @@ int *ioBuffFieldInt; int nz_varid; int ny_varid; int nx_varid; + +/*IO-profiles*/ +int towerIOSelector; +char *towerSpecsFile; +char *towerPath; /* Directory Path where tower files are to be written */ +int nProfs; +ioProfiles_t towerProfiles; + /*######################------------------- IO module function definitions ---------------------#################*/ /*----->>>>> int ioGetParams(); ---------------------------------------------------------------------- @@ -74,6 +92,12 @@ int ioGetParams(){ errorCode = queryPathParameter("outPath", &outPath, PARAM_MANDATORY); errorCode = queryStringParameter("outFileBase", &outFileBase, PARAM_MANDATORY); errorCode = queryIntegerParameter("frqOutput", &frqOutput, 0, INT_MAX, PARAM_MANDATORY); + towerIOSelector = 0; // Default to 0 + errorCode = queryIntegerParameter("towerIOSelector", &towerIOSelector, 0, 1, PARAM_OPTIONAL); + if(towerIOSelector > 0){ + errorCode = queryFileParameter("towerSpecsFile", &towerSpecsFile, PARAM_MANDATORY); + errorCode = queryPathParameter("towerPath", &towerPath, PARAM_MANDATORY); + } return(errorCode); } //end ioGetParams() @@ -92,6 +116,11 @@ int ioInit(){ printParameter("outPath", "Path where output files are to be written"); printParameter("outFileBase", "Base name of the output file series as in (outFileBase).element-in-series"); printParameter("frqOutput", "frequency (in timesteps) at which to produce output; should be an even multiple of NtBatch"); + printParameter("towerIOSelector", "Virtual Tower IO-Selector: 0=off, 1=on "); + if(towerIOSelector > 0){ + printParameter("towerSpecsFile", "netCDF file with virtual tower IO specifications "); + printParameter("towerPath", "Path where tower files are to be written"); + } } //end if(mpi_rank_world == 0) /*Broadcast the parameters across mpi_ranks*/ @@ -166,13 +195,46 @@ int ioInit(){ MPI_Bcast(outFileBase, strLength, MPI_CHARACTER, 0, MPI_COMM_WORLD); /*frqOutput*/ MPI_Bcast(&frqOutput, 1, MPI_INTEGER, 0, MPI_COMM_WORLD); + /*towerIOSelector*/ + MPI_Bcast(&towerIOSelector, 1, MPI_INT, 0, MPI_COMM_WORLD); + if(towerIOSelector > 0){ + strLength = 0; + if(mpi_rank_world == 0){ + if(towerSpecsFile != NULL){ + strLength = strlen(towerSpecsFile)+1; + }else{ + strLength = 0; + } + } //end if(mpi_rank_world == 0) + MPI_Bcast(&strLength, 1, MPI_INTEGER, 0, MPI_COMM_WORLD); + if(strLength > 0){ + if(mpi_rank_world != 0){ + towerSpecsFile = (char *) malloc(strLength*sizeof(char)); + } //if a non-root mpi_rank + MPI_Bcast(towerSpecsFile, strLength, MPI_CHARACTER, 0, MPI_COMM_WORLD); + } + /*towerPath string*/ + strLength = 0; + if(mpi_rank_world == 0){ + if(towerPath != NULL){ + strLength = strlen(towerPath)+1; + }else{ + strLength = 0; + } + } //end if(mpi_rank_world == 0) + MPI_Bcast(&strLength, 1, MPI_INTEGER, 0, MPI_COMM_WORLD); + if(mpi_rank_world != 0){ + towerPath = (char *) malloc(strLength*sizeof(char)); + } //if a non-root mpi_rank + MPI_Bcast(towerPath, strLength, MPI_CHARACTER, 0, MPI_COMM_WORLD); + } /*end-- Broadcast the parameters... */ /*Allocate IO private arrays*/ inFileName = (char *) malloc(3*MAX_LEN*sizeof(char)); /* 3 for each part of path/base.subString */ outSubString = (char *) malloc(MAX_LEN*sizeof(char)); outFileName = (char *) malloc(3*MAX_LEN*sizeof(char)); /* 3 for each part of path/base.subString */ - + return(errorCode); } //end ioInit() @@ -196,11 +258,161 @@ int ioAllocateBuffers(int globalNx, int globalNy, int globalNz){ return(errorCode); } //end ioAllocateBuffers() -// Include the netCDF-centric source code -#include +/*----->>>>> int ioProfilePreparations(); ------------------------------------------------------------ + * Profile Preparations routine for the IO module. Includes counting registered variables, reading + * profiles, allocating associated arrays, + * and... + */ +int ioProfilePreparations(){ + int errorCode = IO_SUCCESS; + int ncid; + int ncfldid; + int dimids[64]; + size_t count[64]; + size_t start[64]; + char fldName[64]; -// Include the unformatted N-to-N binary-centric source code -#include + int iprofile; + + //Count the total number of variables registered with IO + registeredVars = countVarsInList(®istered3dVars,®istered2dVars); + printf("mpi_rank_world--%d/%d ioProfilePreparations(): There are %d total registered variables in the primary ioVarsList.\n", + mpi_rank_world,mpi_size_world,registeredVars); + printf("mpi_rank_world--%d/%d ioProfilePreparations(): There are %d registered 3-dimensional variables in the primary ioVarsList.\n", + mpi_rank_world,mpi_size_world,registered3dVars); + printf("mpi_rank_world--%d/%d ioProfilePreparations(): There are %d registered 2-dimensional variables in the primary ioVarsList.\n", + mpi_rank_world,mpi_size_world,registered2dVars); + fflush(stdout); +#if 0 + if(mpi_rank_world == 0){ + errorCode = printList(); + } + fflush(stdout); +#endif + + /* ------------------------ Read in the towerSpecsFile -----------------------------*/ + //Root-rank should read the netcdf towerSpecsFile + if(mpi_rank_world == 0){ + //Open the netcdf towerSpecsFile + errorCode = ioOpenNetCDFinFile(towerSpecsFile, &ncid); + if(errorCode > 0){ + printf("Failed to open towerSpecsFile = %s EXITING NOW!!!!\n",towerSpecsFile); + fflush(stdout); + exit(0); + } + //Inquire for the dimID of the towers input specification fundamental parameter, nProfs + if((errorCode = nc_inq_dimid(ncid, "nProfs", &dimids[0]))){ + ERR(errorCode); + } + //Inquire for the value of nProfs + if((errorCode = nc_inq_dimlen(ncid, dimids[0], &count[dimids[0]]))){ + ERR(errorCode); + } + //Assign the nProfs to the value of the netCDF file dimension + nProfs = (int)count[dimids[0]]; + }//end if mpi_rank_world == 0 + MPI_Bcast(&nProfs, 1, MPI_INTEGER, 0, MPI_COMM_WORLD); + if(nProfs > 0){ + towerProfiles.profIDs = malloc(nProfs*sizeof(int)); + towerProfiles.mpi_ranks = malloc(nProfs*sizeof(int)); + towerProfiles.coordsSN = malloc(nProfs*sizeof(float)); + towerProfiles.coordsWE = malloc(nProfs*sizeof(float)); + towerProfiles.coordsLat = malloc(nProfs*sizeof(double)); + towerProfiles.coordsLon = malloc(nProfs*sizeof(double)); + } + if(mpi_rank_world == 0){ + //Setup the profIDs (these are common profile indices, shared across all ranks) + for(iprofile = 0; iprofile < nProfs; iprofile++){ + towerProfiles.profIDs[iprofile] = iprofile; + } + start[0] = 0; + sprintf(fldName,"coordType"); + if ( (errorCode = nc_inq_varid(ncid, fldName, &ncfldid)) ){ + ERR(errorCode); + } //if nc_inq_varid + if ((errorCode = nc_get_var_int(ncid, ncfldid, &towerProfiles.coordType )) ){ + ERR(errorCode); + } + if(towerProfiles.coordType == 0){ + sprintf(fldName,"coordsLat"); + if ( (errorCode = nc_inq_varid(ncid, fldName, &ncfldid)) ){ + ERR(errorCode); + } //if nc_inq_varid + if ((errorCode = nc_get_vara_double(ncid, ncfldid, &start[0], &count[dimids[0]], towerProfiles.coordsLat )) ){ + ERR(errorCode); + } + sprintf(fldName,"coordsLon"); + if ( (errorCode = nc_inq_varid(ncid, fldName, &ncfldid)) ){ + ERR(errorCode); + } //if nc_inq_varid + if ((errorCode = nc_get_vara_double(ncid, ncfldid, &start[0], &count[dimids[0]], towerProfiles.coordsLon )) ){ + ERR(errorCode); + } + }else{ + sprintf(fldName,"coordsSN"); + if ( (errorCode = nc_inq_varid(ncid, fldName, &ncfldid)) ){ + ERR(errorCode); + } //if nc_inq_varid + if ((errorCode = nc_get_vara_float(ncid, ncfldid, &start[0], &count[dimids[0]], towerProfiles.coordsSN )) ){ + ERR(errorCode); + } + sprintf(fldName,"coordsWE"); + if ( (errorCode = nc_inq_varid(ncid, fldName, &ncfldid)) ){ + ERR(errorCode); + } //if nc_inq_varid + if ((errorCode = nc_get_vara_float(ncid, ncfldid, &start[0], &count[dimids[0]], towerProfiles.coordsWE )) ){ + ERR(errorCode); + } + }//end if coordType == 0, else + } //end if mpi_rank_world == 0 + MPI_Bcast(towerProfiles.profIDs, nProfs, MPI_INTEGER, 0, MPI_COMM_WORLD); + MPI_Bcast(&towerProfiles.coordType, 1, MPI_INTEGER, 0, MPI_COMM_WORLD); + + if(towerProfiles.coordType == 0){ + MPI_Bcast(towerProfiles.coordsLat, nProfs, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(towerProfiles.coordsLon, nProfs, MPI_DOUBLE, 0, MPI_COMM_WORLD); + for(iprofile = 0; iprofile< nProfs; iprofile++){ + printf("mpi_rank_world--%d/%d: profIDs[%d] = %d-- coordType = %d, coords[%d](lat,lon) = (%f,%f)\n", + mpi_rank_world,mpi_size_world, + iprofile,towerProfiles.profIDs[iprofile], + towerProfiles.coordType, + iprofile,towerProfiles.coordsLat[iprofile],towerProfiles.coordsLon[iprofile]); + } + }else{ + MPI_Bcast(towerProfiles.coordsSN, nProfs, MPI_FLOAT, 0, MPI_COMM_WORLD); + MPI_Bcast(towerProfiles.coordsWE, nProfs, MPI_FLOAT, 0, MPI_COMM_WORLD); + for(iprofile = 0; iprofile< nProfs; iprofile++){ + printf("mpi_rank_world--%d/%d: profIDs[%d] = %d-- coordType = %d, coords[%d](y,x) = (%f,%f)\n", + mpi_rank_world,mpi_size_world, + iprofile,towerProfiles.profIDs[iprofile], + towerProfiles.coordType, + iprofile,towerProfiles.coordsSN[iprofile],towerProfiles.coordsWE[iprofile]); + } + }//end if coordType == 0, else + fflush(stdout); + return(errorCode); +} //end ioProfilePreparations() + +/*----->>>>> int ioCleanupProfiles(); ---------------------------------------------------------------------- +Used to free all malloced memory associated with Priofiles output from the IO module. +*/ +int ioCleanupProfiles(){ + int errorCode = IO_SUCCESS; + + if(towerIOSelector > 0 && nProfs > 0){ + //SOA + free(towerProfiles.profIDs); + free(towerProfiles.mpi_ranks); + free(towerProfiles.coordsSN); + free(towerProfiles.coordsWE); + free(towerProfiles.coordsLat); + free(towerProfiles.coordsLon); + free(towerSpecsFile); + free(towerPath); + } + return(errorCode); + +}//end ioCleanupProfiles() /*----->>>>> int ioCleanup(); ---------------------------------------------------------------------- Used to free all malloced memory by the IO module. @@ -208,6 +420,8 @@ Used to free all malloced memory by the IO module. int ioCleanup(){ int errorCode = IO_SUCCESS; + errorCode = ioCleanupProfiles(); + /*free the io-buffers*/ if(mpi_rank_world == 0){ free(ioBuffField); diff --git a/SRC/IO/io.h b/SRC/IO/io.h index 2f996517..acaf5fda 100644 --- a/SRC/IO/io.h +++ b/SRC/IO/io.h @@ -31,6 +31,16 @@ #include #include +typedef struct _ioProfiles_t { // Structure of Arrays (SOA) + int *profIDs; + int *mpi_ranks; + int coordType; //0 = lat,lon or 1 = y,x + double *coordsLon; + double *coordsLat; + float *coordsSN; + float *coordsWE; +} ioProfiles_t; + /*######################------------------- IO module variable declarations ---------------------#################*/ /* Parameters */ extern int ioOutputMode; /*0: N-to-1 gather and write to a netCDF file, 1: N-to-N writes of FastEddy binary files*/ @@ -44,6 +54,10 @@ extern int frqOutput; /*frequency in timesteps to produce output; should be extern char *outSubString; /*subString portion of outFile holding element-in-series as in path/base.substring */ extern char *outFileName; /*full name instance of outFileName = path/base.substring */ extern char *inFileName; /*full name instance of inFileName = path/infile */ +extern int registeredVars; /* Total number of variables registered with the primary ioVarsList */ +extern int registered3dVars; /* Number of 3-dimensional variables registered with the primary ioVarsList */ +extern int registered2dVars; /* Number of 2-dimensional variables registered with the primary ioVarsList */ + extern int nz_varid; extern int ny_varid; extern int nx_varid; @@ -53,6 +67,14 @@ extern float *ioBuffField; extern float *ioBuffFieldTransposed; extern float *ioBuffFieldRho; extern float *ioBuffFieldTransposed2D; +extern int *ioBuffFieldInt; + +/*IO-profiles*/ +extern int towerIOSelector; +extern char *towerSpecsFile; /* The path+filename to a tower specifications file*/ +extern char *towerPath; /* Directory Path where tower files are to be written */ +extern int nProfs; +extern ioProfiles_t towerProfiles; /*######################------------------- IO module function declarations ---------------------#################*/ @@ -71,6 +93,19 @@ int ioInit(); */ int ioAllocateBuffers(int globalNx, int globalNy, int globalNz); +/*----->>>>> int ioProfilePreparations(); ------------------------------------------------------------ + * Profile Preparations routine for the IO module. Includes counting registered variables, reading + * profiles, slices/planes or other reduced-data-volume output specifications, allocating associated arrays, + * and... + */ +//int ioProfilePreparations(int *numProfiles); +int ioProfilePreparations(); + +/*----->>>>> int ioCleanupProfiles(); ---------------------------------------------------------------------- +Used to free all malloced memory associated with Priofiles output from the IO module. +*/ +int ioCleanupProfiles(); + /*----->>>>> int ioCleanup(); ---------------------------------------------------------------------- Used to free all malloced memory by the IO module. */ diff --git a/SRC/IO/ioVarsList.c b/SRC/IO/ioVarsList.c index 94f8edd3..ab4fafc7 100644 --- a/SRC/IO/ioVarsList.c +++ b/SRC/IO/ioVarsList.c @@ -130,7 +130,7 @@ int addStandardAttrsToVar(char *varName, char *units, char *longName, char *stan } //end addStandardAttrsToVar int printList(){ - int i, j; + int i, j; ioVar_t *tmp; /*print the contents of the list from beginning to end*/ i = 0; @@ -175,6 +175,30 @@ int printList(){ return(i); } //end printList +int countVarsInList(int *reg3dVars, int *reg2dVars){ + int i; + int i3d; + int i2d; + ioVar_t *tmp; + /*Count the registered variables in the list from beginning to end*/ + i = 0; + i3d = 0; + i2d = 0; + tmp = head; + while(tmp != NULL){ + if (tmp->nDims == 4){ + i3d++; + }else if(tmp->nDims == 3){ + i2d++; + } + tmp = tmp->next; + i++; + }// end while + *reg3dVars = i3d; + *reg2dVars = i2d; + return(i); +} //end countVarsInList + void destroyList(){ ioVar_t *tmp; diff --git a/SRC/IO/ioVarsList.h b/SRC/IO/ioVarsList.h index a72b4579..9d5bd178 100644 --- a/SRC/IO/ioVarsList.h +++ b/SRC/IO/ioVarsList.h @@ -50,6 +50,7 @@ ioVar_t *getFirstVarFromList(); ioVar_t *getNamedVarFromList(char* name); int addVarToList(char *name, char *type, int nDims, int *dimids, void *varMemAddress); int printList(); +int countVarsInList(int *reg3dVars, int *reg2dVars); void destroyList(); /* Add a single NetCDF attribute to an existing variable in the list diff --git a/SRC/IO/io_binary.c b/SRC/IO/io_binary.c index 442c18c8..fc6e2dbb 100644 --- a/SRC/IO/io_binary.c +++ b/SRC/IO/io_binary.c @@ -184,3 +184,130 @@ int ioPutBinaryoutFileVars(FILE *outptr, int Nx, int Ny, int Nz, int Nh){ return(errorCode); } //ioPutBinaryoutFileVars() + +/*----->>>>> int ioWriteBinaryTowerFileSingleBatch(); --------------------------------------------------------------- + * Used to have N-ranks write N-binary files of virtual tower data structures for a batch of timesteps. + */ +int ioWriteBinaryTowerFileSingleBatch(int tstep, int batchSize, int Nz, float *batchTimes, float *towersData, float *towersSurfData, + int *towerIDs, int rank_nTowers, int towerInstanceSize, int towerSurfInstanceSize){ + int errorCode = IO_SUCCESS; + FILE *output_ptr; + char towerSubString[64]; + char towerFileName[256]; + int itower; + + for(itower = 0; itower < rank_nTowers; itower++){ + //--------------Tower profile-variables file + /* build the subString tag */ + sprintf(towerSubString, "tower_%d.%d",towerIDs[itower],tstep-batchSize); + /* concatenate the fileName components */ + sprintf(towerFileName, "%s%s",towerPath,towerSubString); + /*Open the output file*/ + output_ptr = fopen(towerFileName,"wb"); + /*Write the batch of tower instances to the output file*/ + fwrite(&towerInstanceSize,sizeof(int),1,output_ptr); + fwrite(&batchSize,sizeof(int),1,output_ptr); + fwrite(&batchTimes[0],batchSize*sizeof(float),1,output_ptr); + fwrite(&towersData[itower*batchSize*towerInstanceSize],batchSize*towerInstanceSize*sizeof(float),1,output_ptr); + /*Close the output file*/ + fclose(output_ptr); + + //--------------Tower surface-variables file + /* build the subString tag */ + sprintf(towerSubString, "tower_sv_%d.%d",towerIDs[itower],tstep-batchSize); + /* concatenate the fileName components */ + sprintf(towerFileName, "%s%s",towerPath,towerSubString); + /*Open the output file*/ + output_ptr = fopen(towerFileName,"wb"); + /*Write the batch of tower instances to the output file*/ + fwrite(&towerSurfInstanceSize,sizeof(int),1,output_ptr); + fwrite(&batchSize,sizeof(int),1,output_ptr); + fwrite(&batchTimes[0],batchSize*sizeof(float),1,output_ptr); + fwrite(&towersSurfData[itower*batchSize*towerSurfInstanceSize],batchSize*towerSurfInstanceSize*sizeof(float),1,output_ptr); + /*Close the output file*/ + fclose(output_ptr); + } //end for itower + + + return(errorCode); +} //end ioWriteBinaryTowerFileSingleBatch() + +/*----->>>>> int ioWriteBinaryTowerInitialFile(); --------------------------------------------------------------- + * Used to have N-ranks write N-binary files of virtual tower data structures for the initial timestep and time-independent fields. + */ +int ioWriteBinaryTowerInitialFile(float dt, int itStart, int Nx, int Ny, int Nz, int Nh, + float *towersData, float *towersSurfData, + int *towerIDs, int rank_nTowers, int *tower_iInds, int *tower_jInds, + int coordType, float *tower_xOffs, float *tower_yOffs, double *tower_LonOffs, double *tower_LatOffs, + int batchSize, int towerInstanceSize, int towerSurfInstanceSize, + float *zCoords, float *yCoords, float *xCoords, float *topoFld, int surflayer_offshore, float *seamask){ + int errorCode = IO_SUCCESS; + FILE *output_ptr; + char towerSubString[64]; + char towerFileName[256]; + int itower; + int i,j,k,ijk,ij; + int iStride,jStride,kStride; + int tmpOne = 1; + float timeStart; + double tmpDbleWE; + double tmpDbleSN; + + iStride = (Ny+2*Nh)*(Nz+2*Nh); + jStride = (Nz+2*Nh); + kStride = 1; + + timeStart = itStart*dt; + + for(itower = 0; itower < rank_nTowers; itower++){ + i = tower_iInds[itower]; + j = tower_jInds[itower]; + k = Nh; + ijk = i*iStride + j*jStride + k*kStride; + ij = i*(Ny+2*Nh) + j; + //--------------Tower profile-variables + /* build the subString tag */ + sprintf(towerSubString, "tower_ic_%d.%d",towerIDs[itower],itStart); + /* concatenate the fileName components */ + sprintf(towerFileName, "%s%s",towerPath,towerSubString); + /*Open the output file*/ + output_ptr = fopen(towerFileName,"wb"); + + //--------------Tower time-independent variables + fwrite(&Nz,sizeof(int),1,output_ptr); + fwrite(&zCoords[ijk],Nz*sizeof(float),1,output_ptr); + fwrite(&yCoords[ijk],sizeof(float),1,output_ptr); + fwrite(&xCoords[ijk],sizeof(float),1,output_ptr); + fwrite(&topoFld[ij],sizeof(float),1,output_ptr); + if(surflayer_offshore > 0){ + fwrite(&seamask[ij],sizeof(float),1,output_ptr); + }//end if surfacelayer_offshore + if(coordType == 0){ + tmpDbleSN = tower_LatOffs[itower]; + tmpDbleWE = tower_LonOffs[itower]; + }else{ + tmpDbleSN = (double) tower_yOffs[itower]; + tmpDbleWE = (double) tower_xOffs[itower]; + } + fwrite(&tmpDbleSN,sizeof(double),1,output_ptr); + fwrite(&tmpDbleWE,sizeof(double),1,output_ptr); + + /*Write the batch of tower instances to the output file*/ + fwrite(&towerInstanceSize,sizeof(int),1,output_ptr); + fwrite(&tmpOne,sizeof(int),1,output_ptr); + fwrite(&timeStart,sizeof(float),1,output_ptr); + fwrite(&towersData[itower*batchSize*towerInstanceSize],1*towerInstanceSize*sizeof(float),1,output_ptr); + + //--------------Tower surface-variables + fwrite(&towerSurfInstanceSize,sizeof(int),1,output_ptr); + fwrite(&tmpOne,sizeof(int),1,output_ptr); + fwrite(&timeStart,sizeof(float),1,output_ptr); + fwrite(&towersSurfData[itower*batchSize*towerSurfInstanceSize],1*towerSurfInstanceSize*sizeof(float),1,output_ptr); + + /*Close the output file*/ + fclose(output_ptr); + } //end for itower + + + return(errorCode); +} //end ioWriteBinaryTowerInitialFile() diff --git a/SRC/IO/io_binary.h b/SRC/IO/io_binary.h index 02528aa4..8e389d4a 100644 --- a/SRC/IO/io_binary.h +++ b/SRC/IO/io_binary.h @@ -16,7 +16,7 @@ //////////*********************** INPUT FUNCTIONS *********************************//////// //////////*********************** OUTPUT FUNCTIONS *********************************//////// -/*----->>>>> int ioWriteiBinaryoutFileSingleTime(); --------------------------------------------------------------- +/*----->>>>> int ioWriteBinaryoutFileSingleTime(); --------------------------------------------------------------- * Used to have N-ranks write N-binary files of registered variables for a single timestep. */ #ifdef GAD_EXT @@ -32,3 +32,18 @@ int ioPutBinaryoutFileVars(FILE *outptr, int Nx, int Ny, int Nz, int Nh, int Ntu #else int ioPutBinaryoutFileVars(FILE *outptr, int Nx, int Ny, int Nz, int Nh); #endif +/*----->>>>> int ioWriteBinaryTowerFileSingleBatch(); --------------------------------------------------------------- + * Used to have N-ranks write N-binary files of virtual tower data structures for a batch of timesteps. + */ +int ioWriteBinaryTowerFileSingleBatch(int tstep, int batchSize, int Nz, float *batchTimes, float *towersData, float *towersSurfData, + int *towerIDs, int rank_nTowers, int towerInstanceSize, int towerSurfInstanceSize); + +/*----->>>>> int ioWriteBinaryTowerInitialFile(); --------------------------------------------------------------- + * Used to have N-ranks write N-binary files of virtual tower data structures for the initial timestep and time-independent fields. + */ +int ioWriteBinaryTowerInitialFile(float dt, int itStart, int Nx, int Ny, int Nz, int Nh, + float *towersData, float *towersSurfData, + int *towerIDs, int rank_nTowers, int *tower_iInds, int *tower_jInds, + int coordType, float *tower_xOffs, float *tower_yOffs, double *tower_LonOffs, double *tower_LatOffs, + int batchSize, int towerInstanceSize, int towerSurfInstanceSize, + float *zCoords, float *yCoords, float *xCoords, float *topoFld, int surflayer_offshore, float *seamask); diff --git a/SRC/TIME_INTEGRATION/CUDA/cuda_timeIntDevice.cu b/SRC/TIME_INTEGRATION/CUDA/cuda_timeIntDevice.cu index c91c4732..897c25d3 100644 --- a/SRC/TIME_INTEGRATION/CUDA/cuda_timeIntDevice.cu +++ b/SRC/TIME_INTEGRATION/CUDA/cuda_timeIntDevice.cu @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include @@ -76,8 +77,8 @@ extern "C" int cuda_timeIntDeviceSetup(){ gpuErrchk( cudaPeekAtLastError() ); /*Check for errors in the cudaMalloc calls*/ - //Ensure secondary time-integration dependent hydro_core parameters get initialized - errorCode = cuda_hydroCoreDeviceSecondaryStageSetup(dt); + //Ensure secondary time-integration dependent, device-level hydro_core submodules get initialized + errorCode = cuda_hydroCoreDeviceSecondaryStageSetup(dt, NtBatch); //Inital Host-to-Device field copies errorCode = cuda_hydroCoreInitFieldsDevice(); //Transfer initial/restart conditions to the device //printf("cuda_timeIntDeviceSetup() complete.\n"); @@ -106,6 +107,7 @@ extern "C" int cuda_timeIntDeviceCommence(int it){ int errorCode = TIME_INTEGRATION_SUCCESS; int itBatch; int RKstage; + int Ntaus; #ifdef TIMERS_LEVEL1 float elapsedTime; cudaEvent_t startE, stopE @@ -153,14 +155,32 @@ extern "C" int cuda_timeIntDeviceCommence(int it){ stopSynchReportDestroyEvent(&startE, &stopE, &elapsedTime); printf("cuda_timeIntCommenceRK3_WS2002() Kernel execution time (ms): %12.8f\n", elapsedTime); #endif - } //end for RKstage + } //end for RKstage + if(towerIOSelector > 0){ + Ntaus = 9; + //Append time-advanced data to tower buffers + cudaDevice_towerAppendBuffers<<>>(itBatch, NtBatch, + towersData_d, towersSurfData_d, tower_iInds_d, tower_jInds_d, + Nhydro, hydroFlds_d, + TKESelector*turbulenceSelector, sgstkeScalars_d, + moistureNvars*moistureSelector, moistScalars_d, + NhydroAuxScalars, hydroAuxScalars_d, + hydroSubGridWrite*Ntaus, hydroTauFlds_d, + hydroSubGridWrite*moistureNvars*3, moistTauFlds_d, + z0m_d, z0t_d, tskin_d, qskin_d, + fricVel_d, invOblen_d, htFlux_d, qFlux_d); + + } } //end if(timeMethod == 0){... simTime_it = simTime_it + 1; //Increment the master simulation time step simTime = simTime_it * dt; /*Increment the master simulation time*/ + simTimeBatch[itBatch] = simTime; /*Store the master simulation time*/ }//end for itBatch... + gpuErrchk( cudaDeviceSynchronize() ); //Retrieve desired HYDRO_CORE fields from device - errorCode = cuda_hydroCoreSynchFieldsFromDevice(); + errorCode = cuda_hydroCoreSynchFieldsFromDevice(NtBatch); + gpuErrchk( cudaDeviceSynchronize() ); return(errorCode); }//end cuda_timeIntDeviceCommence() diff --git a/SRC/TIME_INTEGRATION/time_integration.c b/SRC/TIME_INTEGRATION/time_integration.c index 0544b3d7..14dbee17 100644 --- a/SRC/TIME_INTEGRATION/time_integration.c +++ b/SRC/TIME_INTEGRATION/time_integration.c @@ -41,6 +41,7 @@ int simTime_itRestart; /*Master simulation 'Restart' time step*/ int numRKstages; /* number of stages in the time scheme */ /* array fields */ +float* simTimeBatch; /*Array of master simulation time over NtBatch timesteps*/ /*######################------------------- TIME_INTEGRATION module function definitions ---------------------#################*/ @@ -125,6 +126,9 @@ int timeInit(){ numRKstages = 2; } + //Allocate space for an array of NtBatch master simulation time values + simTimeBatch = malloc(NtBatch*sizeof(float)); + /* Done */ return(errorCode); } //end timeInit() @@ -171,7 +175,7 @@ int timeCleanup(){ int errorCode = TIME_INTEGRATION_SUCCESS; /* Free any TIME_INTEGRATION module arrays */ - //currently none + free(simTimeBatch); return(errorCode); diff --git a/SRC/TIME_INTEGRATION/time_integration.h b/SRC/TIME_INTEGRATION/time_integration.h index fa79f746..d18dd4ae 100644 --- a/SRC/TIME_INTEGRATION/time_integration.h +++ b/SRC/TIME_INTEGRATION/time_integration.h @@ -33,6 +33,7 @@ extern int simTime_itRestart; /*Master simulation 'Restart' time step*/ extern int numRKstages; /* number of stages in the time scheme */ /* array fields */ +extern float* simTimeBatch; /*Array of master simulation time over NtBatch timesteps*/ /*############------------------- TIME_INTEGRATION module function declarations ---------------------############*/ @@ -63,4 +64,4 @@ Used to free all malloced memory by the TIME_INTEGRATION module. */ int timeCleanup(); -#endif // _TIME_INTEGRATION_H \ No newline at end of file +#endif // _TIME_INTEGRATION_H diff --git a/scripts/batch_jobs/run_mpassit.sh b/scripts/batch_jobs/run_mpassit.sh new file mode 100755 index 00000000..24645997 --- /dev/null +++ b/scripts/batch_jobs/run_mpassit.sh @@ -0,0 +1,123 @@ +#!/bin/bash + +#PBS -S /bin/csh +#PBS -N mpassit +#PBS -A P48503002 +#PBS -l walltime=50:00 +#PBS -q main +#PBS -o mpassit.out +#PBS -j oe +#PBS -k oed +#PBS -l select=8:ncpus=32:mpiprocs=32 +#PBS -l job_priority=premium +#PBS -m n +#PBS -V + +##SBATCH -J mpassit +##SBATCH -o logs/mpassit.%j +##SBATCH -e logs/mpassit.%j +##SBATCH -n 1200 +##SBATCH --exclusive +##SBATCH --partition=hera +##SBATCH -t 02:00:00 +##SBATCH -A hmtb + +# +# This script runs MPASSIT +# + +start_init=20240216170000 # YYYYMMDDHHMMSS format +diag_output_interval=300 #interval between lbc files in seconds +FCST_RANGE=5400 #length of forecast in seconds +MPAS_EXPT_DIR=/glade/derecho/scratch/wmayfield/dtc_ncar_mpas/expt_dirs/DTC_NCAR_hrrrIC/conus_3km_fasteddy/mpas_atm/2024021617/ens_1 #directory containing diag, history, init files +MPASSIT_CODE_DIR=/glade/campaign/ral/jntp/mayfield/fasteddy/MPASSIT #path to base MPASSIT code directory (contains ./bin/mpassit) +VARLIST_DIR=/glade/campaign/ral/jntp/mayfield/fasteddy/run_mpassit/varlists_mpassit_fasteddy #directory containing the variable lists +TOOL_DIR=/glade/u/home/schwartz/utils/derecho #directory with compiled WRFDA tools, in order to use "da_advance_time.exe" + +# Load modules: +module --force purge +module use ${MPASSIT_CODE_DIR}/modulefiles +module load build.derecho.intel + +#Start and End dates +DATE=`$TOOL_DIR/da_advance_time.exe ${start_init} 0 -f ccyymmddhhnnss` # +end_time=`$TOOL_DIR/da_advance_time.exe ${DATE} ${FCST_RANGE}s -f ccyymmddhhnnss` +echo "Starting init is $start_init" +echo "End time is $end_time" + +while [[ "$DATE" -le "$end_time" ]] ; do + + # ------------------------------------- + # Get current date into proper format + # ------------------------------------- + date_file_format=`${TOOL_DIR}/da_advance_time.exe $DATE 0 -f ccyy-mm-dd_hh.nn.ss` + echo "Date in mpas format is: ${date_file_format}" + vhr=`echo ${date_file_format}` + + # ---------------------------------- + # Make and go to working directory + # ---------------------------------- + rundir=./${vhr} + mkdir -p $rundir + cd $rundir + + #----------------------------------------------------------- + # Link necessary input files and code and fill namelist + #----------------------------------------------------------- + ln -sf ${MPASSIT_CODE_DIR}/bin/mpassit . + ln -sf ${VARLIST_DIR}/* . + + export grid_file=${MPAS_EXPT_DIR}/init.nc + export hist_file=${MPAS_EXPT_DIR}/history.${date_file_format}.nc + export diag_file=${MPAS_EXPT_DIR}/diag.${date_file_format}.nc + export output_file=../proc.${date_file_format}.nc + + #-------------------------------------------- + # Create the MPASSIT namelist.input + #-------------------------------------------- + + rm -f ./namelist.input + cat > ./namelist.input << EOF +&config +grid_file_input_grid="${grid_file}" +hist_file_input_grid="${hist_file}" +diag_file_input_grid="${diag_file}" +file_target_grid="/this/is/an/uneeded/path" +output_file="${output_file}" +target_grid_type = 'lambert' +interp_diag=.true. +interp_hist=.true. +wrf_mod_vars=.true. +esmf_log=.false. +nx = 1578 +ny = 925 +dx = 3000.0 +dy = 3000.0 +ref_lat = 38.4 +ref_lon = -97.0 +truelat1 = 38.4 +truelat2 = 38.4 +stand_lon = -97.0 / +EOF + + #---------------------------------------------------- + # Run MPASSIT + #---------------------------------------------------- + rm -f ./*.log + rm -f ./core* + rm -f ./*.err + + mpirun ./mpassit namelist.input + + if [[ "$status" -ne "0" ]] ; then + echo "MPASSIT failed. Exit." >> ./FAIL + exit 6 + fi + + # Done with this forecast hour; go to next one + cd .. + DATE=`$TOOL_DIR/da_advance_time.exe ${DATE} ${diag_output_interval}s -f ccyymmddhhnnss` + echo "Date is now ${DATE}" +done # loop over time/initializations + +exit 0 diff --git a/scripts/batch_jobs/varlists_mpassit_fasteddy/diaglist b/scripts/batch_jobs/varlists_mpassit_fasteddy/diaglist new file mode 100644 index 00000000..a2d173bb --- /dev/null +++ b/scripts/batch_jobs/varlists_mpassit_fasteddy/diaglist @@ -0,0 +1,3 @@ +q2 Q2 +z0 Z0 +znt ZNT diff --git a/scripts/batch_jobs/varlists_mpassit_fasteddy/histlist_2d b/scripts/batch_jobs/varlists_mpassit_fasteddy/histlist_2d new file mode 100644 index 00000000..290e8eca --- /dev/null +++ b/scripts/batch_jobs/varlists_mpassit_fasteddy/histlist_2d @@ -0,0 +1,2 @@ +surface_pressure PSFC +skintemp TSK diff --git a/scripts/batch_jobs/varlists_mpassit_fasteddy/histlist_3d b/scripts/batch_jobs/varlists_mpassit_fasteddy/histlist_3d new file mode 100644 index 00000000..21795ef2 --- /dev/null +++ b/scripts/batch_jobs/varlists_mpassit_fasteddy/histlist_3d @@ -0,0 +1,9 @@ +zgrid PHB +w W +theta T +uReconstructZonal U +uReconstructMeridional V +qv QVAPOR +qc QCLOUD +rho MUB +pressure P_HYD diff --git a/scripts/batch_jobs/varlists_mpassit_fasteddy/histlist_soil b/scripts/batch_jobs/varlists_mpassit_fasteddy/histlist_soil new file mode 100644 index 00000000..524768a6 --- /dev/null +++ b/scripts/batch_jobs/varlists_mpassit_fasteddy/histlist_soil @@ -0,0 +1 @@ +tslb TSLB diff --git a/scripts/python_utilities/coupler/GenICBCs.py b/scripts/python_utilities/coupler/GenICBCs.py index b2276507..de748fc6 100644 --- a/scripts/python_utilities/coupler/GenICBCs.py +++ b/scripts/python_utilities/coupler/GenICBCs.py @@ -25,7 +25,7 @@ mpi_name = MPI.Get_processor_name() print("{:d}/{:d}: Hello World! on {:s}.".format(mpi_rank, mpi_size, mpi_name)) -DEBUG_COUPLER = False #True +DEBUG_COUPLER = False ###################################################### ### Parse the command line arguments ### @@ -40,6 +40,7 @@ ICBC_dir = params["ICBC_dir"] FE_simGrid = params["FE_simGrid"] +parent_model = params["parent_model"] WRF_PrntDir = params["WRF_PrntDir"] WRF_PrntOutPrefix = params["WRF_PrntOutPrefix"] dateString = params["date0"] @@ -48,11 +49,23 @@ timeSecond0 = params["timeSecond0"] secMax = params["secMax"] secInc = params["secInc"] +FE_PrntDir = params["FE_PrntDir"] +FE_PrntOutPrefix = params["FE_PrntOutPrefix"] +itMin = params["itMin"] +dt_FE = params["dt_FE"] +outputFrequency = params["outputFrequency"] +timeLengthSec = params["timeLengthSec"] +nest_tke_opt = params["nest_tke_opt"] +ideal_opt = params["ideal_opt"] print(f"{mpi_rank}/{mpi_size}: Writing coupler outputs to {ICBC_dir}") print(f"{mpi_rank}/{mpi_size}: Interpolating to FE-domain from {FE_simGrid}") -print(f"{mpi_rank}/{mpi_size}: Processing of WRF-files {WRF_PrntDir}{WRF_PrntOutPrefix}*") -print(f"{mpi_rank}/{mpi_size}: Date and times of WRF-files to process: {dateString}_{timeHour0:02}:{timeMinute0:02}:*, every {secInc} s for {secMax} total seconds.") +if (parent_model == 0): + print(f"{mpi_rank}/{mpi_size}: Processing of WRF-files {WRF_PrntDir}{WRF_PrntOutPrefix}*") + print(f"{mpi_rank}/{mpi_size}: Date and times of WRF-files to process: {dateString}_{timeHour0:02}:{timeMinute0:02}:*, every {secInc} s for {secMax} total seconds.") +elif (parent_model == 1): + print(f"{mpi_rank}/{mpi_size}: Processing of FastEddy-files {FE_PrntDir}{FE_PrntOutPrefix}*") + print(f"{mpi_rank}/{mpi_size}: FastEddy-files to process: start at timestep {itMin}, every {outputFrequency} timesteps for {timeLengthSec} total seconds.") ################################################################################################ ### Create a coupler output directory for initial and boundary conditions if necessary @@ -67,19 +80,49 @@ files_list=[] times=[] -year0 = int(dateString[0:4]) -month0 = int(dateString[5:7]) -day0 = int(dateString[8:10]) - -date_it = dt.datetime(year0,month0,day0,timeHour0,timeMinute0,timeSecond0) -for it in range(0,secMax,secInc): - dateString_it = str(date_it.year) + '-' + "{:02d}".format(date_it.month) + '-' + "{:02d}".format(date_it.day) + '_' - thistime = "{:s}{:02d}:{:02d}:{:02d}".format(dateString_it,date_it.hour,date_it.minute,date_it.second) - file_tmp = f'{WRF_PrntDir}{WRF_PrntOutPrefix}{thistime}' - files_list.append(file_tmp) - if(mpi_rank == 0): - print(file_tmp) - date_it = date_it + dt.timedelta(seconds=secInc) +if (parent_model == 0): # WRF + + name_lat = 'XLAT' + name_lon = 'XLONG' + name_concat_dim = 'Time' + + year0 = int(dateString[0:4]) + month0 = int(dateString[5:7]) + day0 = int(dateString[8:10]) + + date_it = dt.datetime(year0,month0,day0,timeHour0,timeMinute0,timeSecond0) + for it in range(0,secMax,secInc): + dateString_it = str(date_it.year) + '-' + "{:02d}".format(date_it.month) + '-' + "{:02d}".format(date_it.day) + '_' + thistime = "{:s}{:02d}:{:02d}:{:02d}".format(dateString_it,date_it.hour,date_it.minute,date_it.second) + file_tmp = f'{WRF_PrntDir}{WRF_PrntOutPrefix}{thistime}' + files_list.append(file_tmp) + if(mpi_rank == 0): + print(file_tmp) + date_it = date_it + dt.timedelta(seconds=secInc) + +elif (parent_model == 1): # FastEddy + + name_lat = 'lat' + name_lon = 'lon' + name_concat_dim = 'time' + + itInc=np.int32(np.floor(outputFrequency/dt_FE)) + itMax = itMin + itInc*np.int32(np.floor(timeLengthSec/outputFrequency)) + print(f"Creating file list spanning timestep {itMin} to {itMax} in increments of {itInc} timesteps.") + timeMin = itMin*dt_FE + timeMax = itMax*dt_FE + timeInc = itInc*dt_FE + print(f"This corresponds to spanning time = {timeMin} [s] to {timeMax} [s] in increments of {timeInc} [s].") + + for it in range(itMin,itMax+(itInc-1),itInc): + thistime = f'{it}' + times.append(thistime) + # print(f"thistime={thistime}") + for eachtime in times: + file_tmp = f"{FE_PrntDir}{FE_PrntOutPrefix}.{eachtime}" + files_list.append(file_tmp) + if(mpi_rank == 0): + print(file_tmp) ################################################################################################ ### Setup mpi task decomposition over the set of files to process @@ -140,22 +183,64 @@ ## Find the the WRF d02 profiler locations itargs=[] jtargs=[] -#print('WRF:') -#print('(j,i):') - -latFE = ds_FEGrid.lat.values -lonFE = ds_FEGrid.lon.values - -corners_lat = np.asarray([latFE[0,0], latFE[0,-1],latFE[-1,-1], latFE[-1,0]]) -corners_lon = np.asarray([lonFE[0,0], lonFE[0,-1],lonFE[-1,-1], lonFE[-1,0]]) -print('corners_lat.shape=',corners_lat.shape) -print('corners_lon.shape=',corners_lon.shape) -print('corners_lat=',corners_lat) -print('corners_lon=',corners_lon) - -for indx in range(len(corners_lat)): - blah3=np.sqrt( (ds_WRFRef['XLAT'][0,:,:].values-corners_lat[indx])**2 - +(ds_WRFRef['XLONG'][0,:,:].values-corners_lon[indx])**2) + +if (parent_model == 0): + dx_parent = ds_WRFRef.attrs['DX'] + dy_parent = ds_WRFRef.attrs['DY'] + Ngpx = ds_WRFRef.sizes['west_east'] + Ngpy = ds_WRFRef.sizes['south_north'] +elif (parent_model == 1): + dx_parent = (ds_WRFRef['xPos'][0,0,0,1]-ds_WRFRef['xPos'][0,0,0,0]).values + dy_parent = (ds_WRFRef['yPos'][0,0,1,0]-ds_WRFRef['yPos'][0,0,0,0]).values + Ngpx = ds_WRFRef.sizes['xIndex'] + Ngpy = ds_WRFRef.sizes['yIndex'] + +if (not ideal_opt): + + latFE = ds_FEGrid.lat.values + lonFE = ds_FEGrid.lon.values + + corners_lat = np.asarray([latFE[0,0], latFE[0,-1],latFE[-1,-1], latFE[-1,0]]) + corners_lon = np.asarray([lonFE[0,0], lonFE[0,-1],lonFE[-1,-1], lonFE[-1,0]]) + print('corners_lat.shape=',corners_lat.shape) + print('corners_lon.shape=',corners_lon.shape) + print('corners_lat=',corners_lat) + print('corners_lon=',corners_lon) + + len_corners = len(corners_lat) + + corners_var_y = corners_lat + corners_var_x = corners_lon + +else: + + name_lat = 'yPos' + name_lon = 'xPos' + + xcoordFE = ds_FEGrid.xPos.isel(zIndex=0).values + ycoordFE = ds_FEGrid.yPos.isel(zIndex=0).values + print("xcoordFE.shape=",xcoordFE.shape) + + corners_x = np.asarray([xcoordFE[0,0], xcoordFE[0,-1], xcoordFE[-1,-1], xcoordFE[-1,0]]) + corners_y = np.asarray([ycoordFE[0,0], ycoordFE[0,-1], ycoordFE[-1,-1], ycoordFE[-1,0]]) + print('corners_x.shape=',corners_x.shape) + print('corners_y.shape=',corners_y.shape) + print('corners_x=',corners_x) + print('corners_y=',corners_y) + + len_corners = len(corners_y) + + corners_var_y = corners_x + corners_var_x = corners_y + +for indx in range(len_corners): + if (not ideal_opt): + blah3=np.sqrt( (ds_WRFRef[name_lat][0,:,:].values-corners_lat[indx])**2 + +(ds_WRFRef[name_lon][0,:,:].values-corners_lon[indx])**2) + else: + blah3=np.sqrt( (ds_WRFRef[name_lat][0,0,:,:].values-corners_y[indx])**2 + +(ds_WRFRef[name_lon][0,0,:,:].values-corners_x[indx])**2) + locCount=0 for jtarg, itarg in np.argwhere(blah3 == np.min(blah3,axis=(0,1))): if locCount < 1: @@ -166,74 +251,142 @@ else: #skip this redundant location of minimum distance if(mpi_rank == 0): - print('Skipping redundant closest corner location: ',jtarg,itarg) + print('Skipping redundant closest corner location: ',jtarg,itarg) #### Append the corner index pairs to the WRFref dataset dFE_jindxs and dFE_iindxs ds_WRFRef['dFE_jindxs']=xr.DataArray(np.asarray(jtargs,dtype=np.int32),dims=["corners"]) ds_WRFRef['dFE_iindxs']=xr.DataArray(np.asarray(itargs,dtype=np.int32),dims=["corners"]) -for indx in range(len(corners_lat)): - if(mpi_rank == 0): - print(f"corner({indx}) @ WRF({ds_WRFRef['dFE_jindxs'][indx].values},{ds_WRFRef['dFE_iindxs'][indx].values})") - print('[WRF,corner({:d})]: lats = [{:f},{:f}], lons = [{:f},{:f}]'.format(indx,ds_WRFRef['XLAT'][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]].values, - corners_lat[indx], - ds_WRFRef['XLONG'][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]].values, - corners_lon[indx])) - + +print("ds_WRFRef['dFE_jindxs'].values=",ds_WRFRef['dFE_jindxs'].values) +print("ds_WRFRef['dFE_iindxs'].values=",ds_WRFRef['dFE_iindxs'].values) + +if (not ideal_opt): + + for indx in range(len_corners): + if(mpi_rank == 0): + print(f"corner({indx}) @ WRF({ds_WRFRef['dFE_jindxs'][indx].values},{ds_WRFRef['dFE_iindxs'][indx].values})") + print('[WRF,corner({:d})]: lats = [{:f},{:f}], lons = [{:f},{:f}]'.format(indx,ds_WRFRef[name_lat][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]].values, + corners_lat[indx], + ds_WRFRef[name_lon][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]].values, + corners_lon[indx])) -### Compute FE-domain corner lat/lon offsets from closest dsWRFRef cell-centered lat/lons -print('\t Pre-correction offsets:') -for indx in range(len(corners_lat)): - latOff = ds_WRFRef['XLAT'][0,ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values]-corners_lat[indx] - lonOff = ds_WRFRef['XLONG'][0,ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values]-corners_lon[indx] - if(mpi_rank == 0): - print('\t corner({:d}): latOff = {:f}, lonOff = {:f}'.format(indx,latOff,lonOff)) - yoffset = -(ds_WRFRef['XLAT'][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]]-corners_lat[indx])\ - *(ds_WRFRef.attrs['DY']/( ds_WRFRef['XLAT'][0,ds_WRFRef['dFE_jindxs'][indx]+1,ds_WRFRef['dFE_iindxs'][indx]] - -ds_WRFRef['XLAT'][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]])) - xoffset = -(ds_WRFRef['XLONG'][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]]-corners_lon[indx])\ - *(ds_WRFRef.attrs['DX']/( ds_WRFRef['XLONG'][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]+1] - -ds_WRFRef['XLONG'][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]])) + ### Compute FE-domain corner lat/lon offsets from closest dsWRFRef cell-centered lat/lons + print('\t Pre-correction offsets:') + for indx in range(len_corners): + latOff = ds_WRFRef[name_lat][0,ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values]-corners_lat[indx] + lonOff = ds_WRFRef[name_lon][0,ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values]-corners_lon[indx] + if(mpi_rank == 0): + print('\t corner({:d}): latOff = {:f}, lonOff = {:f}'.format(indx,latOff,lonOff)) + yoffset = -(ds_WRFRef[name_lat][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]]-corners_lat[indx])\ + *(dy_parent/( ds_WRFRef[name_lat][0,ds_WRFRef['dFE_jindxs'][indx]+1,ds_WRFRef['dFE_iindxs'][indx]] + -ds_WRFRef[name_lat][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]])) + xoffset = -(ds_WRFRef[name_lon][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]]-corners_lon[indx])\ + *(dx_parent/( ds_WRFRef[name_lon][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]+1] + -ds_WRFRef[name_lon][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]])) + if(mpi_rank == 0): + print('\t corner({:d}): yOff = {:f}, xOff = {:f}'.format(indx,yoffset.values,xoffset.values)) + if indx<2: ## 0=southwest, or 1=southeast corner + if yoffset < 0.0: #FE domain SW/SE corner is south of the closest wrf cell, decrement the bounding jindx + ds_WRFRef['dFE_jindxs'][indx]-=1 + if indx < 1: + if xoffset < 0.0: #FE domain SW corner is west of the closest wrf cell, decrement the bounding iindx + ds_WRFRef['dFE_iindxs'][indx]-=1 + else: + if xoffset > 0.0: #FE domain SE corner is east of the closest wrf cell, increment the bounding iindx + ds_WRFRef['dFE_iindxs'][indx]+=1 + else: ## 3=northwest, or 2=northeast corner + if yoffset > 0.0: #FE domain NE/NW corner is north of the closest wrf cell, increment the bounding jindx + ds_WRFRef['dFE_jindxs'][indx]+=1 + if indx < 3: + if xoffset > 0.0: #FE domain NE corner is east of the closest wrf cell, increment the bounding iindx + ds_WRFRef['dFE_iindxs'][indx]+=1 + else: + if xoffset < 0.0: #FE domain NW corner is west of the closest wrf cell, decrement the bounding iindx + ds_WRFRef['dFE_iindxs'][indx]-=1 if(mpi_rank == 0): - print('\t corner({:d}): yOff = {:f}, xOff = {:f}'.format(indx,yoffset.values,xoffset.values)) - if indx<2: ## 0=southwest, or 1=southeast corner - if yoffset < 0.0: #FE domain SW/SE corner is south of the closest wrf cell, decrement the bounding jindx - ds_WRFRef['dFE_jindxs'][indx]-=1 - if indx < 1: - if xoffset < 0.0: #FE domain SW corner is west of the closest wrf cell, decrement the bounding iindx - ds_WRFRef['dFE_iindxs'][indx]-=1 - else: - if xoffset > 0.0: #FE domain SE corner is east of the closest wrf cell, increment the bounding iindx - ds_WRFRef['dFE_iindxs'][indx]+=1 - else: ## 3=northwest, or 2=northeast corner - if yoffset > 0.0: #FE domain NE/NW corner is north of the closest wrf cell, increment the bounding jindx - ds_WRFRef['dFE_jindxs'][indx]+=1 - if indx < 3: - if xoffset > 0.0: #FE domain NE corner is east of the closest wrf cell, increment the bounding iindx - ds_WRFRef['dFE_iindxs'][indx]+=1 - else: - if xoffset < 0.0: #FE domain NW corner is west of the closest wrf cell, decrement the bounding iindx - ds_WRFRef['dFE_iindxs'][indx]-=1 -if(mpi_rank == 0): - print('Bounding-box corrected offsets:') -for indx in range(len(corners_lat)): - latOff = ds_WRFRef['XLAT'][0,ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values]-corners_lat[indx] - lonOff = ds_WRFRef['XLONG'][0,ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values]-corners_lon[indx] - yoffset = -(ds_WRFRef['XLAT'][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]]-corners_lat[indx])\ - *(ds_WRFRef.attrs['DY']/( ds_WRFRef['XLAT'][0,ds_WRFRef['dFE_jindxs'][indx]+1,ds_WRFRef['dFE_iindxs'][indx]] - -ds_WRFRef['XLAT'][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]])) - xoffset = -(ds_WRFRef['XLONG'][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]]-corners_lon[indx])\ - *(ds_WRFRef.attrs['DX']/( ds_WRFRef['XLONG'][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]+1] - -ds_WRFRef['XLONG'][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]])) - if(mpi_rank == 0): - print('corner({:d}):latOff = {:f}, lonOff = {:f} -- yOff = {:f}, xOff = {:f}'.format(indx,latOff,lonOff,yoffset.values,xoffset.values)) - if indx == 0: - ll_yoffset=yoffset.values - ll_xoffset=xoffset.values -#print('WRF:') -#print('(j,i):') -for indx in range(len(corners_lat)): + print('Bounding-box corrected offsets:') + + for indx in range(len_corners): + latOff = ds_WRFRef[name_lat][0,ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values]-corners_lat[indx] + lonOff = ds_WRFRef[name_lon][0,ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values]-corners_lon[indx] + yoffset = -(ds_WRFRef[name_lat][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]]-corners_lat[indx])\ + *(dy_parent/( ds_WRFRef[name_lat][0,ds_WRFRef['dFE_jindxs'][indx]+1,ds_WRFRef['dFE_iindxs'][indx]] + -ds_WRFRef[name_lat][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]])) + xoffset = -(ds_WRFRef[name_lon][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]]-corners_lon[indx])\ + *(dx_parent/( ds_WRFRef[name_lon][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]+1] + -ds_WRFRef[name_lon][0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]])) + if(mpi_rank == 0): + print('corner({:d}):latOff = {:f}, lonOff = {:f} -- yOff = {:f}, xOff = {:f}'.format(indx,latOff,lonOff,yoffset.values,xoffset.values)) + if indx == 0: + ll_yoffset=yoffset.values + ll_xoffset=xoffset.values + for indx in range(len_corners): + if(mpi_rank == 0): + print('{:d},{:d}'.format(ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values)) + +else: + + for indx in range(len_corners): + if(mpi_rank == 0): + print(f"corner({indx}) @ WRF({ds_WRFRef['dFE_jindxs'][indx].values},{ds_WRFRef['dFE_iindxs'][indx].values})") + print('[WRF,corner({:d})]: yPos = [{:f},{:f}], xPos = [{:f},{:f}]'.format(indx,ds_WRFRef[name_lat][0,0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]].values, + corners_y[indx], + ds_WRFRef[name_lon][0,0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]].values, + corners_y[indx])) + + ### Compute FE-domain corner lat/lon offsets from closest dsWRFRef cell-centered lat/lons + print('\t Pre-correction offsets:') + for indx in range(len_corners): + xcoordOff = ds_WRFRef[name_lat][0,0,ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values]-corners_y[indx] + ycoordOff = ds_WRFRef[name_lon][0,0,ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values]-corners_x[indx] + if(mpi_rank == 0): + print('\t corner({:d}): ycoordOff = {:f}, xcoordOff = {:f}'.format(indx,ycoordOff,xcoordOff)) + yoffset = -(ds_WRFRef[name_lat][0,0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]]-corners_y[indx])\ + *(dy_parent/( ds_WRFRef[name_lat][0,0,ds_WRFRef['dFE_jindxs'][indx]+1,ds_WRFRef['dFE_iindxs'][indx]] + -ds_WRFRef[name_lat][0,0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]])) + xoffset = -(ds_WRFRef[name_lon][0,0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]]-corners_x[indx])\ + *(dx_parent/( ds_WRFRef[name_lon][0,0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]+1] + -ds_WRFRef[name_lon][0,0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]])) + if(mpi_rank == 0): + print('\t corner({:d}): yOff = {:f}, xOff = {:f}'.format(indx,yoffset.values,xoffset.values)) + if indx<2: ## 0=southwest, or 1=southeast corner + if yoffset < 0.0: #FE domain SW/SE corner is south of the closest wrf cell, decrement the bounding jindx + ds_WRFRef['dFE_jindxs'][indx]-=1 + if indx < 1: + if xoffset < 0.0: #FE domain SW corner is west of the closest wrf cell, decrement the bounding iindx + ds_WRFRef['dFE_iindxs'][indx]-=1 + else: + if xoffset > 0.0: #FE domain SE corner is east of the closest wrf cell, increment the bounding iindx + ds_WRFRef['dFE_iindxs'][indx]+=1 + else: ## 3=northwest, or 2=northeast corner + if yoffset > 0.0: #FE domain NE/NW corner is north of the closest wrf cell, increment the bounding jindx + ds_WRFRef['dFE_jindxs'][indx]+=1 + if indx < 3: + if xoffset > 0.0: #FE domain NE corner is east of the closest wrf cell, increment the bounding iindx + ds_WRFRef['dFE_iindxs'][indx]+=1 + else: + if xoffset < 0.0: #FE domain NW corner is west of the closest wrf cell, decrement the bounding iindx + ds_WRFRef['dFE_iindxs'][indx]-=1 if(mpi_rank == 0): - print('{:d},{:d}'.format(ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values)) + print('Bounding-box corrected offsets:') + for indx in range(len_corners): + xcoordOff = ds_WRFRef[name_lat][0,0,ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values]-corners_y[indx] + ycoordOff = ds_WRFRef[name_lon][0,0,ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values]-corners_x[indx] + yoffset = -(ds_WRFRef[name_lat][0,0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]]-corners_y[indx])\ + *(dy_parent/( ds_WRFRef[name_lat][0,0,ds_WRFRef['dFE_jindxs'][indx]+1,ds_WRFRef['dFE_iindxs'][indx]] + -ds_WRFRef[name_lat][0,0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]])) + xoffset = -(ds_WRFRef[name_lon][0,0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]]-corners_x[indx])\ + *(dx_parent/( ds_WRFRef[name_lon][0,0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]+1] + -ds_WRFRef[name_lon][0,0,ds_WRFRef['dFE_jindxs'][indx],ds_WRFRef['dFE_iindxs'][indx]])) + for indx in range(len_corners): + if(mpi_rank == 0): + print('corner({:d}):ycoordOff = {:f}, xcoordOff = {:f} -- yOff = {:f}, xOff = {:f}'.format(indx,ycoordOff,xcoordOff,yoffset.values,xoffset.values)) + if indx == 0: + ll_yoffset=yoffset.values + ll_xoffset=xoffset.values + for indx in range(len_corners): + if(mpi_rank == 0): + print('{:d},{:d}'.format(ds_WRFRef['dFE_jindxs'][indx].values,ds_WRFRef['dFE_iindxs'][indx].values)) #Nesting configuration parameters ll_jindx=ds_WRFRef['dFE_jindxs'].min(dim='corners').values @@ -242,33 +395,33 @@ i_extent=ds_WRFRef['dFE_iindxs'].max(dim='corners').values-ll_iindx ##Ensure WRF interpolation area extents will entirely encompass FE target x & y domain -y_distWRF = (j_extent-1)*ds_WRFRef.attrs['DY']-ll_yoffset -x_distWRF = (i_extent-1)*ds_WRFRef.attrs['DX']-ll_xoffset +y_distWRF = (j_extent-1)*dy_parent-ll_yoffset +x_distWRF = (i_extent-1)*dx_parent-ll_xoffset dxFE=(ds_FEGrid['xPos'][0,0,1]-ds_FEGrid['xPos'][0,0,0]).values dyFE=(ds_FEGrid['yPos'][0,1,0]-ds_FEGrid['yPos'][0,0,0]).values x_distFE = (ds_FEGrid.sizes['xIndex']-1)*dxFE y_distFE = (ds_FEGrid.sizes['yIndex']-1)*dyFE while x_distWRF <= x_distFE: i_extent += 1 - x_distWRF = (i_extent-1)*ds_WRFRef.attrs['DX']-ll_xoffset + x_distWRF = (i_extent-1)*dx_parent-ll_xoffset while y_distWRF <= y_distFE: j_extent += 1 - y_distWRF = (j_extent-1)*ds_WRFRef.attrs['DY']-ll_yoffset + y_distWRF = (j_extent-1)*dy_parent-ll_yoffset if(mpi_rank == 0): if (ll_jindx < 0): - print(f"Southern FE domain boundary coordinate falls outside of the provided WRF domain, exiting.") + print(f"Southern FE nested domain boundary coordinate falls outside of the provided parent domain, exiting.") exit() elif (ll_iindx < 0): - print(f"Western Ft domain boundary coordinate falls outside of the provided WRF domain, exiting.") + print(f"Western FE nested domain boundary coordinate falls outside of the provided parent domain, exiting.") exit() - elif (ll_jindx+j_extent > ds_WRFRef.sizes['south_north']): - print(f"Northern Ft domain boundary coordinate falls outside of the provided WRF domain, exiting.") + elif (ll_jindx+j_extent > Ngpy): + print(f"Northern FE nested domain boundary coordinate falls outside of the provided parent domain, exiting.") exit() - elif (ll_iindx+i_extent > ds_WRFRef.sizes['west_east']): - print(f"Eastern Ft domain boundary coordinate falls outside of the provided WRF domain, exiting.") + elif (ll_iindx+i_extent > Ngpx): + print(f"Eastern FE nested domain boundary coordinate falls outside of the provided parent domain, exiting.") exit() - else: #All set to perform strictly interpolation in the horizontal of WRF outputs to FE domain + else: #All set to perform strictly interpolation in the horizontal of parent outputs to nested FE domain print('ll: ({:d},{:d})'.format(ll_jindx,ll_iindx)) print('extents: ({:d},{:d})'.format(j_extent,i_extent)) print('y,x offsets: ({:f},{:f})'.format(ll_yoffset,ll_xoffset)) @@ -276,8 +429,8 @@ ###################################################################################################################### ### Define the Cartesian southwest corner origin (x,y) WRF coordinate system for the horizontal FE-bounding domain ### ###################################################################################################################### -xWRF,stepX=np.linspace((ll_iindx+0.5)*ds_WRFRef.attrs['DX'],(ll_iindx+0.5+i_extent)*ds_WRFRef.attrs['DX'],i_extent,endpoint=False,retstep=True) -yWRF,stepY=np.linspace((ll_jindx+0.5)*ds_WRFRef.attrs['DY'],(ll_jindx+0.5+j_extent)*ds_WRFRef.attrs['DY'],j_extent,endpoint=False,retstep=True) +xWRF,stepX=np.linspace((ll_iindx+0.5)*dx_parent,(ll_iindx+0.5+i_extent)*dx_parent,i_extent,endpoint=False,retstep=True) +yWRF,stepY=np.linspace((ll_jindx+0.5)*dy_parent,(ll_jindx+0.5+j_extent)*dy_parent,j_extent,endpoint=False,retstep=True) print(xWRF,'\n',yWRF,'\n') print(stepX,stepY) @@ -285,9 +438,6 @@ YvWRF,XvWRF=np.meshgrid(yWRF,xWRF, sparse=False, indexing='ij') print(XvWRF.shape,XvWRF.shape) - -print(ds_WRFRef['HGT'][0,ll_jindx:ll_jindx+j_extent,ll_iindx:ll_iindx+i_extent].values) - #################################################################################### ### Map the target FE domain into the WRF bounding-grid relative x,y coordinates ### #################################################################################### @@ -330,15 +480,34 @@ NzWRFInterp = 200 #275 zRect = np.linspace(zBottom,zTop,NzWRFInterp) if(mpi_rank == 0) and DEBUG_COUPLER: + print(f"zBottom,zTop,NzWRFInterp={zBottom},{zTop},{NzWRFInterp}") print(zRect[0],zRect[-1]) +## Establish a kMaxPrnt that minimizes the length of the interpolating function (for performance) +if parent_model == 1: + jPs = ds_WRFRef['dFE_jindxs'].values[0] + jPe = ds_WRFRef['dFE_jindxs'].values[-1] + iPs = ds_WRFRef['dFE_iindxs'].values[0] + iPe = ds_WRFRef['dFE_iindxs'].values[1] + if np.min(ds_WRFRef['zPos'][0,-1,jPs:jPe,iPs:iPe].values, axis=(0,1)) > zTop and zTop > np.max(ds_WRFRef['zPos'][0,0,jPs:jPe,iPs:iPe].values,axis=(0,1)): + kMaxPrnt = np.min(np.where(np.min(ds_WRFRef['zPos'][0,:,jPs:jPe,iPs:iPe].values,axis=(1,2))>zTop))+1 + print(f"Established kMaxPrnt = {kMaxPrnt} of {ds_WRFRef.sizes['zIndex']} total k-levels...") + else: + print(f"Error zTop = {zTop} is not within parent domain vertical bounds.\n Exiting Now!") + exit() ############################################################################## ### Create lists of relevant variable names in the WRF-FE coupling process ### ############################################################################## -varsList = ['Z','ALT','U','V','W','T','QVAPOR','QCLOUD'] -surfVarsList = ['TSK','Q2','HGT','PSFC'] #Note: Q2 in absence of QVG (which is not in wrfout by deafult) from WRF -FEvarsList = ['rho','u','v','w','theta','qv','ql'] -FEsurfVarsList = ['tskin','qskin','topoWRF','psfc','SeaMask'] +fe_low_tke = 1.0e-10; +FEvarsList = ['rho','u','v','w','theta','qv','ql','TKE_0'] +FEsurfVarsList = ['tskin','qskin'] +FEvar_mp = ['qv','ql'] +if (parent_model == 0): + varsList = ['Z','ALT','U','V','W','T','QVAPOR','QCLOUD','QKE'] + surfVarsList = ['TSK','Q2','HGT','PSFC'] #Note: Q2 in absence of QVG (which is not in wrfout by deafult) from WRF +elif (parent_model == 1): + varsList = ['zPos','rho','u','v','w','theta','qv','ql','TKE_0'] + surfVarsList = ['tskin','qskin','topoPos'] ####################################################################### ### Finally go ahead and create the initial and boundary conditions ### @@ -352,36 +521,57 @@ if not(os.path.isfile(bdyFileName)): print('{:d}{:d}: {:s} does not exist, creating it...'.format(mpi_rank, mpi_size, bdyFileName)) print("{:d}{:d}: Working on file {:s}".format(mpi_rank, mpi_size, files_list[Bdy_file_num])) - ds_ref = xr.open_mfdataset(files_list[Bdy_file_num],combine='nested',concat_dim='Time') + #ds_ref = xr.open_mfdataset(files_list[Bdy_file_num],combine='nested',concat_dim=name_concat_dim) + ds_ref = xr.open_dataset(files_list[Bdy_file_num]) t0s = time.perf_counter() - dsWRF=interpWRFToGrids(ds_ref,it0,varsList,surfVarsList,zRect,ll_iindx,i_extent,ll_jindx,j_extent) + if parent_model == 0: + dsWRF=interpWRFToGrids(ds_ref,it0,varsList,surfVarsList,zRect,ll_iindx,i_extent,ll_jindx,j_extent) + else: + dsWRF=interpFEToGrids(ds_ref,it0,varsList,surfVarsList,zRect,ll_iindx,i_extent,ll_jindx,j_extent,kMaxPrnt) + t0e = time.perf_counter() print('{:d}/{:d}: t0_elapsed = {:f} (s)'.format(mpi_rank, mpi_size, t0e-t0s)) t1s = time.perf_counter() - ds=copyAndTranspose(dsWRF) + if parent_model == 0: ### Only needed if parent model is WRF + ds=copyAndTranspose(dsWRF) + else: + ds=dsWRF t1e = time.perf_counter() print('{:d}/{:d}: t1_elapsed = {:f} (s)'.format(mpi_rank, mpi_size, t1e-t1s)) t2s = time.perf_counter() - dsFENew=interp2DForFE(ds,ds_FEGrid,XvWRF,YvWRF,xVec,yVec) + #dsFENew=interp2DForFE(ds,ds_FEGrid,XvWRF,YvWRF,xVec,yVec) + dsFENew=interp2DForFE(ds,ds_FEGrid,XvWRF,YvWRF,xVec,yVec,parent_model) t2e = time.perf_counter() print('{:d}/{:d}: t2_elapsed = {:f} (s)'.format(mpi_rank, mpi_size, t2e-t2s)) t3s = time.perf_counter() - dsFEFinal=create_dsFEFinal(ds_FEGrid) - verticalInterpFinal(ds_FEGrid,dsFENew,dsFEFinal,zRect) + dsFEFinal=create_dsFEFinal(ds_FEGrid,parent_model) + verticalInterpFinal(ds_FEGrid,dsFENew,dsFEFinal,zRect,parent_model) if 'BuildingMask' in list(dsFEFinal.variables): for var in ['u','v','w','ql','TKE_0']: if var in list(dsFEFinal.variables): dsFEFinal[var][:,:,:]=dsFEFinal[var][:,:,:]*np.where((dsFEFinal['BuildingMask'][:,:,:]>1e-3),0.0,1.0) + if (not nest_tke_opt): + print(f"Zeroing out TKE_0 since nest_tke_opt={nest_tke_opt}") + dsFEFinal['TKE_0'][:,:,:] = fe_low_tke; + else: # clip TKE_0 to avoid very small and negative values + dsFEFinal['TKE_0'][:,:,:] = np.clip(dsFEFinal['TKE_0'][:,:,:],a_min=fe_low_tke,a_max=None) + # ensure moisture and hydrometeors are not negative + for var_mp in FEvar_mp: + if var_mp in FEvarsList: + dsFEFinal[var_mp][:,:,:] = np.clip(dsFEFinal[var_mp][:,:,:],a_min=0.0,a_max=None) t3e = time.perf_counter() print('{:d}/{:d}: t3_elapsed = {:f} (s)'.format(mpi_rank, mpi_size, t3e-t3s)) addTimeDim_FEfinal(dsFEFinal) if Bdy_file_num == 0: - timeLabel="{:02d}{:02d}{:02d}UTC".format(timeHour0,timeMinute0,timeSecond0) + if parent_model == 0: + timeLabel="{:02d}{:02d}{:02d}UTC".format(timeHour0,timeMinute0,timeSecond0) + elif parent_model == 1: + timeLabel=f"{itMin}" dsFEFinal.to_netcdf(ICBC_dir+'FE_interp_{:s}.{:d}'.format(timeLabel,0),format='NETCDF4', encoding={'xIndex': {'dtype': 'i4'},'yIndex': {'dtype': 'i4'},'zIndex': {'dtype': 'i4'}}) t4s = time.perf_counter() - ds_Bdy=create_dsBdy(ds_FEGrid,FEvarsList,FEsurfVarsList) + ds_Bdy=create_dsBdy(ds_FEGrid,FEvarsList,FEsurfVarsList,parent_model) t4e = time.perf_counter() print('{:d}/{:d}: t4_elapsed = {:f} (s)'.format(mpi_rank, mpi_size, t4e-t4s)) t5s = time.perf_counter() diff --git a/scripts/python_utilities/coupler/GeoSpec.py b/scripts/python_utilities/coupler/GeoSpec.py index 60117634..5b1525af 100644 --- a/scripts/python_utilities/coupler/GeoSpec.py +++ b/scripts/python_utilities/coupler/GeoSpec.py @@ -21,7 +21,7 @@ name_dom = params["name_dom"] gis_root = params["gis_root"] gis_file = params["gis_file"] -nlcd_name = params["nlcd_name"] +landcover_table = params["landcover_table"] water_cats = params["water_cats"] urban_opt = params["urban_opt"] @@ -34,7 +34,7 @@ # derived paths -file_nlcd = gis_root + nlcd_name +file_nlcd = gis_root + landcover_table FE_new_nc = FE_dataset_path + name_dom + name_dom_add + '.nc' FE_plot = FE_dataset_path + name_dom + name_dom_add + '_geospec.png' print('FE_new_nc:', FE_new_nc) diff --git a/scripts/python_utilities/coupler/LandCoverMetadata_NLCD16.csv b/scripts/python_utilities/coupler/LandCoverMetadata_NLCD16.csv index cb2f2e36..fc79f29a 100755 --- a/scripts/python_utilities/coupler/LandCoverMetadata_NLCD16.csv +++ b/scripts/python_utilities/coupler/LandCoverMetadata_NLCD16.csv @@ -1,18 +1,18 @@ -Value,Type,z0 -0,Unclassified,0.03 -11,Open Water,0.001 -12,Perennial Snow/Ice,0.012 -21,"Developed, Open Space",0.05 -22,"Developed, Low Intensity",0.33 -23,"Developed, Medium Intensity",0.5 -24,Developed High Intensity,0.39 -31,Barren Land,0.09 -41,Deciduous Forest,0.65 -42,Evergreen Forest,0.72 -43,Mixed Forest,0.71 -52,Shrub/Scrub,0.12 -71,Herbaceous,0.04 -81,Hay/Pasture,0.06 -82,Cultivated Crops,0.05 -90,Woody Wetlands,0.55 -95,Emergent Herbaceous Wetlands,0.11 +Value,Type,z0,z0urbanLES +0,Unclassified,0.03,0.0 +11,Open Water,0.001,0.0 +12,Perennial Snow/Ice,0.012,0.0 +21,"Developed, Open Space",0.05,0.05 +22,"Developed, Low Intensity",0.33,0.05 +23,"Developed, Medium Intensity",0.5,0.05 +24,Developed High Intensity,0.39,0.05 +31,Barren Land,0.09,0.0 +41,Deciduous Forest,0.65,0.0 +42,Evergreen Forest,0.72,0.0 +43,Mixed Forest,0.71,0.0 +52,Shrub/Scrub,0.12,0.0 +71,Herbaceous,0.04,0.0 +81,Hay/Pasture,0.06,0.0 +82,Cultivated Crops,0.05,0.0 +90,Woody Wetlands,0.55,0.0 +95,Emergent Herbaceous Wetlands,0.11,0.0 diff --git a/scripts/python_utilities/coupler/SimGrid.py b/scripts/python_utilities/coupler/SimGrid.py index 3061ec7d..fef5d447 100644 --- a/scripts/python_utilities/coupler/SimGrid.py +++ b/scripts/python_utilities/coupler/SimGrid.py @@ -11,6 +11,7 @@ import math import struct from scipy.interpolate import RectBivariateSpline, NearestNDInterpolator +from skimage.measure import block_reduce from couplingUtils import * ####################################### @@ -29,6 +30,9 @@ urban_opt = params["urban_opt"] FE_new_nc_path = params["FE_new_nc_path"] name_dom_add = params["name_dom_add"] +urban_heatRedis_opt = params["urban_heatRedis_opt"] +landcover_table = params["landcover_table"] +topo_average_opt = params["topo_average_opt"] save_plot_opt = params["save_plot_opt"] ####################################### @@ -85,16 +89,16 @@ ## -npx_inc = int(d_xi/dx_inter) -npy_inc = int(d_eta/dy_inter) -if (npx_inc==0): - npx_inc = d_xi/dx_inter - npy_inc = d_eta/dy_inter - x_e = x_s + int(np.ceil(Nx*npx_inc)) - y_e = y_s + int(np.ceil(Ny*npy_inc)) -else: +if ((d_xi % dx_inter) == 0.0 and d_xi >= dx_inter): + npx_inc = int(d_xi/dx_inter) + npy_inc = int(d_eta/dy_inter) x_e = x_s + Nx*npx_inc y_e = y_s + Ny*npy_inc +else: + npx_inc = d_xi/dx_inter + npy_inc = d_eta/dy_inter + x_e = x_s + int(np.ceil((Nx-1)*npx_inc)) + 1 + y_e = y_s + int(np.ceil((Ny-1)*npy_inc)) + 1 print('x_s,x_e,y_s,y_e=',x_s,x_e,y_s,y_e) box_indx = [x_s,x_e,x_e,x_s,x_s] @@ -113,14 +117,16 @@ print('Grid is not an even factor of GIS resolution, use interpolation') interp_flag = 1 print('interp_flag=',interp_flag) - +if (topo_average_opt==1): + print('Using block averaging for topography') + verticalDeformSwitch = int(str(FE_params['verticalDeformSwitch'][0])) print('verticalDeformSwitch=',verticalDeformSwitch) if (verticalDeformSwitch==1): c1 = float(str(FE_params['verticalDeformFactor'][0])) fCoeff = float(str(FE_params['verticalDeformQuadCoeff'][0])) else: - c1 = 0.0 + c1 = 1.0 fCoeff = 0.0 print('c1,fCoeff=',c1,fCoeff) @@ -130,24 +136,28 @@ ## Read in terrain elevation array topo = ds_GIS.topoPos.values if (interp_flag==0): - data_topo0 = topo[y_s:y_e:npy_inc,x_s:x_e:npx_inc] + if (topo_average_opt==0): + data_topo0 = topo[y_s:y_e:npy_inc,x_s:x_e:npx_inc] + else: + data_topo0 = block_reduce(topo[y_s:y_e,x_s:x_e], block_size=(int(d_eta/dy_inter),int(d_xi/dx_inter)), func=np.mean) else: xPos_2d_dom_ori = xPos_2d[y_s:y_e,x_s:x_e] yPos_2d_dom_ori = yPos_2d[y_s:y_e,x_s:x_e] - topo_dom_ori = topo[y_s:y_e,x_s:x_e] print('xPos_2d_dom_ori.shape=',xPos_2d_dom_ori.shape) - f_topo = RectBivariateSpline(xPos_2d_dom_ori[0,:], yPos_2d_dom_ori[:,0], topo_dom_ori.T, kx=3, ky=3) - xPos_1d_new = np.arange(x_box_corners[0],x_box_corners[1],d_xi) yPos_1d_new = np.arange(y_box_corners[0],y_box_corners[2],d_eta) - - data_topo0_b = f_topo(xPos_1d_new, yPos_1d_new).T xPos_2d_new_b, yPos_2d_new_b = np.meshgrid(xPos_1d_new, yPos_1d_new) - - data_topo0 = data_topo0_b[0:Ny,0:Nx] xPos_2d_new = xPos_2d_new_b[0:Ny,0:Nx] yPos_2d_new = yPos_2d_new_b[0:Ny,0:Nx] + topo_dom_ori = topo[y_s:y_e,x_s:x_e] + if (topo_average_opt==0): + f_topo = RectBivariateSpline(xPos_2d_dom_ori[0,:], yPos_2d_dom_ori[:,0], topo_dom_ori.T, kx=3, ky=3) + data_topo0_b = f_topo(xPos_1d_new, yPos_1d_new).T + data_topo0 = data_topo0_b[0:Ny,0:Nx] + else: + data_topo0 = block_average_topo(topo_dom_ori,dx_inter,dy_inter,d_xi,d_eta,Nx,Ny) + data_topo = smoothTerrain(data_topo0,d_xi) topoPos_min = np.min(data_topo,axis=(1,0)) @@ -176,7 +186,10 @@ for kk in range(0,Nz): zPos_uni = kk*d_zeta + 0.5*d_zeta - zPos_str[kk] = zDeform(zPos_uni,zbot,ztop,c1,fCoeff) + if (verticalDeformSwitch==1): + zPos_str[kk] = zDeform(zPos_uni,zbot,ztop,c1,fCoeff) + else: + zPos_str[kk] = zPos_uni if (kk==0): print('kk,zPos_str,dz=',kk,',',zPos_str[kk],', -') else: @@ -189,7 +202,10 @@ for i in range(Nx): zbot = data_topo[j,i] zPos_uni = np.linspace(0.5*d_zeta,(Nz-0.5)*d_zeta,Nz) - zPos_str = zDeform(zPos_uni,zbot,ztop,c1,fCoeff) + if (verticalDeformSwitch==1): + zPos_str = zDeform(zPos_uni,zbot,ztop,c1,fCoeff) + else: + zPos_str = zPos_uni zarr[:,j,i] = zPos_str if (j==0) and (i==0): for k in range(Nz): @@ -236,7 +252,7 @@ if (interp_flag==0): data_bmask = bdg_heights[y_s:y_e:npy_inc,x_s:x_e:npx_inc] else: - f_bdg = NearestNDInterpolator(list(zip(xPos_2d_dom_ori.flatten(), yPos_2d_dom_ori.flatten())), data_bmask[y_s:y_e,x_s:x_e].flatten()) + f_bdg = NearestNDInterpolator(list(zip(xPos_2d_dom_ori.flatten(), yPos_2d_dom_ori.flatten())), bdg_heights[y_s:y_e,x_s:x_e].flatten()) data_bmask = f_bdg(xPos_2d_new, yPos_2d_new) bdg3d_tmp = np.zeros((Nz,Ny,Nx),dtype=np.float32) @@ -300,6 +316,13 @@ lat_dom = lat_dom_b[0:Ny,0:Nx] lon_dom = lon_dom_b[0:Ny,0:Nx] +# Surface heat flux redistribution + +if (urban_opt == 1 and urban_heatRedis_opt == 1): + z0_original, z0_modified = read_lc_table(landcover_table) + z1 = zarr[0,:,:]-data_topo + shfr = SHFR_process_polygons(data_landc,data_bmask,z1,z0_original,z0_modified) + # Save to netCDF file ds_data = xr.Dataset() @@ -315,6 +338,8 @@ if (urban_opt == 1): ds_data['BuildingMask']= xr.DataArray(bdg3d_tmp.astype(dtype=np.float32),dims=(['zIndex','yIndex','xIndex'])) ds_data['BuildingHeights']= xr.DataArray(data_bmask.astype(dtype=np.float32),dims=(['yIndex','xIndex'])) + if (urban_heatRedis_opt == 1): + ds_data['UrbanHeatRedis']= xr.DataArray(shfr.astype(dtype=np.float32),dims=(['yIndex','xIndex'])) ds_data['lat']= xr.DataArray(lat_dom.astype(dtype=np.float64),dims=(['yIndex','xIndex'])) ds_data['lon']= xr.DataArray(lon_dom.astype(dtype=np.float64),dims=(['yIndex','xIndex'])) ds_data['xIndex']= xr.DataArray(np.arange(0,xarr.shape[2],dtype=np.int32),dims='xIndex') diff --git a/scripts/python_utilities/coupler/couplingUtils.py b/scripts/python_utilities/coupler/couplingUtils.py index 5d6bd2ae..d2c0d886 100644 --- a/scripts/python_utilities/coupler/couplingUtils.py +++ b/scripts/python_utilities/coupler/couplingUtils.py @@ -1,11 +1,13 @@ import argparse import math import time -from scipy.ndimage import gaussian_filter +import scipy.ndimage as ndimage import numpy as np import xarray as xr +import pandas as pd from scipy import interpolate from scipy.interpolate import RectBivariateSpline +from scipy.interpolate import BSpline, make_interp_spline def parse_args(): """ parse the command line arguments """ @@ -43,6 +45,42 @@ def zDeform(zRect, zGround, zCeiling, c1, fCoeff): return zStretch +def block_average_topo(topo_dom_ori,dx_inter,dy_inter,d_xi,d_eta,Nx,Ny): + ny_src, nx_src = topo_dom_ori.shape + data_topo0 = np.zeros((Ny, Nx), dtype=float) + x_edges_src = np.arange(nx_src + 1) * dx_inter + y_edges_src = np.arange(ny_src + 1) * dy_inter + x_edges_dst = np.arange(Nx + 1) * d_xi + y_edges_dst = np.arange(Ny + 1) * d_eta + for j in range(Ny): + y0 = y_edges_dst[j] + y1 = y_edges_dst[j + 1] + iy0 = np.searchsorted(y_edges_src, y0, side="right") - 1 + iy1 = np.searchsorted(y_edges_src, y1, side="left") + for i in range(Nx): + x0 = x_edges_dst[i] + x1 = x_edges_dst[i + 1] + ix0 = np.searchsorted(x_edges_src, x0, side="right") - 1 + ix1 = np.searchsorted(x_edges_src, x1, side="left") + val = 0.0 + wsum = 0.0 + for iy in range(max(iy0, 0), min(iy1 + 1, ny_src)): + ys0 = y_edges_src[iy] + ys1 = y_edges_src[iy + 1] + overlap_y = min(y1, ys1) - max(y0, ys0) + for ix in range(max(ix0, 0), min(ix1 + 1, nx_src)): + xs0 = x_edges_src[ix] + xs1 = x_edges_src[ix + 1] + overlap_x = min(x1, xs1) - max(x0, xs0) + if overlap_x > 0.0 and overlap_y > 0.0: + w = overlap_x * overlap_y + val += topo_dom_ori[iy, ix] * w + wsum += w + if wsum == 0.0: + raise ValueError(f'No overlap found for cell j={j}, i={i}.') + data_topo0[j, i] = val / wsum + return data_topo0 + def smoothTerrain(tPos0,dx): start = time.time() slopeThresh=math.tan(math.radians(35.0)) @@ -71,7 +109,7 @@ def smoothTerrain(tPos0,dx): min_col = max(0, max_idx[1] - n//2) max_col = min(tPos1.shape[1], max_idx[1] + n//2+1) local_area = tPos1[min_row:max_row, min_col:max_col] - local_blur_center = gaussian_filter(local_area, sigma=sigTry) + local_blur_center = ndimage.gaussian_filter(local_area, sigma=sigTry) for i in range(min_row, max_row): for j in range(min_col, max_col): if ((i, j) == max_idx): @@ -95,6 +133,31 @@ def smoothTerrain(tPos0,dx): print(f'Elapsed time [s]: {np.round(end-start,3)}') return tPos1 +def interpFEToGrids(ds_ref,it0,varsList,surfVarsList,zRect,ll_iindx,i_extent,ll_jindx,j_extent,kMaxPrnt): + verbose = False + surfVarDict = {'tskin':'tskin','qskin':'qskin','topoPos':'topoParent'} + dsFE = xr.Dataset() + zPrntArray = ds_ref['zPos'][it0,:kMaxPrnt,ll_jindx:ll_jindx+j_extent,ll_iindx:ll_iindx+i_extent].values + for var in varsList: + tmpArray=np.zeros((zRect.shape[0],j_extent,i_extent)) + if var != 'zPos': + varPrntArray=ds_ref[var][it0,:kMaxPrnt,ll_jindx:ll_jindx+j_extent,ll_iindx:ll_iindx+i_extent].values + for i in range(i_extent): + for j in range(j_extent): + f1=make_interp_spline(zPrntArray[:,j,i],varPrntArray[:,j,i], + k=3, bc_type='natural') + tmpArray[:,j,i] = f1(zRect) + dsFE[var] = xr.DataArray(tmpArray,dims=('zIndex','yIndex','xIndex')) + if verbose: + print(f"interpFEToGrids: {var}--({dsFE[var].values.shape})") + for surfVar in surfVarsList: + tmpArray=ds_ref[surfVar][it0,ll_jindx:ll_jindx+j_extent,ll_iindx:ll_iindx+i_extent].values + dsFE[surfVarDict[surfVar]] = xr.DataArray(tmpArray,dims=('yIndex','xIndex')) + if verbose: + print(f"interpFEToGrids: {surfVarDict[surfVar]}--({dsFE[surfVarDict[surfVar]].values.shape})") + + return dsFE + def interpWRFToGrids(ds_ref,it0,varsList,surfVarsList,zRect,ll_iindx,i_extent,ll_jindx,j_extent): dsWRF = xr.Dataset() for i in range(ll_iindx,ll_iindx+i_extent): @@ -112,42 +175,55 @@ def interpWRFToGrids(ds_ref,it0,varsList,surfVarsList,zRect,ll_iindx,i_extent,ll return dsWRF def get_dsWRFStandardZprof(it0,j0,i0,ds_ref,varsList,surfVarsList,zProf): - ds_ij=getFEProfileDS(getWRFProfileDS(it0,j0,i0,ds_ref,varsList,surfVarsList),varsList,surfVarsList,zProf) + ds_ij=getFEProfileDS(getWRFProfileDS(it0,j0,i0,ds_ref,varsList,surfVarsList,zProf[-1]),varsList,surfVarsList,zProf) return ds_ij def getFEProfileDS(dsWrf,varsList,surfVarsList,zFE): ##### Map (Interp/Extrap-olate) a collected set of WRF vertical profiles from ##### the WRF vertical coordinate to a specified cartesian z-coord (zFE) ds_ret=xr.Dataset() - fromRestart = False - if fromRestart: - varDict = {'Z':'zPos','U_1':'u','V_1':'v','W_1':'w','T':'theta','QVAPOR':'qv','QCLOUD':'ql','ALT':'rho','ALB':'BS_0','PB':'BS_4'} - else: - varDict = {'Z':'zPos','U':'u','V':'v','W':'w','T':'theta','QVAPOR':'qv','QCLOUD':'ql','ALT':'rho'} + varDict = {'Z':'zPos','U':'u','V':'v','W':'w','T':'theta','QVAPOR':'qv','QCLOUD':'ql','ALT':'rho','QKE':'TKE_0'} surfVarDict = {'TSK':'tskin','Q2':'qskin','HGT':'topoWRF','T2':'t2','PSFC':'psfc'} #Note using Q2 instead of QVG since QVG not default in wrfout files + #print(f"kMaxPrnt = {dsWrf.sizes['bottom_top']}") for var in varsList: if var != 'Z': f1=interpolate.interp1d(dsWrf['Z'],dsWrf[var],kind='linear',fill_value='extrapolate') if var != 'ALT' and var != 'ALB': - ds_ret[varDict[var]] = xr.DataArray(f1(zFE),dims=['zIndex']) #,coords={'zIndex':np.array0:zFE.size} + ds_ret[varDict[var]] = xr.DataArray(f1(zFE),dims=['zIndex']) else: - ds_ret[varDict[var]] = xr.DataArray(1.0/f1(zFE),dims=['zIndex']) #,coords={'zIndex':np.array0:zFE.size} + ds_ret[varDict[var]] = xr.DataArray(1.0/f1(zFE),dims=['zIndex']) for surfVar in surfVarsList: ds_ret[surfVarDict[surfVar]] = xr.DataArray(dsWrf[surfVar]) + return ds_ret -def getWRFProfileDS(it,j,i,dsWrf,varsList,surfVarsList): ##### Destagger and collect a set of required WRF vertical profiles from a given i,j location in WRF +def getWRFProfileDS(it,j,i,dsWrf,varsList,surfVarsList,zTargTop): ##### Destagger and collect a set of required WRF vertical profiles from a given i,j location in WRF ds_ret=xr.Dataset() - fromRestart = False for var in varsList: if var == 'Z': - if fromRestart: - ds_ret[var] = xr.DataArray((0.5*(dsWrf.PH_1[it,0:-1,j,i]+dsWrf.PH_1[it,1:,j,i]) - +0.5*(dsWrf.PHB[it,0:-1,j,i]+dsWrf.PHB[it,1:,j,i]))/9.81, + ds_ret[var] = xr.DataArray((0.5*(dsWrf.PH[it,0:-1,j,i]+dsWrf.PH[it,1:,j,i]) + +0.5*(dsWrf.PHB[it,0:-1,j,i]+dsWrf.PHB[it,1:,j,i]))/9.81, + dims=(['bottom_top'])) + ## Determine a kMaxPrnt for this i,j to be used for the interpolation step in the calling function + kMaxPrnt = np.min(np.where(ds_ret[var].values>zTargTop))+1 + elif var == 'QKE': + tke_infile = dsWrf.get(var) + if tke_infile is not None: + tke_var = 1 + else: # try with TKE_PBL alternatively + tke_infile = dsWrf.get('TKE_PBL') + if tke_infile is not None: + tke_var = 2 + else: + tke_var = 0 + if (tke_var == 1): + ds_ret[var] = xr.DataArray(0.5*dsWrf[var][it,:,j,i], # QKE is 2.0*TKE dims=(['bottom_top'])) - else: - ds_ret[var] = xr.DataArray((0.5*(dsWrf.PH[it,0:-1,j,i]+dsWrf.PH[it,1:,j,i]) - +0.5*(dsWrf.PHB[it,0:-1,j,i]+dsWrf.PHB[it,1:,j,i]))/9.81, + elif (tke_var == 2): # TKE_PBL is bottom_top_stag + ds_ret[var] = xr.DataArray(0.5*(dsWrf['TKE_PBL'][it,0:-1,j,i]+dsWrf['TKE_PBL'][it,1:,j,i]), + dims=(['bottom_top'])) + else: # tke_var == 0 (no tke variable present) + ds_ret[var] = xr.DataArray(dsWrf['T'][it,:,j,i]*0.0+1e-10, dims=(['bottom_top'])) elif 'west_east_stag' in dsWrf[var].dims: ds_ret[var] = xr.DataArray(0.5*(dsWrf[var][it,:,j,i]+dsWrf[var][it,:,j,i+1]), @@ -166,7 +242,10 @@ def getWRFProfileDS(it,j,i,dsWrf,varsList,surfVarsList): ##### Destagger and col for surfVar in surfVarsList: ds_ret[surfVar] = xr.DataArray(dsWrf[surfVar][it,j,i])#, - return ds_ret + ## Determine a kMaxPrnt for this i,j to be used for the interpolation step in the calling function + + return ds_ret.isel(bottom_top=slice(0,kMaxPrnt)) + def copyAndTranspose(dsWRF): ds=xr.Dataset() @@ -177,7 +256,8 @@ def copyAndTranspose(dsWRF): ds[var]=xr.DataArray(dsWRF[var].values,dims=('yIndex','xIndex')) #Note no transpose needed here... return ds -def interp2DForFE(ds,ds_FE,XvWRF,YvWRF,xVec,yVec): +#def interp2DForFE(ds,ds_FE,XvWRF,YvWRF,xVec,yVec): +def interp2DForFE(ds,ds_FE,XvWRF,YvWRF,xVec,yVec,parent_model): k_val = 1 k_val_surf = 3 dsFENew=xr.Dataset() @@ -187,12 +267,19 @@ def interp2DForFE(ds,ds_FE,XvWRF,YvWRF,xVec,yVec): tmpVar3d = np.zeros((ds.sizes['zIndex'],ds_FE.sizes['yIndex'],ds_FE.sizes['xIndex'])) print(tmpVar3d.shape) for k in range(ds.sizes['zIndex']): - fInterp2 = RectBivariateSpline(YvWRF[:,0],XvWRF[0,:],ds[var][k,:,:].values.transpose(), kx=k_val, ky=k_val) + if parent_model == 0: + fInterp2 = RectBivariateSpline(YvWRF[:,0],XvWRF[0,:],ds[var][k,:,:].values.transpose(), kx=k_val, ky=k_val) + else: + fInterp2 = RectBivariateSpline(YvWRF[:,0],XvWRF[0,:],ds[var][k,:,:].values, kx=k_val, ky=k_val) tmp=fInterp2(yVec,xVec) tmpVar3d[k,:,:]=tmp dsFENew[var]=xr.DataArray(tmpVar3d,dims=('zIndex','yIndex','xIndex')) elif len(ds[var].sizes) == 2: - fInterp = RectBivariateSpline(YvWRF[:,0],XvWRF[0,:],ds[var].values.transpose(), kx=k_val_surf, ky=k_val_surf) + if parent_model == 0: + fInterp = RectBivariateSpline(YvWRF[:,0],XvWRF[0,:],ds[var].values.transpose(), kx=k_val_surf, ky=k_val_surf) + else: + fInterp = RectBivariateSpline(YvWRF[:,0],XvWRF[0,:],ds[var].values, kx=k_val_surf, ky=k_val_surf) + tmp=fInterp(yVec,xVec) print(tmp.shape) dsFENew[var]=xr.DataArray(tmp,dims=('yIndex','xIndex')) @@ -200,22 +287,21 @@ def interp2DForFE(ds,ds_FE,XvWRF,YvWRF,xVec,yVec): print('{:s} required {:f} s for ij-interpolation'.format(var,t1e-t1s)) return dsFENew -def create_dsFEFinal(ds_FE): +def create_dsFEFinal(ds_FE,parent_model): dsFEFinal=ds_FE.copy(deep=True) dsFEFinal.load() - for var in ['rho', 'u', 'v', 'w', 'theta', 'TKE_0', 'qv', 'pressure']: + for var in ['rho', 'u', 'v', 'w', 'theta', 'TKE_0', 'qv', 'ql', 'pressure']: dsFEFinal[var]=0.0*dsFEFinal['xPos'] for var in ['fricVel','htFlux','invOblen','qFlux']: - dsFEFinal[var]=0.0*dsFEFinal['z0m'] - for var in ['ql', 'XLAT','XLONG','topoWRF','t2','psfc']: - if var in list(dsFEFinal.variables): - if var in ['ql']: - dsFEFinal[var]=0.0*dsFEFinal['rho'] - elif var in ['XLAT','XLONG','topoWRF','t2','psfc']: - dsFEFinal[var]=0.0*dsFEFinal['tskin'] + dsFEFinal[var]=0.0*dsFEFinal['z0m'] + if (parent_model == 0): + for var in ['XLAT','XLONG','topoWRF','t2','psfc']: + if var in list(dsFEFinal.variables): + if var in ['XLAT','XLONG','topoWRF','t2','psfc']: + dsFEFinal[var]=0.0*dsFEFinal['tskin'] return dsFEFinal -def verticalInterpFinal(ds_FE,dsFENew,dsFEFinal,zRect): +def verticalInterpFinal(ds_FE,dsFENew,dsFEFinal,zRect,parent_model): z3d=ds_FE['zPos'][:,:,:].values.squeeze() for var in dsFENew.variables: print(var) @@ -244,13 +330,13 @@ def verticalInterpFinal(ds_FE,dsFENew,dsFEFinal,zRect): dsFEFinal[var]=xr.DataArray(tmp,dims=('yIndex','xIndex')) t1e = time.perf_counter() print('{:s} required {:f} s for vertical interpolation of the i,j-set'.format(var,t1e-t1s)) - if 'qv' in list(dsFEFinal.variables): - #Scale the water vapor mixing ratio from kg/kg (WRF) to g/kg (FE) - dsFEFinal['qv'] = 1e3*dsFEFinal['qv'] - if 'ql' in list(dsFEFinal.variables): - dsFEFinal['ql'] = 1e3*dsFEFinal['ql'] - if 'qskin' in list(dsFEFinal.variables): - dsFEFinal['qskin'] = 1e3*dsFEFinal['qskin'] + if (parent_model == 0): #Scale the water vapor mixing ratio from kg/kg (WRF) to g/kg (FE) + if 'qv' in list(dsFEFinal.variables): + dsFEFinal['qv'] = 1e3*dsFEFinal['qv'] + if 'ql' in list(dsFEFinal.variables): + dsFEFinal['ql'] = 1e3*dsFEFinal['ql'] + if 'qskin' in list(dsFEFinal.variables): + dsFEFinal['qskin'] = 1e3*dsFEFinal['qskin'] def interpolateIrregularVertical(zRect,fld3dRect,z3d): NzR,NyR,NxR = fld3dRect.shape @@ -262,7 +348,7 @@ def interpolateIrregularVertical(zRect,fld3dRect,z3d): tmpVar3d[:,j,i]=tmp return tmpVar3d -def create_dsBdy(ds_FE,FEvarsList,FEsurfVarsList): +def create_dsBdy(ds_FE,FEvarsList,FEsurfVarsList,parent_model): ds_Bdy=xr.Dataset() notit=0 for var in FEvarsList: @@ -274,12 +360,21 @@ def create_dsBdy(ds_FE,FEvarsList,FEsurfVarsList): ds_Bdy[var+'_XZH']=xr.DataArray(ds_FE[var][:,:,ds_FE.sizes['yIndex']-1,:]) ds_Bdy[var+'_XYL']=xr.DataArray(ds_FE[var][:,0,:,:]) ds_Bdy[var+'_XYH']=xr.DataArray(ds_FE[var][:,ds_FE.sizes['zIndex']-1,:,:]) - for surfVar in FEsurfVarsList: - print('{:s}: notit={:d}'.format(surfVar,notit)) - if surfVar in ['topoWRF','t2','psfc','tskin','qskin']: - notit +=1 - else: - ds_Bdy[surfVar]=xr.DataArray(ds_FE[surfVar][:,:]) + if (parent_model == 0): + for surfVar in FEsurfVarsList: + print('{:s}: notit={:d}'.format(surfVar,notit)) + if surfVar in ['topoWRF','topoParent','t2','psfc','tskin','qskin']: + notit +=1 + else: + ds_Bdy[surfVar]=xr.DataArray(ds_FE[surfVar][:,:]) + elif (parent_model == 1): + for surfVar in FEsurfVarsList: + print('{:s}: notit={:d}'.format(surfVar,notit)) + if surfVar in ['tskin','qskin']: + notit +=1 + else: + ds_Bdy[surfVar]=xr.DataArray(ds_FE[surfVar][:,:]) + return ds_Bdy def createBdysFrom3D(ds_Bdy,ds3D,FEvarsList,FEsurfVarsList): @@ -317,3 +412,65 @@ def addTimeDim_FEfinal(dsFEFinal): if len(dsFEFinal[var].values.shape) != 1: dsFEFinal[var] = dsFEFinal[var].expand_dims(dim={'time':1},axis=0) return + +def read_lc_table(filepath): + df = pd.read_csv(filepath) + z0_original = {int(k): float(v) for k, v in zip(df.iloc[:,0], df.iloc[:,2])} + z0_modified = {int(k): float(v) for k, v in zip(df.iloc[:,0], df.iloc[:,3])} + return z0_original, z0_modified + +def SHFR_process_polygons(landcover, buildings, z1, z0_original, z0_modified, nodata=0, N0 = 10, Nmin = 25, fmin = 0.10): + result = np.zeros_like(landcover, dtype=np.float32) + labels = np.zeros_like(landcover, dtype=np.int32) + current_label = 1 + #------- + r0 = 0 + r1 = 0 + r2 = 0 + r3 = 0 + r4 = 0 + #------- + for category in np.unique(landcover): + category_mask = landcover == category + category_labels, num_labels = ndimage.label(category_mask, structure = ndimage.generate_binary_structure(2,2)) + for i in range(1, num_labels + 1): + labels[category_labels == i] = current_label + current_label += 1 + num_polygons = current_label - 1 + print(f'Number of land cover polygons: {num_polygons}') + for label in range(1, num_polygons + 1): + polygon_mask = labels == label + total_area = np.sum(polygon_mask) + building_mask = (buildings > nodata) & polygon_mask + building_area = np.sum(building_mask) + no_building_area = total_area - building_area + Nreq = np.maximum(Nmin, total_area*fmin) + lc = landcover[polygon_mask][0] + if building_area == 0: + no_building_value = 1.0 + r0 += 1 + else: + if no_building_area <= N0: + no_building_value = 1.0 + r1 += 1 + else: + z = z1[polygon_mask].mean() + z0lc = z0_original[lc] + z0st = z0_modified[lc] if z0_modified[lc] > 0.0 else z0lc + factor_z0 = (np.log(z/z0st+1)*np.log(z/(0.1*z0st)+1)) / (np.log(z/z0lc+1)*np.log(z/(0.1*z0lc)+1)) + if no_building_area < Nreq: + w = (no_building_area-N0)/(Nreq-N0) + no_building_value = 1.0 + w * np.minimum( 4, (building_area/no_building_area)*factor_z0 ) + r2 += 1 + else: + if (building_area / no_building_area)*factor_z0 > 4: + no_building_value = 5.0 + r4 += 1 + else: + no_building_value = 1.0 + (building_area / no_building_area)*factor_z0 + r3 += 1 + if (z0_modified[lc] == 0.0 and no_building_value < 1.1): + no_building_value = 1.0 + result[polygon_mask & ~building_mask] = no_building_value + print(f'R0 = {r0}, R1 = {r1}, R2 = {r2}, R3 = {r3}, R4 = {r4}') + return result diff --git a/scripts/python_utilities/coupler/genicbcs.json b/scripts/python_utilities/coupler/genicbcs.json index 0ca9302b..04dbfe68 100644 --- a/scripts/python_utilities/coupler/genicbcs.json +++ b/scripts/python_utilities/coupler/genicbcs.json @@ -1,6 +1,11 @@ { "ICBC_dir": "/path_FEsim/ICBC/", "FE_simGrid": "/path_simgrid/simgrid.0", + "parent_model": 0, + "nest_tke_opt": true, + "_comment": "#######################################################", + "_comment": "### WRF as parent domain section (parent_model = 0) ###", + "_comment": "#######################################################", "WRF_PrntDir": "/path_wrf_data/", "WRF_PrntOutPrefix": "wrf_fasteddy_d01_", "date0": "2024-02-16", @@ -8,5 +13,15 @@ "timeMinute0": 0, "timeSecond0": 0, "secMax": 5401, - "secInc": 300 + "secInc": 300, + "_comment": "############################################################", + "_comment": "### FastEddy as parent domain section (parent_model = 1) ###", + "_comment": "############################################################", + "FE_PrntDir": "/path_fe_parent/", + "FE_PrntOutPrefix": "FE_parent_name", + "itMin": 108000, + "dt_FE": 0.025, + "outputFrequency": 10.0, + "timeLengthSec": 3600.0, + "ideal_opt": false } diff --git a/scripts/python_utilities/coupler/geospec.json b/scripts/python_utilities/coupler/geospec.json index 4969eb13..580feb0c 100644 --- a/scripts/python_utilities/coupler/geospec.json +++ b/scripts/python_utilities/coupler/geospec.json @@ -2,7 +2,7 @@ "name_dom": "DomainName", "gis_root": "/path_gis/", "gis_file": "inputs_gis.nc", - "nlcd_name": "LandCoverMetadata_NLCD16.csv", + "landcover_table": "LandCoverMetadata_NLCD16.csv", "water_cats": [11], "urban_opt": 0, "FE_dataset_path": "/path_geospec/", diff --git a/scripts/python_utilities/coupler/mpassit_to_fasteddy.json b/scripts/python_utilities/coupler/mpassit_to_fasteddy.json new file mode 100644 index 00000000..af3e4662 --- /dev/null +++ b/scripts/python_utilities/coupler/mpassit_to_fasteddy.json @@ -0,0 +1,12 @@ +{ + "MPASSIT_Dir": "./", + "MPASSIT_Prefix": "proc.", + "FE_PrntOutDir": "./", + "FE_PrntOutPrefix": "mpassit_", + "date0": "2024-02-16", + "timeHour0": 17, + "timeMinute0": 0, + "timeSecond0": 0, + "secMax": 5401, + "secInc": 300 +} diff --git a/scripts/python_utilities/coupler/mpassit_to_fasteddy.py b/scripts/python_utilities/coupler/mpassit_to_fasteddy.py new file mode 100644 index 00000000..b4ae3793 --- /dev/null +++ b/scripts/python_utilities/coupler/mpassit_to_fasteddy.py @@ -0,0 +1,65 @@ +import os, sys +import numpy as np +import xarray as xr +import argparse +import json +from netCDF4 import Dataset +import datetime as dt +import xarray as xr + +from couplingUtils import * + +# Read run parameters from json +args = parse_args() +with open(args.file) as file: + params = json.loads(file.read()) + +MPASSIT_Dir = params["MPASSIT_Dir"] +MPASSIT_Prefix = params["MPASSIT_Prefix"] +FE_PrntOutDir = params["FE_PrntOutDir"] +FE_PrntOutPrefix = params["FE_PrntOutPrefix"] +dateString = params["date0"] +timeHour0 = params["timeHour0"] +timeMinute0 = params["timeMinute0"] +timeSecond0 = params["timeSecond0"] +secMax = params["secMax"] +secInc = params["secInc"] + +# Define Constants +grav=9.80665 +p1000mb=100000.0 +rv=461.6 +rd=287.0 +cp=7.0*rd/2.0 +cv=cp-rd +rvovrd=rv/rd +cvpm=-1.0*(cv/cp) + +# Define a list of MPAS files to process from the specified coupler parameters +files_list_mpassit=[] +files_list_fe=[] + +year0 = int(dateString[0:4]) +month0 = int(dateString[5:7]) +day0 = int(dateString[8:10]) + +date_it = dt.datetime(year0,month0,day0,timeHour0,timeMinute0,timeSecond0) +for it in range(0,secMax,secInc): + dateString_it = str(date_it.year) + '-' + "{:02d}".format(date_it.month) + '-' + "{:02d}".format(date_it.day) + '_' + thistime_mpassit = "{:s}{:02d}.{:02d}.{:02d}".format(dateString_it,date_it.hour,date_it.minute,date_it.second) + thistime_fe = "{:s}{:02d}:{:02d}:{:02d}".format(dateString_it,date_it.hour,date_it.minute,date_it.second) + file_tmp_mpassit = f'{MPASSIT_Dir}{MPASSIT_Prefix}{thistime_mpassit}.nc' + file_tmp_fe = f'{FE_PrntOutDir}{FE_PrntOutPrefix}{thistime_fe}' + files_list_mpassit.append(file_tmp_mpassit) + files_list_fe.append(file_tmp_fe) + date_it = date_it + dt.timedelta(seconds=secInc) + +# Perform a set of conversions on the new netcdf files +for idx,file in enumerate(files_list_mpassit): + print("Processing "+file) + ds=xr.open_dataset(file) + ds['PHB']=grav*ds['PHB'] + ds['PH']=grav*ds['PH'] + ds['ALT']=(rd/p1000mb)*(300.0+ds['T'])*(1.0+rvovrd*ds['QVAPOR'])*(((ds['P']+ds['PB'])/p1000mb)**cvpm) + ds.to_netcdf(files_list_fe[idx]) # rewrite to netcdf + ds.close() diff --git a/scripts/python_utilities/coupler/simgrid.json b/scripts/python_utilities/coupler/simgrid.json index 697270e6..ddaabe22 100644 --- a/scripts/python_utilities/coupler/simgrid.json +++ b/scripts/python_utilities/coupler/simgrid.json @@ -1,5 +1,5 @@ { - "name_dom": "FortCollinsCO", + "name_dom": "FEdomainName", "FE_ref_GIS_nc": "geospec_file.nc", "FE_params_file": "FE_parameters_file.in", "center_lat": 40.5948, @@ -7,5 +7,8 @@ "urban_opt": 0, "FE_new_nc_path": "/path_simgrid/", "name_dom_add": "", + "urban_heatRedis_opt": 0, + "landcover_table": "/path/to/LandCoverMetadata_NLCD16.csv", + "topo_average_opt": 0, "save_plot_opt": 1 } diff --git a/scripts/python_utilities/post-processing/FEtowersToNetCDF.py b/scripts/python_utilities/post-processing/FEtowersToNetCDF.py new file mode 100644 index 00000000..2db3d525 --- /dev/null +++ b/scripts/python_utilities/post-processing/FEtowersToNetCDF.py @@ -0,0 +1,269 @@ +import os, sys +import struct +import numpy as np +import numpy.matlib +import xarray as xr +import pandas as pd +import time +import warnings +import gc +import json +import argparse +from pathlib import Path +def parse_args(): + """ parse the command line arguments """ + + parser = argparse.ArgumentParser() + parser.add_argument("-f", "--file", required=True, help="JSON file with coupler parameter settings") + args = parser.parse_args() + return args + +def get_params_FE(FE_params_file): + + FE_params_dict = {} + n_header = 1 + f = open(FE_params_file,'r') + data = f.readlines() + f.close + row_len = len(data) - n_header + col_len = len(data[n_header].split()) + for rr in range(n_header,row_len+n_header): + row_rr = data[rr] + if (row_rr[0]=='#'): + continue + varname_rr = row_rr.split('=')[0].split() + varval_rr = row_rr.split('=')[1].split('#')[0].split() + FE_params_dict[varname_rr[0]] = varval_rr + + return FE_params_dict + +################## main() ################################################################################ + +######################################## +### Parse the command line arguments ### +######################################## +args = parse_args() + +######################################################### +### Read the json file of converter script parameters ### +######################################################### +with open(args.file) as file: + params = json.loads(file.read()) + +runPath = params["runPath"] +FEparamsFile = params["FEparamsFile"] +outputFileName = params["outputFileName"] +startStep = params["startStep"] +endStep = params["endStep"] + +#Open and parse the FEparamsFile if it exists, else eit +if(Path(f"{runPath}/{FEparamsFile}").exists()): + FE_params = get_params_FE(f"{runPath}/{FEparamsFile}") +elif(Path(f"{FEparamsFile}").exists()): + FE_params = get_params_FE(f"{FEparamsFile}") +else: + sys.exit(f"ERROR: Could not find either {FEparamsFile} or {runPath}/{FEparamsFile}.\nExiting Now!") # Exit with error message + +#Gather configuration parameters from the FEtowerSpecsFile and FEparamsFile +batchSteps=int(FE_params['NtBatch'][0]) +Nz=int(FE_params['Nz'][0]) +towerPath=f"{runPath}/{FE_params['towerPath'][0]}" + +FEtowerSpecsFile = FE_params["towerSpecsFile"][0] +#Open the FEtowerSpecsFile if it exists, else exit +if(Path(f"{runPath}/{FEtowerSpecsFile}").exists()): + ds_towSpecs = xr.open_dataset(f"{runPath}/{FEtowerSpecsFile}") +elif(Path(f"{FEtowerSpecsFile}").exists()): + ds_towSpecs = xr.open_dataset(f"{FEtowerSpecsFile}") +else: + sys.exit(f"ERROR: Could not find either {FEtowerSpecsFile} or {runPath}/{FEtowerSpecsFile}.\nExiting Now!") # Exit with error message +#Read the towerSpecs to determine the number of tower in the run +nTowers = ds_towSpecs.sizes['nProfs'] + +#Build the list of towerFlds and towerSurfFlds contained in the raw binary files by parsing the controlling switches in the FEparamsFile +nTowerVars = 5 #The minimum number of tower variable profiles that should have been written +towerFldNames = ['rho', 'u', 'v', 'w', 'theta'] +nSurfVars = 6 #The minimum number of surf variables that should have been written +towerSurfFldNames = ['z0m', 'z0t', 'tskin', 'fricVel', 'invObLen', 'htFlux'] +nTKE=int(FE_params['TKESelector'][0])*int(FE_params['turbulenceSelector'][0]) +for iTKE in range(nTKE): + nTowerVars += 1 + towerFldNames.append(f"TKE_{iTKE}") +nmoist=int(FE_params['moistureNvars'][0])*int(FE_params['moistureSelector'][0]) +for imoist in range(nmoist): + nTowerVars += 1 + if imoist == 0: + towerFldNames.append("qv") + nSurfVars = nSurfVars + 2 + towerSurfFldNames.extend(['qskin', 'qFlux']) + else: + towerFldNames.append("ql") +if "NhydroAuxScalars" in FE_params: + NhydroAuxScalars = int(FE_params['NhydroAuxScalars'][0]) + for iAuxSc in range(NhydroAuxScalars): + nTowerVars += 1 + towerFldNames.append(f"AuxScalar_{iAuxSc}") +if(int(FE_params['hydroSubGridWrite'][0]) > 0): + nTaus = 9 + nTowerVars = nTowerVars + nTaus + towerFldNames.extend(['Tau11', 'Tau21', 'Tau31', 'Tau32', 'Tau22', 'Tau33', 'TauTH1', 'TauTH2', 'TauTH3']) + for imoist in range(nmoist): + nTowerVars = nTowerVars+3 + if imoist == 0: + towerFldNames.extend(['TauQv1','TauQv2','TauQv3']) + else: + towerFldNames.extend(['TauQl1','TauQl2','TauQl3']) + +#Determine the number of batches the user requested to consolidate into a NetCDF file +nBatches = np.int32((endStep-startStep)/batchSteps) + +#Summarize the intended consolidation parameters +print(f"Attempting to read raw tower files from {towerPath}.") +print(f"Consolidating {nBatches} of {batchSteps} timestep-instances into a single timeseries.") +print(f"Expecting nTowerVars = {nTowerVars}, nSurfVars = {nSurfVars}") + +#Preallocate numpy arrays for all the binary data +towerZ = np.zeros((nTowers,Nz),dtype=np.float32) +towerY = np.zeros((nTowers),dtype=np.float32) +towerX = np.zeros((nTowers),dtype=np.float32) +towerElev = np.zeros((nTowers),dtype=np.float32) +towerSeaMask = np.zeros((nTowers),dtype=np.float32) +towerYoffset = np.zeros((nTowers),dtype=np.float64) +towerXoffset = np.zeros((nTowers),dtype=np.float64) +towerTimes = np.zeros((nBatches*batchSteps+1),dtype=np.float32) +towerData = np.zeros((nTowers,nBatches*batchSteps+1,nTowerVars,Nz),dtype=np.float32) +towerSurfData = np.zeros((nTowers,nBatches*batchSteps+1,nSurfVars),dtype=np.float32) + +########################################################################### +### Read the tower_ic_*.0 files (Initial/static conditions) +########################################################################### +iStep=0 +for itower in range(0,nTowers): + thisFile=f"{towerPath}/tower_ic_{itower}.{iStep}" + flength = os.stat(thisFile).st_size + try: + with open(thisFile, mode='rb') as f: + while(f.tell() < flength): #while the filepointer is not at the end of the binary file + #----- Tower static data + towerNz=struct.unpack("i", f.read(4))[0] + towerZ[itower,:]=np.frombuffer(f.read(towerNz*4),dtype=np.float32) + towerY[itower]=np.frombuffer(f.read(4),dtype=np.float32)[0] + towerX[itower]=np.frombuffer(f.read(4),dtype=np.float32)[0] + towerElev[itower]=np.frombuffer(f.read(4),dtype=np.float32)[0] + if "surflayer_offshore" not in FE_params: + towerSeaMask[itower]=np.frombuffer(f.read(4),dtype=np.float32)[0] + else: + if(int(FE_params['surflayer_offshore'][0]) > 0): + towerSeaMask[itower]=np.frombuffer(f.read(4),dtype=np.float32)[0] + towerYoffset[itower]=np.frombuffer(f.read(8),dtype=np.float64)[0] + towerXoffset[itower]=np.frombuffer(f.read(8),dtype=np.float64)[0] + + #----- Tower profile data + ## Read and parse the number elements per tower instance (single timestep) + towerInstanceSize=struct.unpack("i", f.read(4))[0] + ## Read and parse the number instances in this file (batch of timesteps) + batchSize=struct.unpack("i", f.read(4))[0] + ## Read the full set of batch tower time values in this file + if itower == 0: + towerTimes[iStep]=np.frombuffer(f.read(4),dtype=np.float32)[0] + else: + towerTmpTimes=np.frombuffer(f.read(4),dtype=np.float32) + ## Read the full set of batch tower instances in this file + towerData[itower,iStep,:,:]=np.frombuffer(f.read(towerInstanceSize*4),dtype=np.float32).reshape((nTowerVars,Nz)) + + #---- Tower surf data + ## Read and parse the number elements per tower instance (single timestep) + towerSurfInstanceSize=struct.unpack("i", f.read(4))[0] + ## Read and parse the number instances in this file (batch of timesteps) + batchSize=struct.unpack("i", f.read(4))[0] + ## Read the full set of batch tower time values in this file + towerSurfTimes=np.frombuffer(f.read(batchSize*4),dtype=np.float32) + ## Read the full set of batch tower instances in this file + towerSurfData[itower,iStep,:]=np.frombuffer(f.read(towerSurfInstanceSize*4),dtype=np.float32).reshape((towerSurfInstanceSize)) + + except IOError: + print(f"Error While Opening the file: {thisFile}") + + finally: + if f: # Check if f was successfully assigned a file object + f.close() + # The file is closed here + +########################################################################### +### Reads all of the subsequent batches of timesteps +########################################################################### +if(startStep == 0): + iStep = (startStep+1) +else: + iStep = (startStep) +while iStep < (endStep+1)-batchSteps+1: + print(f"Reading at iStep = {iStep}...") + for itower in range(nTowers): + thisFile=f"{towerPath}/tower_{itower}.{(iStep-1)}" + thisSurfFile=f"{towerPath}/tower_sv_{itower}.{(iStep-1)}" + flength = os.stat(thisFile).st_size + try: + with open(thisFile, mode='rb') as f: + while(f.tell() < flength): #while the filepointer is not at the end of the binary file + ## Read and parse the number elements per tower instance (single timestep) + towerInstanceSize=struct.unpack("i", f.read(4))[0] + ## Read and parse the number instances in this file (batch of timesteps) + batchSize=struct.unpack("i", f.read(4))[0] + ## Read the full set of batch tower time values in this file + if itower == 0: + towerTimes[iStep:(iStep+batchSteps)]=np.frombuffer(f.read(batchSize*4),dtype=np.float32) + else: + towerTmpTimes=np.frombuffer(f.read(batchSize*4),dtype=np.float32) + ## Read the full set of batch tower instances in this file + towerData[itower,iStep:(iStep+batchSteps),:,:]=np.frombuffer(f.read(batchSize*towerInstanceSize*4),dtype=np.float32).reshape((batchSize,nTowerVars,Nz)) + + except IOError: + print('Error While Opening the file: {:s}'.format(thisFile)) + finally: + if f: # Check if f was successfully assigned a file object + f.close() + # The file is closed here + + flength = os.stat(thisSurfFile).st_size + try: + with open(thisSurfFile, mode='rb') as f: + while(f.tell() < flength): #while the filepointer is not at the end of the binary file + ## Read and parse the number elements per tower instance (single timestep) + towerSurfInstanceSize=struct.unpack("i", f.read(4))[0] + ## Read and parse the number instances in this file (batch of timesteps) + batchSize=struct.unpack("i", f.read(4))[0] + ## Read the full set of batch tower time values in this file + towerSurfTimes=np.frombuffer(f.read(batchSize*4),dtype=np.float32) + ## Read the full set of batch tower instances in this file + towerSurfData[itower,iStep:(iStep+batchSteps),:]=np.frombuffer(f.read(batchSize*towerSurfInstanceSize*4),dtype=np.float32).reshape((batchSize,towerSurfInstanceSize)) + + except IOError: + print('Error While Opening the file: {:s}'.format(thisFile)) + finally: + if f: # Check if f was successfully assigned a file object + f.close() + # The file is closed here + iStep = iStep + batchSteps + +ds = xr.Dataset() +ds['time'] = xr.DataArray(towerTimes,dims=['time']) +for iVar in range(nTowerVars): + ds[towerFldNames[iVar]] = xr.DataArray(towerData[:,:,iVar,:],dims=['towerID','time','zIndex']) + if towerFldNames[iVar] in ['u', 'v', 'w', 'theta', 'TKE_0', 'TKE_1', 'qv', 'ql', 'qr', + 'Tau11','Tau21','Tau31','Tau32','Tau22','Tau33', + 'TauTH1','TauTH2','TauTH3', + 'TauQv1','TauQv2','TauQv3', + 'TauQv1','TauQl2','TauQl3',]: + ds[towerFldNames[iVar]] = ds[towerFldNames[iVar]]/ds['rho'] +for iVar in range(nSurfVars): + ds[towerSurfFldNames[iVar]] = xr.DataArray(towerSurfData[:,:,iVar],dims=['towerID','time']) + +ds['z'] = xr.DataArray(towerZ,dims=['towerID','zIndex']) +ds['y'] = xr.DataArray(towerY,dims=['towerID']) +ds['x'] = xr.DataArray(towerX,dims=['towerID']) +ds['elevation'] = xr.DataArray(towerElev,dims=['towerID']) +ds['SeaMask'] = xr.DataArray(towerSeaMask,dims=['towerID']) +ds['yOffset'] = xr.DataArray(towerYoffset,dims=['towerID']) +ds['xOffset'] = xr.DataArray(towerXoffset,dims=['towerID']) + +ds.to_netcdf(f"{towerPath}/{outputFileName}") diff --git a/scripts/python_utilities/post-processing/field_attributes.json b/scripts/python_utilities/post-processing/field_attributes.json index 8fffaab6..03fa39a2 100644 --- a/scripts/python_utilities/post-processing/field_attributes.json +++ b/scripts/python_utilities/post-processing/field_attributes.json @@ -24,12 +24,12 @@ "ql": ["g kg-1", "Cloud liquid water mixing ratio", "cloud_liquid_water_mixing_ratio"], "fricVel": ["m s-1", "Surface friction velocity", "surface_friction_velocity"], "htFlux": ["K m s-1", "Surface sensible heat flux", "surface_upward_sensible_heat_flux"], - "qFlux": ["kg kg-1 m s-1", "Surface latent heat flux", "surface_upward_latent_heat_flux"], + "qFlux": ["g kg-1 m s-1", "Surface latent heat flux", "surface_upward_latent_heat_flux"], "tskin": ["K", "Surface skin temperature", "surface_temperature"], - "qskin": ["kg kg-1", "Surface skin water vapor mixing ratio", null], + "qskin": ["g kg-1", "Surface skin water vapor mixing ratio", null], "z0m": ["m", "Roughness length for momentum", "surface_roughness_length_for_momentum_in_air"], "z0t": ["m", "Roughness length for heat", "surface_roughness_length_for_heat_in_air"], - "invOblen": ["m-1", "Inverse Obukhov length", "atmosphere_boundary_layer_thickness"], + "invOblen": ["m-1", "Inverse Obukhov length", null], "cellpert_amp": ["K", "Cell perturbation amplitude", null], "cellpert_nts": ["-", "Cell perturbation number of time steps", null], "cellpert_ktop": ["-", "Cell perturbation top grid level", null], diff --git a/scripts/python_utilities/post-processing/towers.json b/scripts/python_utilities/post-processing/towers.json new file mode 100644 index 00000000..80b21999 --- /dev/null +++ b/scripts/python_utilities/post-processing/towers.json @@ -0,0 +1,7 @@ +{ + "runPath": "INSERT_PATH_TO_YOUR_RUN_DIRECTORY/", + "FEparamsFile": "FE_params.in", + "outputFileName": "Alltowers.nc", + "startStep": 0, + "endStep": 300000 +}