Commit ef0aba59 authored by Nicole Marsaglia's avatar Nicole Marsaglia

Kripke miniapp that utilizes the conduit data adaptor.

parent 571a82ed
......@@ -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}")
......@@ -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()
###############################################################################
# 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()
/**
* 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
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)
#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;
}
}
}
}
/*--------------------------------------------------------------------------
* 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
This diff is collapsed.
#ifndef KRIPKE_GRID_DATA_H__
#define KRIPKE_GRID_DATA_H__
#include <Kripke/Directions.h>
#include <Kripke/Kernel.h>
#include <Kripke/Subdomain.h>
#include <Kripke/Timing.h>
#include <mpi.h>
#include <vector>
// Foreward Decl
struct Input_Variables;
struct Grid_Data;
struct SubTVec;
/**
* Contains all grid parameters and variables.
*/
struct Grid_Data {
public:
explicit Grid_Data(Input_Variables *input_vars);
~Grid_Data();
void randomizeData(void);
void copy(Grid_Data const &b);
bool compare(Grid_Data const &b, double tol, bool verbose);
double particleEdit(void);
#ifdef KRIPKE_USE_SILO
void writeSilo(std::string const &fname);
#endif
Timing timing;
int niter;
double source_value;
std::vector<double> sigma_tot; // Cross section data
int num_group_sets; // Number of group-sets
int num_groups_per_set; // How many groups in each set
int num_direction_sets; // Number of direction-sets
int num_directions_per_set; // Number of directions per dir set
int num_zone_sets; // Number of zone sets
int legendre_order; // Legendra expansion order ( >= 0 )
int total_num_moments; // Number of spherical harmonic moments
std::vector<int> moment_to_coeff; // Map from harmonic moments to legendre coefficients
std::vector<Directions> directions; // Quadrature point data, for all directions
Kernel *kernel; // Layout-specific math kernels
std::vector<Subdomain> subdomains; // Group/Angle/Zone set data
std::vector<int> zs_to_sdomid; // map of zonesets to subdomains with ds=gs=0
// Variables:
std::vector<SubTVec *> sigs; // scattering lookup table for each material
// Per directionset ell and ell_plus matrices (Subdomain point into these arrays)
std::vector<SubTVec *> ell; // L matrix in nm_offset coordinates
std::vector<SubTVec *> ell_plus; // L+ matrix in nm_offset coordinates
// Per zoneset phi and phi_out (Subdomains point into these arrays)
std::vector<SubTVec *> phi; // Moments of psi
std::vector<SubTVec *> phi_out; // Scattering source
};
#endif
/*--------------------------------------------------------------------------
* Header file for the Input_Variables structure
*--------------------------------------------------------------------------*/
#ifndef KRIPKE_INPUT_VARIABLES_H__
#define KRIPKE_INPUT_VARIABLES_H__
#include<Kripke.h>
/**
* This structure defines the input parameters to setup a problem.
*/
struct Input_Variables {
int npx, npy, npz; // The number of processors in x,y,z
int nx, ny, nz; // Number of spatial zones in x,y,z
int num_dirsets_per_octant;
int num_dirs_per_dirset;
int num_groupsets;
int num_groups_per_groupset; //
int num_zonesets_dim[3]; // number of zoneset in x, y, z
int niter; // number of solver iterations to run
int legendre_order; // Scattering order (number Legendre coeff's - 1)
int layout_pattern; // Which subdomain/task layout to use
int quad_num_polar; // Number of polar quadrature points
int quad_num_azimuthal; // Number of azimuthal quadrature points
ParallelMethod parallel_method;
double sigt[3]; // total cross section for 3 materials
double sigs[3]; // total scattering cross section for 3 materials
#ifdef KRIPKE_USE_SILO
std::string silo_basename; // name prefix for silo output files
#endif
Nesting_Order nesting; // Data layout and loop ordering (of Psi)
};
#endif
#include<Kripke/Kernel.h>
#include<Kripke/Grid.h>
#include<Kripke/SubTVec.h>
#include<Kripke/Kernel/Kernel_3d_GDZ.h>
#include<Kripke/Kernel/Kernel_3d_DGZ.h>
#include<Kripke/Kernel/Kernel_3d_ZDG.h>
#include<Kripke/Kernel/Kernel_3d_DZG.h>
#include<Kripke/Kernel/Kernel_3d_ZGD.h>
#include<Kripke/Kernel/Kernel_3d_GZD.h>
/**
* Factory to create a kernel object for the specified nesting
*/
Kernel *createKernel(Nesting_Order nest, int num_dims){
if(num_dims == 3){
switch(nest){
case NEST_GDZ:
return new Kernel_3d_GDZ();
case NEST_DGZ:
return new Kernel_3d_DGZ();
case NEST_ZDG:
return new Kernel_3d_ZDG();
case NEST_DZG:
return new Kernel_3d_DZG();
case NEST_ZGD:
return new Kernel_3d_ZGD();
case NEST_GZD:
return new Kernel_3d_GZD();
}
}
MPI_Abort(MPI_COMM_WORLD, 1);
}
/*--------------------------------------------------------------------------
* Header file for the Grid_Data data structures
*--------------------------------------------------------------------------*/
#ifndef KRIPKE_KERNEL_H__
#define KRIPKE_KERNEL_H__
#include <Kripke.h>
struct Grid_Data;
struct SubTVec;
struct Subdomain;
/**
* This is the Kernel base-class and interface definition.
* This abstracts the storage of Psi, Phi, L, L+ from the rest of the code,
* providing data-layout specific routines.
*/
class Kernel {
public:
virtual Nesting_Order nestingPsi(void) const = 0;
virtual Nesting_Order nestingPhi(void) const = 0;
virtual Nesting_Order nestingSigt(void) const = 0;
virtual Nesting_Order nestingEll(void) const = 0;
virtual Nesting_Order nestingEllPlus(void) const = 0;
virtual Nesting_Order nestingSigs(void) const = 0;
// Computational Kernels
virtual void LTimes(Grid_Data *grid_data) = 0;
virtual void LPlusTimes(Grid_Data *grid_data) = 0;
virtual void scattering(Grid_Data *grid_data) = 0;
virtual void source(Grid_Data *grid_data) = 0;
virtual void sweep(Subdomain *ga_set) = 0;
};
// Factory to create correct kernel object
Kernel *createKernel(Nesting_Order, int num_dims);
#endif
This diff is collapsed.
#ifndef KRIPKE_KERNEL_3D_DGZ_H__
#define KRIPKE_KERNEL_3D_DGZ_H__
#include<Kripke/Kernel.h>
class Kernel_3d_DGZ : public Kernel {
public: