Updates to read land/sea mask (#1656)

* changes to read land/sea mask * changes to read land/sea mask
erf-model · Jun 21, 2024 · 5ce79bb · 5ce79bb
1 parent 966cdb7
commit 5ce79bb
Showing 1 changed file with 21 additions and 30 deletions.
diff --git a/Source/ERF_read_waves.cpp b/Source/ERF_read_waves.cpp
@@ -12,12 +12,10 @@ ERF::read_waves (int lev)
 {
     for ( MFIter mfi(*Hwave_onegrid[lev],false); mfi.isValid(); ++mfi)
     {
-         //auto my_H_ptr = Hwave[lev]->array(mfi);
-         //auto my_L_ptr = Lwave[lev]->array(mfi);
+
          const auto & bx = mfi.validbox();
 
-        amrex::Print() <<  " HERE " << bx << std::endl;
-         // How to declare my_H_ptr directly?
+         amrex::Print() <<  " HERE " << bx << std::endl;
          amrex::Array4<Real> my_H_arr = Hwave_onegrid[lev]->array(mfi);
          amrex::Array4<Real> my_L_arr = Lwave_onegrid[lev]->array(mfi);
 
@@ -35,17 +33,17 @@ ERF::read_waves (int lev)
          }
 
 
-         int nx=2147483647; // sanity check
+         int nx=2147483647;
          int ny=2147483647; // sanity check
 
-         //JUST RECV
+         //JUST RECEIVED
          if (amrex::MPMD::MyProc() == this_root) {
-             if (rank_offset == 0) // the first program
+             if (rank_offset == 0) // First program
              {
                  MPI_Recv(&nx, 1, MPI_INT, other_root, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                  MPI_Recv(&ny, 1, MPI_INT, other_root, 7, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
              }
-             else // the second program
+             else // Second program
              {
                  MPI_Recv(&nx, 1, MPI_INT, other_root, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                  MPI_Recv(&ny, 1, MPI_INT, other_root, 6, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
@@ -54,17 +52,9 @@ ERF::read_waves (int lev)
              ParallelDescriptor::Bcast(&nx, 1);
              ParallelDescriptor::Bcast(&ny, 1);
          }
+
          if((nx)*(ny) > 0) {
              int nsealm = (nx)*ny;
-             Print()<<nsealm<<std::endl;
-             Print()<<" NX = " << nx<<std::endl;
-             Print()<<" Ny = " << ny<<std::endl;
-             Print()<<" BX  NPTS = " << bx.numPts() <<std::endl;
-             Print()<<" BX       = " << bx <<std::endl;
-             // Print()<<" FAB NPTS = " << Hwave_onegrid[0]->array(mfi).size() <<std::endl;
-             Print()<<" FAB NPTS = " << (*Hwave_onegrid[0])[0].box().numPts()<<std::endl;
-             Print()<<" FAB      = " << (*Hwave_onegrid[0])[0].box() <<std::endl;
-             // AMREX_ALWAYS_ASSERT_WITH_MESSAGE((nx)*ny <= bx.numPts() || bx.numPts() == 0, "total number of points being filled exceeds the size of the current box\n");
 
              if (amrex::MPMD::MyProc() == this_root) {
                  if (rank_offset == 0) // the first program
@@ -78,30 +68,31 @@ ERF::read_waves (int lev)
                      MPI_Recv(my_L_ptr, nsealm, MPI_DOUBLE, other_root, 4, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                  }
              }
-             amrex::Print()<<"Just recieved "<<nsealm<<"as a double*"<<std::endl;
 
-             if(bx.contains(IntVect(192,2,0))) {
-                 std::cout<<my_H_arr(192,92,0)<<std::endl;
-                 std::cout<<my_L_arr(192,92,0)<<std::endl;
-                 std::cout<<my_H_ptr[192-0+(92-0)*193]<<std::endl;
-                 std::cout<<my_L_ptr[192-0+(92-0)*193]<<std::endl;
-             }
-             amrex::AllPrintToFile("output_HS_cpp.txt")<<FArrayBox(my_H_arr)<<std::endl;
-             amrex::AllPrintToFile("output_L_cpp.txt")<<FArrayBox(my_L_arr)<<std::endl;
-         } else {
-             finished_wave = true;
          }
     }
+
+
     //May need to be Redistribute
     //    ParallelCopy(Hwave[lev],Hwave_onegrid[lev],0,0,1,0);
     Hwave[lev]->ParallelCopy(*Hwave_onegrid[lev]);
     Hwave[lev]->FillBoundary(geom[lev].periodicity());
     Lwave[lev]->ParallelCopy(*Lwave_onegrid[lev]);
     Lwave[lev]->FillBoundary(geom[lev].periodicity());
     amrex::Print() << "HWAVE BOX " << (*Hwave[lev])[0].box() << std::endl;
+    for (MFIter mfi(*Hwave[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+        Box bx = mfi.tilebox();
+        const Array4<Real const>& Hwave_arr = Hwave[lev]->const_array(mfi);
+        const Array4<int>& Lmask_arr = lmask_lev[lev][0]->array(mfi);
+        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k){
+            if (Hwave_arr(i,j,k)<0) {
+                Lmask_arr(i,j,k) = 1;
+             } else {
+                Lmask_arr(i,j,k) = 0;
+            }
+        });
+    }
 
-    print_state(*Hwave[lev],IntVect(103,-3,0),0,IntVect(3,3,0));
-    print_state(*Hwave[lev],IntVect(103,88,0),0,IntVect(3,3,0));
 }
 #endif