C++ and Fortran implementation of Legion and Kokkos plus with Cuda driver as well

src/cpp/CanopyHydrology_Interception.cpp

@@ -10,7 +10,6 @@
 using namespace std;
 
 namespace ELM {
-
 void CanopyHydrology_Interception(double dtime,
         const double& forc_rain,
         const double& forc_snow,

tests/tests_c_legion/Makefile

-OBJECT        = ../../src/
-KERNEL_LANG	  = cc
-SRCDIR	      =	$(OBJECT)$(KERNEL_LANG)
-
-include $(OBJECT)config/Makefile.config
-
-INC_FLAGS = -I$(AMANZI_TPLS_DIR)/include -I../../src/cc/ -I../test_c_legion
-
-TESTS = test_CanopyHydrology_kern1_single \
-        test_CanopyHydrology_kern1_multiple \
-        test_CanopyHydrology_module
-
-EXEC_TESTS = CanopyHydrology_kern1_single \
-             CanopyHydrology_kern1_multiple \
-             CanopyHydrology_module
-
-
-.PHONY: links library test
-
-default: all
-
-all: links library $(TESTS)
-
-test: $(EXEC_TESTS)
-	python ../compare_to_gold.py $(TESTS)
-
-
-CanopyHydrology_kern1_single: test_CanopyHydrology_kern1_single
-	./test_CanopyHydrology_kern1_single &> test_CanopyHydrology_kern1_single.stdout
-
-CanopyHydrology_kern1_multiple: test_CanopyHydrology_kern1_multiple
-	./test_CanopyHydrology_kern1_multiple &> test_CanopyHydrology_kern1_multiple.stdout
-
-CanopyHydrology_module: test_CanopyHydrology_module
-	./test_CanopyHydrology_module &> test_CanopyHydrology_module.stdout
-
-
-test_%: %.cc.o readers.hh utils.hh library
-	$(CXX) -o $@ $< $(CC_LD_FLAGS)
-
-
-clean:
-	@$(ELM_CLEAN)
-	$(RM) test_*
-
-allclean:
-	@$(ELM_CLEAN)
-	$(RM) test_*
-	$(MAKE) -C $(OBJECT) allclean
-
-links:
-	@echo "making in links"
-	$(MAKE) -C ../links links
-
-library:
-	$(MAKE) -C $(OBJECT) all
\ No newline at end of file

tests/tests_kokkos_c/CanopyHydrology_kern1_multiple.cpp

@@ -13,9 +13,153 @@
 #include <Kokkos_Core.hpp>
 #include "utils.hh"
 #include "readers.hh"
-#include "CanopyHydrology_cpp.hh" 
+#include "landunit_varcon.h"
+#include "column_varcon.h" 
 
 
+namespace ELM {
+void CanopyHydrology_Interception(double dtime,
+        const double& forc_rain,
+        const double& forc_snow,
+        const double& irrig_rate,
+        const int& ltype, const int& ctype,
+        const bool& urbpoi, const bool& do_capsnow,
+        const double& elai, const double& esai,
+        const double& dewmx, const int& frac_veg_nosno,
+        double& h2ocan,
+        int n_irrig_steps_left, //fix it
+        double& qflx_prec_intr,
+        double& qflx_irrig,
+        double& qflx_prec_grnd,
+        double& qflx_snwcp_liq,
+        double& qflx_snwcp_ice,
+        double& qflx_snow_grnd_patch,
+        double& qflx_rain_grnd)
+
+ {  
+
+  
+      double  fpi, xrun, h2ocanmx   ;
+      double  qflx_candrip, qflx_through_snow, qflx_through_rain ;
+      double  qflx_prec_grnd_snow;
+      double  qflx_prec_grnd_rain ;
+      double  fracsnow ;
+      double  fracrain , forc_irrig;
+
+
+      if (ltype==istsoil || ltype==istwet || urbpoi || ltype==istcrop) {
+
+         qflx_candrip = 0.0      ;
+         qflx_through_snow = 0.0 ;
+         qflx_through_rain = 0.0 ;
+         qflx_prec_intr = 0.0    ;
+         fracsnow = 0.0          ;
+         fracrain = 0.0          ;
+         forc_irrig = 0.0;
+
+
+         if (ctype != icol_sunwall && ctype != icol_shadewall) {
+            if (frac_veg_nosno == 1 && (forc_rain + forc_snow) > 0.0) {
+
+              
+               fracsnow = forc_snow/(forc_snow + forc_rain);
+               fracrain = forc_rain/(forc_snow + forc_rain);
+
+               
+               h2ocanmx = dewmx * (elai + esai);
+
+               
+               fpi = 0.250*(1.0 - exp(-0.50*(elai + esai)));
+
+              
+               qflx_through_snow = forc_snow * (1.0-fpi);
+               qflx_through_rain = forc_rain * (1.0-fpi);
+
+               
+               qflx_prec_intr = (forc_snow + forc_rain) * fpi;
+               
+
+
+               
+               h2ocan = fmax(0.0, h2ocan + dtime*qflx_prec_intr);
+
+               
+               qflx_candrip = 0.0;
+
+               
+               xrun = (h2ocan - h2ocanmx)/dtime;
+
+               
+               if (xrun > 0.0) {
+                  qflx_candrip = xrun;
+                  h2ocan = h2ocanmx;
+               }
+
+            }
+         }
+
+      else if (ltype==istice || ltype==istice_mec) {
+         
+         h2ocan            = 0.0;
+         qflx_candrip      = 0.0;
+         qflx_through_snow = 0.0;
+         qflx_through_rain = 0.0;
+         qflx_prec_intr    = 0.0;
+         fracsnow          = 0.0;
+         fracrain          = 0.0;
+
+      }
+
+      
+
+      if (ctype != icol_sunwall && ctype != icol_shadewall) {
+         if (frac_veg_nosno == 0) {
+            qflx_prec_grnd_snow = forc_snow;
+            qflx_prec_grnd_rain = forc_rain;  }
+         else{
+            qflx_prec_grnd_snow = qflx_through_snow + (qflx_candrip * fracsnow);
+            qflx_prec_grnd_rain = qflx_through_rain + (qflx_candrip * fracrain);
+          }
+      }   
+      else{
+         qflx_prec_grnd_snow = 0.;
+         qflx_prec_grnd_rain = 0.;
+        }
+
+      
+      if (n_irrig_steps_left > 0) {
+         qflx_irrig         = forc_irrig;
+         n_irrig_steps_left = n_irrig_steps_left - 1; }
+      else{
+         qflx_irrig = 0.0;
+        }
+
+      
+      qflx_prec_grnd_rain = qflx_prec_grnd_rain + qflx_irrig;
+
+      
+
+      qflx_prec_grnd = qflx_prec_grnd_snow + qflx_prec_grnd_rain;
+
+      if (do_capsnow) {
+         qflx_snwcp_liq = qflx_prec_grnd_rain;
+         qflx_snwcp_ice = qflx_prec_grnd_snow;
+
+         qflx_snow_grnd_patch = 0.0;
+         qflx_rain_grnd = 0.0;  }
+      else{
+
+         qflx_snwcp_liq = 0.0;
+         qflx_snwcp_ice = 0.0;
+         qflx_snow_grnd_patch = qflx_prec_grnd_snow   ;      //ice onto ground (mm/s)
+         qflx_rain_grnd     = qflx_prec_grnd_rain      ;   //liquid water onto ground (mm/s)
+        }
+
+    }
+  }
+
+}
+
 namespace ELM {
 namespace Utils {
 
@@ -51,7 +195,8 @@ int main(int argc, char ** argv)
   const double dewmx = 0.1;
   double dtime = 1800.0;
 
-  Kokkos::initialize( argc, argv );
+  Kokkos::HostSpace::execution_space::initialize();
+  Kokkos::Cuda::initialize( );//argc, argv );
   {
 
   typedef Kokkos::View<double**>  ViewMatrixType;
@@ -82,8 +227,9 @@ int main(int argc, char ** argv)
   ViewMatrixType::HostMirror h_forc_snow = Kokkos::create_mirror_view( forc_snow );
   ViewMatrixType::HostMirror h_forc_air_temp = Kokkos::create_mirror_view( forc_air_temp );
   const int n_times = ELM::Utils::read_forcing("../links/forcing", n_max_times, 0, n_grid_cells, h_forc_rain, h_forc_snow, h_forc_air_temp);
-
-  ELM::Utils::MatrixForc forc_irrig; forc_irrig = 0.;
+  ViewMatrixType forc_irrig( "forc_irrig", n_max_times,n_grid_cells );
+  ViewMatrixType::HostMirror h_forc_irrig = Kokkos::create_mirror_view( forc_irrig );
+  //ELM::Utils::MatrixForc forc_irrig; forc_irrig = 0.;
   
   // output state by the grid cell
   // auto qflx_prec_intr = std::array<double,n_grid_cells>();
@@ -120,7 +266,7 @@ int main(int argc, char ** argv)
   //auto h2o_can = ELM::Utils::MatrixState(); 
   ViewMatrixType::HostMirror h_h2o_can = Kokkos::create_mirror_view( h2o_can );
   //h_h2o_can = 0.;
-  auto h2o_can1 = ELM::Utils::MatrixState(); 
+  //auto h2o_can1 = ELM::Utils::MatrixState(); 
   // Array<int64_t, 2> a = h_h2o_can;
   //const size_t n0 = h2o_can.extent_0 ();
 
@@ -130,6 +276,7 @@ int main(int argc, char ** argv)
   Kokkos::deep_copy( esai, h_esai);
   Kokkos::deep_copy( forc_rain, h_forc_rain);
   Kokkos::deep_copy( forc_snow, h_forc_snow);
+  Kokkos::deep_copy( forc_irrig, h_forc_irrig);
   Kokkos::deep_copy( forc_air_temp, h_forc_air_temp);
   Kokkos::deep_copy( qflx_prec_intr, h_qflx_prec_intr);
   Kokkos::deep_copy( qflx_irrig, h_qflx_irrig);
@@ -196,6 +343,8 @@ int main(int argc, char ** argv)
 
   }
   }
-  Kokkos::finalize();
+  //Kokkos::finalize();
+  Kokkos::Cuda::finalize();
+  Kokkos::HostSpace::execution_space::finalize();
   return 0;
 }

tests/tests_kokkos_c/CanopyHydrology_module.cpp

@@ -10,11 +10,12 @@
 #include <iomanip>
 #include <numeric>
 #include <Kokkos_Core.hpp>
+#include "landunit_varcon.h"    
+#include "column_varcon.h" 
+#include "clm_varpar.h"         
+#include "clm_varctl.h"   
 #include "utils.hh"
 #include "readers.hh"
-#include "CanopyHydrology_cpp.hh"   
-#include "CanopyHydrology_SnowWater_impl.hh" 
-
 
 namespace ELM {
 namespace Utils {
@@ -35,6 +36,503 @@ using VectorColumnInt = VectorStatic<n_grid_cells,int>;
 } // namespace
 } // namespace
 
+namespace ELM {
+void CanopyHydrology_Interception(double dtime,
+        const double& forc_rain,
+        const double& forc_snow,
+        const double& irrig_rate,
+        const int& ltype, const int& ctype,
+        const bool& urbpoi, const bool& do_capsnow,
+        const double& elai, const double& esai,
+        const double& dewmx, const int& frac_veg_nosno,
+        double& h2ocan,
+        int n_irrig_steps_left, //fix it
+        double& qflx_prec_intr,
+        double& qflx_irrig,
+        double& qflx_prec_grnd,
+        double& qflx_snwcp_liq,
+        double& qflx_snwcp_ice,
+        double& qflx_snow_grnd_patch,
+        double& qflx_rain_grnd)
+
+ {  
+
+  
+      double  fpi, xrun, h2ocanmx   ;
+      double  qflx_candrip, qflx_through_snow, qflx_through_rain ;
+      double  qflx_prec_grnd_snow;
+      double  qflx_prec_grnd_rain ;
+      double  fracsnow ;
+      double  fracrain , forc_irrig;
+
+
+      if (ltype==istsoil || ltype==istwet || urbpoi || ltype==istcrop) {
+
+         qflx_candrip = 0.0      ;
+         qflx_through_snow = 0.0 ;
+         qflx_through_rain = 0.0 ;
+         qflx_prec_intr = 0.0    ;
+         fracsnow = 0.0          ;
+         fracrain = 0.0          ;
+         forc_irrig = 0.0;
+
+
+         if (ctype != icol_sunwall && ctype != icol_shadewall) {
+            if (frac_veg_nosno == 1 && (forc_rain + forc_snow) > 0.0) {
+
+              
+               fracsnow = forc_snow/(forc_snow + forc_rain);
+               fracrain = forc_rain/(forc_snow + forc_rain);
+
+               
+               h2ocanmx = dewmx * (elai + esai);
+
+               
+               fpi = 0.250*(1.0 - exp(-0.50*(elai + esai)));
+
+              
+               qflx_through_snow = forc_snow * (1.0-fpi);
+               qflx_through_rain = forc_rain * (1.0-fpi);
+
+               
+               qflx_prec_intr = (forc_snow + forc_rain) * fpi;
+               
+
+
+               
+               h2ocan = std::max(0.0, h2ocan + dtime*qflx_prec_intr);
+
+               
+               qflx_candrip = 0.0;
+
+               
+               xrun = (h2ocan - h2ocanmx)/dtime;
+
+               
+               if (xrun > 0.0) {
+                  qflx_candrip = xrun;
+                  h2ocan = h2ocanmx;
+               }
+
+            }
+         }
+
+      else if (ltype==istice || ltype==istice_mec) {
+         
+         h2ocan            = 0.0;
+         qflx_candrip      = 0.0;
+         qflx_through_snow = 0.0;
+         qflx_through_rain = 0.0;
+         qflx_prec_intr    = 0.0;
+         fracsnow          = 0.0;
+         fracrain          = 0.0;
+
+      }
+
+      
+
+      if (ctype != icol_sunwall && ctype != icol_shadewall) {
+         if (frac_veg_nosno == 0) {
+            qflx_prec_grnd_snow = forc_snow;
+            qflx_prec_grnd_rain = forc_rain;  }
+         else{
+            qflx_prec_grnd_snow = qflx_through_snow + (qflx_candrip * fracsnow);
+            qflx_prec_grnd_rain = qflx_through_rain + (qflx_candrip * fracrain);
+          }
+      }   
+      else{
+         qflx_prec_grnd_snow = 0.;
+         qflx_prec_grnd_rain = 0.;
+        }
+
+      
+      if (n_irrig_steps_left > 0) {
+         qflx_irrig         = forc_irrig;
+         n_irrig_steps_left = n_irrig_steps_left - 1; }
+      else{
+         qflx_irrig = 0.0;
+        }
+
+      
+      qflx_prec_grnd_rain = qflx_prec_grnd_rain + qflx_irrig;
+
+      
+
+      qflx_prec_grnd = qflx_prec_grnd_snow + qflx_prec_grnd_rain;
+
+      if (do_capsnow) {
+         qflx_snwcp_liq = qflx_prec_grnd_rain;
+         qflx_snwcp_ice = qflx_prec_grnd_snow;
+
+         qflx_snow_grnd_patch = 0.0;
+         qflx_rain_grnd = 0.0;  }
+      else{
+
+         qflx_snwcp_liq = 0.0;
+         qflx_snwcp_ice = 0.0;
+         qflx_snow_grnd_patch = qflx_prec_grnd_snow   ;      //ice onto ground (mm/s)
+         qflx_rain_grnd     = qflx_prec_grnd_rain      ;   //liquid water onto ground (mm/s)
+        }
+
+    }
+  }
+
+}
+
+namespace ELM {
+
+void CanopyHydrology_FracWet(const int& frac_veg_nosno,
+        const double& h2ocan,
+        const double& elai, 
+        const double& esai,
+        const double& dewmx,
+        double& fwet,
+        double& fdry)
+{  
+
+  double vegt, dewmxi ;    
+
+  if (frac_veg_nosno == 1) {
+    if (h2ocan > 0.0) {
+        vegt    = frac_veg_nosno*(elai + esai);
+        dewmxi  = 1.00/dewmx;
+        fwet = std::pow(((dewmxi/vegt)*h2ocan), 2.0/3);
+        fwet = std::min(fwet,1.00);   
+        fdry = (1.0-fwet)*elai/(elai+esai);
+      }
+    else{
+        fwet = 0.0;
+        fdry = 0.0 ;
+      } }
+  else{
+     fwet = 0.0;
+     fdry = 0.0;
+  }
+
+}
+
+}
+
+
+namespace ELM {
+
+template<typename Array_d>
+void CanopyHydrology_SnowWater(const double& dtime,
+        const double& qflx_floodg,
+        const int& ltype,
+        const int& ctype,
+        const bool& urbpoi,
+        const bool& do_capsnow,                            
+        const int& oldfflag,
+        const double& forc_air_temp,
+        const double& t_grnd,
+        const double& qflx_snow_grnd_col,
+        const double& qflx_snow_melt,
+        const double& n_melt,
+        const double& frac_h2osfc,
+        double& snow_depth,
+        double& h2osno,
+        double& integrated_snow,
+        Array_d swe_old,
+        Array_d h2osoi_liq,
+        Array_d h2osoi_ice,
+        Array_d t_soisno,
+        Array_d frac_iceold,
+        int& snow_level,
+        Array_d dz,
+        Array_d z,
+        Array_d zi,
+        int& newnode,
+        double& qflx_floodc,
+        double& qflx_snow_h2osfc,
+        double& frac_sno_eff,
+        double& frac_sno)
+{       
+  
+  
+//parameters
+  double rpi=4.0e0*atan(1.0e0)  ;
+  double tfrz=273.15;
+  double zlnd = 0.010;
+
+  
+  // real(r8), intent(inout), dimension(-nlevsno+1:0)  :: swe_old 
+  // real(r8), intent(inout), dimension(-nlevsno+1:0) :: h2osoi_liq, h2osoi_ice
+  // real(r8), intent(inout), dimension(-nlevsno+1:0)  :: t_soisno, frac_iceold
+  // real(r8), intent(inout), dimension(-nlevsno+1:0)  :: dz, z, zi
+  
+//local variables 
+  double  temp_intsnow, temp_snow_depth, z_avg, fmelt, dz_snowf, snowmelt ;
+  double  newsnow, bifall, accum_factor, fsno_new, smr ;
+  int j ;
+
+//apply gridcell flood water flux to non-lake columns
+  if (ctype != icol_sunwall && ctype != icol_shadewall) {      
+     qflx_floodc = qflx_floodg; }
+  else{
+     qflx_floodc = 0.0;
+  }
+
+//Determine snow height and snow water
+
+//Use Alta relationship, Anderson(1976); LaChapelle(1961),
+//U.S.Department of Agriculture Forest Service, Project F,
+//Progress Rep. 1, Alta Avalanche Study Center:Snow Layer Densification.
+
+  qflx_snow_h2osfc = 0.0;
+//set temporary variables prior to updating
+  temp_snow_depth=snow_depth;
+//save initial snow content
+  for(j = -nlevsno+1; j < snow_level; j++) {
+     swe_old[j] = 0.00;
+  }
+  for(j = snow_level+1; j < 0; j++) {
+     swe_old[j]=h2osoi_liq[j]+h2osoi_ice[j];
+  }
+
+  if (do_capsnow) {
+     dz_snowf = 0.;
+     newsnow = (1. - frac_h2osfc) * qflx_snow_grnd_col * dtime;
+     frac_sno=1.;
+     integrated_snow = 5.e2; 
+  } else {
+
+     if (forc_air_temp > tfrz + 2.) {
+        bifall=50. + 1.7*std::pow((17.0),1.5);
+     } else if (forc_air_temp > tfrz - 15.) {
+        bifall=50. + 1.7*std::pow((forc_air_temp - tfrz + 15.),1.5);
+     } else {
+        bifall=50.;
+     }
+
+     // newsnow is all snow that doesn't fall on h2osfc
+     newsnow = (1. - frac_h2osfc) * qflx_snow_grnd_col * dtime;
+
+     // update integrated_snow
+     integrated_snow = std::max(integrated_snow,h2osno) ; //h2osno could be larger due to frost
+
+     // snowmelt from previous time step * dtime
+     snowmelt = qflx_snow_melt * dtime;
+
+     // set shape factor for accumulation of snow
+     accum_factor=0.1;
+
+     if (h2osno > 0.0) {
+
+        //======================  FSCA PARAMETERIZATIONS  ======================
+        // fsca parameterization based on *changes* in swe
+        // first compute change from melt during previous time step
+        if(snowmelt > 0.) {
+
+           smr=std::min(1.,(h2osno)/(integrated_snow));
+
+           frac_sno = 1. - std::pow((acos(fmin(1.,(2.*smr - 1.)))/rpi),(n_melt)) ;
+
+        }
+
+        // update fsca by new snow event, add to previous fsca
+        if (newsnow > 0.) {
+           fsno_new = 1. - (1. - tanh(accum_factor*newsnow))*(1. - frac_sno);
+           frac_sno = fsno_new;
+
+           // reset integrated_snow after accumulation events
+           temp_intsnow= (h2osno + newsnow) / (0.5*(cos(rpi*std::pow((1.0-std::max(frac_sno,1e-6)),(1.0/n_melt))+1.0))) ;
+           integrated_snow = std::min(1.e8,temp_intsnow) ;
+        }
+
+        //====================================================================
+
+        // for subgrid fluxes
+        if (subgridflag ==1 && ! urbpoi) {