Merge branch 'kripke_miniapp' into 'master'

Kripke miniapp that utilizes the conduit data adaptor.

See merge request !107

CMake/options.cmake

@@ -64,6 +64,7 @@ cmake_dependent_option(ENABLE_VTKM_RENDERING
 option(ENABLE_PARALLEL3D "Enable Parallel3D miniapp" ON)
 option(ENABLE_OSCILLATORS "Enable Oscillators miniapp" ON)
 option(ENABLE_CONDUITTEST "Enable Conduit miniapp" OFF)
+option(ENABLE_KRIPKE "Enable Kripke miniapp" OFF)
 
 message(STATUS "ENABLE_SENSEI=${ENABLE_SENSEI}")
 message(STATUS "ENABLE_PYTHON=${ENABLE_PYTHON}")
@@ -82,3 +83,4 @@ message(STATUS "ENABLE_VTKM_RENDERING=${ENABLE_VTKM_RENDERING}")
 message(STATUS "ENABLE_PARALLEL3D=${ENABLE_PARALLEL3D}")
 message(STATUS "ENABLE_OSCILLATORS=${ENABLE_OSCILLATORS}")
 message(STATUS "ENABLE_CONDUITTEST=${ENABLE_CONDUITTEST}")
+message(STATUS "ENABLE_KRIPKE=${ENABLE_KRIPKE}")

miniapps/CMakeLists.txt

@@ -19,3 +19,10 @@ else()
   message(STATUS "Disabled: Oscillators miniapp.")
 endif()
 
+
+if(ENABLE_KRIPKE)
+  message(STATUS "Enabled: Kripke miniapp")
+  add_subdirectory(kripke)
+else()
+  message(STATUS "Disabled: Kripke miniapp")
+endif()

miniapps/kripke/CMakeLists.txt

+###############################################################################
+# Copyright (c) 2015-2018, Lawrence Livermore National Security, LLC.
+# 
+# Produced at the Lawrence Livermore National Laboratory
+# 
+# LLNL-CODE-716457
+# 
+# All rights reserved.
+# 
+# This file is part of Ascent. 
+# 
+# For details, see: http://ascent.readthedocs.io/.
+# 
+# Please also read ascent/LICENSE
+# 
+# Redistribution and use in source and binary forms, with or without 
+# modification, are permitted provided that the following conditions are met:
+# 
+# * Redistributions of source code must retain the above copyright notice, 
+#   this list of conditions and the disclaimer below.
+# 
+# * Redistributions in binary form must reproduce the above copyright notice,
+#   this list of conditions and the disclaimer (as noted below) in the
+#   documentation and/or other materials provided with the distribution.
+# 
+# * Neither the name of the LLNS/LLNL nor the names of its contributors may
+#   be used to endorse or promote products derived from this software without
+#   specific prior written permission.
+# 
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED. IN NO EVENT SHALL LAWRENCE LIVERMORE NATIONAL SECURITY,
+# LLC, THE U.S. DEPARTMENT OF ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY
+# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
+# DAMAGES  (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+# OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, 
+# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
+# POSSIBILITY OF SUCH DAMAGE.
+# 
+###############################################################################
+
+
+###############################################################################
+#
+# Kripke CMake Build for Ascent
+#
+###############################################################################
+
+
+set(KRIPKE_SOURCES
+    Kripke/Directions.cpp  
+    Kripke/Grid.cpp  
+    Kripke/Kernel.cpp  
+    Kripke/Layout.cpp
+    Kripke/Subdomain.cpp  
+    Kripke/Sweep_Solver.cpp
+    Kripke/ParallelComm.cpp  
+    Kripke/Timing.cpp  
+    
+    Kripke/Kernel/Kernel_3d_GDZ.cpp
+    Kripke/Kernel/Kernel_3d_DGZ.cpp
+    Kripke/Kernel/Kernel_3d_ZDG.cpp
+    Kripke/Kernel/Kernel_3d_DZG.cpp
+    Kripke/Kernel/Kernel_3d_ZGD.cpp
+    Kripke/Kernel/Kernel_3d_GZD.cpp
+  
+    Kripke/ParallelComm/BlockJacobiComm.cpp
+    Kripke/ParallelComm/SweepComm.cpp 
+    tools/kripke.cpp
+    tools/testKernels.cpp
+
+    bridge.cpp
+)
+
+# kripke reqs mpi
+if(MPI_FOUND)
+    # copy over example command line
+    configure_file(run_kripke_simple_example.sh 
+                   ${CMAKE_CURRENT_BINARY_DIR}/run_kripke_simple_example.sh
+                   COPYONLY)
+
+    if(ENABLE_OPENMP)
+       set(kripke_openmp_flags "-DKRIPKE_USE_OPENMP")
+    else()
+       set(kripke_openmp_flags "")
+    endif()
+
+    include_directories(.)
+    include_directories("tools")
+
+
+    set(libs mpi sconduit sensei)
+    if(ENABLE_OPENMP)
+           list(APPEND libs openmp)
+    endif()
+    
+    add_executable(kripke_par ${KRIPKE_SOURCES})
+    target_link_libraries(kripke_par ${libs})
+
+
+    # install target for kripke mpi
+    install(TARGETS kripke_par
+            RUNTIME DESTINATION bin
+    )
+
+#    install(FILES ASCENT_README.md 
+#                  ascent_actions.json
+#                  ascent_options.json
+#            DESTINATION examples/proxies/kripke)
+
+endif()
+
+
+
+
+
+

miniapps/kripke/Kripke.h

+/**
+ * This is the main header file for the Kripke Mini-App.
+ */
+#ifndef KRIPKE_H__
+#define KRIPKE_H__
+
+#include<string>
+#include<vector>
+#include<stdio.h>
+#include<cmath>
+#include<strings.h>
+
+// Make sure that there's openmp support, otherwise error out
+#if KRIPKE_USE_OPENMP
+#ifndef _OPENMP
+#error "OpenMP selected for build, but OpenMP is not available"
+#endif
+#endif
+
+// Forward Decl
+struct Grid_Data;
+
+#define KRESTRICT __restrict__
+
+
+// In Kripke/Sweep_Solver.cpp
+int SweepSolver(Grid_Data *grid_data, bool block_jacobi, const std::string& file);
+void SweepSubdomains (std::vector<int> subdomain_list, Grid_Data *grid_data, bool block_jacobi);
+
+/**
+ * Tags for choosing which data nesting to be chosen
+ */
+enum Nesting_Order {
+  // Nestings for Psi and Phi
+  // D referes to directions OR moments, depending on context
+  NEST_DGZ,
+  NEST_DZG,
+  NEST_GDZ,
+  NEST_GZD,
+  NEST_ZDG,
+  NEST_ZGD
+};
+
+
+/**
+  Tags for which parallel algorithm to use.
+*/
+enum ParallelMethod {
+  PMETHOD_SWEEP,
+  PMETHOD_BJ
+};
+
+/**
+ * Converts a nesting tag to a human-readable string.
+ */
+inline std::string nestingString(Nesting_Order nesting){
+  switch(nesting){
+    case NEST_DGZ: return("DGZ");
+    case NEST_DZG: return("DZG");
+    case NEST_GDZ: return("GDZ");
+    case NEST_GZD: return("GZD");
+    case NEST_ZDG: return("ZDG");
+    case NEST_ZGD: return("ZGD");
+  }
+  return("UNKNOWN");
+}
+
+/**
+ * Converts a string (eg. from command line) to a nesting tag.
+ */
+inline Nesting_Order nestingFromString(std::string const &str){
+  for(int i = 0;i < 6;++ i){
+    if(!strcasecmp(str.c_str(), nestingString((Nesting_Order)i).c_str())){
+      return (Nesting_Order)i;
+  }
+ }
+  return (Nesting_Order)-1;
+}
+
+
+/**
+ * Compares two vectors for differences.
+ * Used in testing suite.
+ */
+inline bool compareVector(std::string const &name,
+    std::vector<double> const &a,
+    std::vector<double> const &b, double tol, bool verbose){
+
+  if(a.size() != b.size()){
+    if(verbose){
+      printf("Vectors are different lengths: %ld, %ld\n",
+          (long)a.size(), (long)b.size());
+    }
+    return true;
+  }
+
+  bool is_diff = false;
+  for(size_t i = 0;i < a.size();++i){
+    if(std::abs(a[i]-b[i]) > tol){
+      is_diff = true;
+      if(verbose){
+        printf("%s[%d]:%e, %e [%e]\n",
+            name.c_str(), (int)i,
+            a[i], b[i], std::abs(a[i]-b[i]));
+        is_diff = true;
+      }
+      else{
+        break;
+      }
+    }
+  }
+
+  return is_diff;
+}
+
+/**
+ * Compares two scalars for differences.
+ * Used in testing suite.
+ */
+inline bool compareScalar(std::string const &name,
+    double a, double b, double tol, bool verbose){
+
+  if(std::abs(a-b) > tol){
+    if(verbose){
+      printf("%s:%e, %e [%e]\n",
+          name.c_str(),
+          a, b, std::abs(a-b));
+    }
+    return true;
+  }
+  return false;
+}
+
+#endif
+

miniapps/kripke/Kripke/CMakeLists.txt

+
+set(KRIPKE_SOURCE ${KRIPKE_SOURCE}
+  Kripke/Directions.cpp  
+  Kripke/Grid.cpp  
+  Kripke/Kernel.cpp  
+  Kripke/Layout.cpp
+  Kripke/Subdomain.cpp  
+  Kripke/Sweep_Solver.cpp
+  Kripke/ParallelComm.cpp  
+  Kripke/Timing.cpp  
+    
+  Kripke/Kernel/Kernel_3d_GDZ.cpp
+  Kripke/Kernel/Kernel_3d_DGZ.cpp
+  Kripke/Kernel/Kernel_3d_ZDG.cpp
+  Kripke/Kernel/Kernel_3d_DZG.cpp
+  Kripke/Kernel/Kernel_3d_ZGD.cpp
+  Kripke/Kernel/Kernel_3d_GZD.cpp
+  
+  Kripke/ParallelComm/BlockJacobiComm.cpp
+  Kripke/ParallelComm/SweepComm.cpp 
+)
+                                                                                                          
+set(KRIPKE_SOURCE ${KRIPKE_SOURCE} PARENT_SCOPE) 

miniapps/kripke/Kripke/Directions.cpp

+#include <Kripke/Directions.h>
+#include <Kripke/Grid.h>
+#include <Kripke/Input_Variables.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <float.h>
+#include <algorithm>
+
+namespace {
+  /*
+    GaussLegendre returns the n point Gauss-Legendre quadrature rule for
+    the integral between x1 and x2.
+  */
+  void GaussLegendre(double x1, double x2, std::vector<double> &x,
+      std::vector<double> &w, double eps)
+  {
+    int n = x.size();
+    int m, j, i;
+    double z1, z, xm, xl, pp, p3, p2, p1;
+
+    m=(n+1)/2;
+    xm=0.5*(x2+x1);
+    xl=0.5*(x2-x1);
+    for(i=1; i<=m; i++){
+      z=cos(M_PI*(i-0.25)/(n+0.5));
+      do {
+        p1=1.0;
+        p2=0.0;
+        for(j=1; j<=n; j++){
+          p3=p2;
+          p2=p1;
+          p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
+        }
+        pp=n*(z*p1-p2)/(z*z-1.0);
+        z1=z;
+        z=z1-p1/pp;
+      } while(fabs(z-z1) > eps);
+      x[i-1]=xm-xl*z;
+      x[n-i]=xm+xl*z;
+      w[i-1]=2.0*xl/((1.0-z*z)*pp*pp);
+
+      w[n-i]=w[i-1];
+    }
+  }
+
+
+  bool dirSortFcn(Directions const &a, Directions const &b){
+    return b.octant < a.octant;
+  }
+}
+
+/**
+ * Initializes the quadrature set information for a Grid_Data object.
+ * This guarantees that each <GS,DS> pair have a single originating octant.
+ */
+void InitDirections(Grid_Data *grid_data, Input_Variables *input_vars)
+{
+  std::vector<Directions> &directions = grid_data->directions;
+
+  // Get set description from user
+  int num_directions_per_octant = input_vars->num_dirs_per_dirset *
+                                  input_vars->num_dirsets_per_octant;
+  int num_directions = 8*num_directions_per_octant;
+
+  // allocate storage
+  directions.resize(num_directions);
+
+  // Are we running a REAL quadrature set?
+  int num_polar = input_vars->quad_num_polar;
+  int num_azimuth = input_vars->quad_num_azimuthal;
+
+  std::vector<double> polar_cos;
+  std::vector<double> polar_weight;
+  if(num_polar > 0){
+    // make sure the user specified the correct number of quadrature points
+    if(num_polar % 4 != 0){
+      printf("Must have number of polar angles be a multiple of 4\n");
+      MPI_Abort(MPI_COMM_WORLD, 1);
+    }
+    if(num_azimuth % 2 != 0){
+      printf("Must have number of azimuthal angles be a multiple of 2\n");
+      MPI_Abort(MPI_COMM_WORLD, 1);
+    }
+    if(num_polar*num_azimuth != num_directions){
+      printf("You need to specify %d total directions, not %d\n",
+          num_polar*num_azimuth, num_directions);
+      MPI_Abort(MPI_COMM_WORLD, 1);
+    }
+
+    // Compute gauss legendre weights
+    polar_cos.resize(num_polar);
+    polar_weight.resize(num_polar);
+    GaussLegendre(-1.0, 1.0, polar_cos, polar_weight, DBL_EPSILON);
+
+    // compute azmuhtal angles and weights
+    std::vector<double> az_angle(num_azimuth);
+    std::vector<double> az_weight(num_azimuth);
+    double dangle = 2.0*M_PI/((double) num_azimuth);
+
+    for(int i=0; i<num_azimuth; i++){
+      if(i == 0){
+        az_angle[0] = dangle/2.0;
+      }
+      else{
+        az_angle[i] = az_angle[i-1] + dangle;
+      }
+      az_weight[i] = dangle;
+    }
+
+
+    // Loop over polar 'octants
+    int d = 0;
+    for(int i=0; i< num_polar; i++){
+      for(int j=0; j< num_azimuth; j++){
+        double xcos = sqrt(1.0-polar_cos[i]*polar_cos[i]) * cos(az_angle[j]);
+        double ycos = sqrt(1.0-polar_cos[i]*polar_cos[i]) * sin(az_angle[j]);
+        double zcos = polar_cos[i];
+        double w = polar_weight[i]*az_weight[j];
+
+        directions[d].id = (xcos > 0.) ? 1 : -1;
+        directions[d].jd = (ycos > 0.) ? 1 : -1;
+        directions[d].kd = (zcos > 0.) ? 1 : -1;
+
+        directions[d].octant = 0;
+        if(directions[d].id == -1){
+          directions[d].octant += 1;
+        }
+        if(directions[d].jd == -1){
+          directions[d].octant += 2;
+        }
+        if(directions[d].kd == -1){
+          directions[d].octant += 4;
+        }
+
+        directions[d].xcos = std::abs(xcos);
+        directions[d].ycos = std::abs(ycos);
+        directions[d].zcos = std::abs(zcos);
+        directions[d].w = w;
+
+        ++ d;
+      }
+    }
+
+    // Sort by octant.. so each set has same directions
+    std::sort(directions.begin(), directions.end(), dirSortFcn);
+  }
+  else{
+    // Do (essentialy) an S2 quadrature.. but with repeated directions
+
+    // Compute x,y,z cosine values
+    double mu  = cos(M_PI/4);
+    double eta = sqrt(1-mu*mu) * cos(M_PI/4);
+    double xi  = sqrt(1-mu*mu) * sin(M_PI/4);
+    int d = 0;
+    for(int octant = 0;octant < 8;++ octant){
+      double omegas[3];
+      omegas[0] = octant & 0x1;
+      omegas[1] = (octant>>1) & 0x1;
+      omegas[2] = (octant>>2) & 0x1;
+
+      for(int sd=0; sd<num_directions_per_octant; sd++, d++){
+        // Store which logical direction of travel we have
+        directions[d].id = (omegas[0] > 0.) ? 1 : -1;
+        directions[d].jd = (omegas[1] > 0.) ? 1 : -1;
+        directions[d].kd = (omegas[2] > 0.) ? 1 : -1;
+
+        // Store quadrature point's weight
+        directions[d].w = 4.0*M_PI / (double)num_directions;
+        directions[d].xcos = mu;
+        directions[d].ycos = eta;
+        directions[d].zcos = xi;
+      }
+    }
+  }
+}
+
+
+
+

miniapps/kripke/Kripke/Directions.h

+/*--------------------------------------------------------------------------
+ * Header file for the Directions data structures
+ *--------------------------------------------------------------------------*/
+
+#ifndef KRIPKE_DIRECTIONS_H__
+#define KRIPKE_DIRECTIONS_H__
+
+#include <vector>
+
+class Grid_Data;
+struct Input_Variables;
+
+/**
+ * Contains information needed for one quadrature set direction.
+ */
+struct Directions{
+  double xcos;              /* Absolute value of the x-direction cosine. */
+  double ycos;              /* Absolute value of the y-direction cosine. */
+  double zcos;              /* Absolute value of the z-direction cosine. */
+  double w;                 /* weight for the quadrature rule.*/
+  int id;                   /* direction flag (= 1 if x-direction
+                            cosine is positive; = -1 if not). */
+  int jd;                   /* direction flag (= 1 if y-direction
+                            cosine is positive; = -1 if not). */
+  int kd;                   /* direction flag (= 1 if z-direction
+                            cosine is positive; = -1 if not). */
+  int octant;
+};
+
+
+void InitDirections(Grid_Data *grid_data, Input_Variables *input_vars);
+
+#endif

miniapps/kripke/Kripke/Grid.cpp

+#include <Kripke/Grid.h>
+
+#include <Kripke/Input_Variables.h>
+#include <Kripke/Layout.h>
+#include <Kripke/SubTVec.h>
+#include <cmath>
+#include <sstream>
+#include <mpi.h>
+
+#ifdef KRIPKE_USE_SILO
+#include <sys/stat.h>
+#include <silo.h>
+#include <string.h>
+#endif
+
+/**
+ * Grid_Data constructor
+*/
+Grid_Data::Grid_Data(Input_Variables *input_vars)
+{
+  // Create object to describe processor and subdomain layout in space
+  // and their adjacencies
+  Layout *layout = createLayout(input_vars);
+
+  // create the kernel object based on nesting
+  kernel = createKernel(input_vars->nesting, 3);
+
+  // Create quadrature set (for all directions)
+  int total_num_directions = input_vars->num_dirsets_per_octant * input_vars->num_dirs_per_dirset;
+  InitDirections(this, input_vars);
+
+  num_direction_sets = 8*input_vars->num_dirsets_per_octant;
+  num_directions_per_set = input_vars->num_dirs_per_dirset;
+  num_group_sets = input_vars->num_groupsets;
+  num_groups_per_set = input_vars->num_groups_per_groupset;
+  num_zone_sets = 1;
+  for(int dim = 0;dim < 3;++ dim){
+    num_zone_sets *= input_vars->num_zonesets_dim[dim];
+  }
+
+  legendre_order = input_vars->legendre_order;
+  total_num_moments = (legendre_order+1)*(legendre_order+1);
+
+  int num_subdomains = num_direction_sets*num_group_sets*num_zone_sets;
+
+  Nesting_Order nest = input_vars->nesting;
+
+  /* Set ncalls */
+  niter = input_vars->niter;
+
+  // setup mapping of moments to legendre coefficients
+  moment_to_coeff.resize(total_num_moments);
+  int nm = 0;
+  for(int n = 0;n < legendre_order+1;++ n){
+    for(int m = -n;m <= n; ++ m){
+      moment_to_coeff[nm] = n;
+      ++ nm;
+    }
+  }
+
+  // setup cross-sections
+  int total_num_groups = num_group_sets*num_groups_per_set;
+  sigma_tot.resize(total_num_groups, 0.0);
+
+  // Setup scattering transfer matrix for 3 materials
+  sigs.resize(3);
+  for(int mat = 0;mat < 3;++ mat){
+    // allocate transfer matrix
+    sigs[mat] = new SubTVec(kernel->nestingSigs(), total_num_groups, legendre_order+1, total_num_groups);
+
+    // Set to isotropic scattering given user inputs
+    sigs[mat]->clear(0.0);
+    for(int g = 0;g < total_num_groups;++ g){
+      (*sigs[mat])(g, 0, g) = input_vars->sigs[mat];
+    }
+  }
+
+  // just allocate pointer vectors, we will allocate them below
+  ell.resize(num_direction_sets, NULL);