diff --git a/INSTALL b/INSTALL index be19f667..962ba9f2 100644 --- a/INSTALL +++ b/INSTALL @@ -37,9 +37,12 @@ Step 0: Build any necessary third-party library as needed -AMPE requires the third-party packages hypre and SAMRAI. These libaries -are typically not pre-installed on most systems and need to be installed. +AMPE requires the third-party packages hypre, Sundials and SAMRAI. These +libaries are typically not pre-installed on most systems and need to be +installed. + SAMRAI: compatible with version 4.0.0 or later, built without Sundials. +Sundials: compatible with version 5.4.0 or later. To build SAMRAI4.0.0, do something similar to this: -------------------------------------------------- @@ -64,79 +67,21 @@ make install //////////////////////////////////////////////////////////////////////////// -Step 1: -Run autoconf: - - autoconf - -This should generate the 'configure' file - -//////////////////////////////////////////////////////////////////////////// - -Step 2: Configure for a specific combination of platform, compiler and - optimization level - -After checking out from the repository, cd into the top level -directory and execute the configure script for the -platform-compiler-optimization you want to build. For example, to -configure for syrah (LLNL) using the Intel C++ and Fortran compilers with -optimization, execute - - scripts/doconfig_peloton_icpc_ifort_opt - -If a script does not exist for the platform-compiler_optimization -combination you want, you'll have to create your own using one of the -existing scripts as an example. Note that all of the scripts -currently in the "scripts" subdirectory assume that they are being -executed from the top level directory (as indicated above), which is -the location of of the configure script the "doconfig" scripts depend -upon. - -//////////////////////////////////////////////////////////////////////////// - -Step 3: Install the base (i.e., third-party) libraries for the - new platform-compiler-optimization combination - -The configuration step just performed will create a subdirectory -"build" of the top level directory (if it didn't already exist) and a -subdirectory of the "build" directory labeled by the associated -platform-compiler-optimization combination. Underneath that -subdirectory are three more subdirectories: base, config, and objects. -In other words, for the current example, the following directory tree -is created if ruuning the script on syrah: - -(top level directory) - | - | - build - | - | - syrah-mpic++-mpif77-opt - | | | - | | | - base config objects - -To install the base libraries, cd into the base subdirectory and -execute the doinstall script located there: - - cd build/syrah-mpic++-mpif77-opt/base - ./doinstall - -This will compile and install the base libraries, which currently -consist of Sundials (including CPODES). - -//////////////////////////////////////////////////////////////////////////// +Step 1: Build the AMPE code with cmake -Step 4: Build the AMPE code with cmake +mkdir build +cd build From the build directory, run a "cmake script". For example on syrah: - ../scripts/peloton_cmake +../scripts/peloton_cmake where the solver dimension is specified (2d or 3d). If a script does not exist for the platform-compiler_optimization combination you want, you'll have to create your own using one of the existing "cmake" scripts as an example. +Note that the option "-DUSE_CVODE=ON" is required for that build to work. + Then run: make @@ -145,19 +90,3 @@ as executable unit tests. To run the test suite, run: make test - -//////////////////////////////////////////////////////////////////////////// - -Miscellaneous notes: - -1. The configuration performed in Step 1 generates the Makefile in -build/syrah-mpic++-mpif77-opt/objects from the template Makefile.in in -source. This includes the generation of the file Makefile.depend in -source, which describes the build dependencies. This means that if -any changes are made to the source code that could affect the -dependencies, the configuration step should be redone. - -2. For the "CVODE" build (no projections, not quaternions propagation), -one should link with an external (Sundials v3.0.0 or later) using cmake -options like: --DUSE_CVODE=ON -DSUNDIALS_DIR=$HOME/sundials/sundials-3.0.0/install diff --git a/README.md b/README.md index c95e862a..1c0474fb 100644 --- a/README.md +++ b/README.md @@ -17,8 +17,6 @@ Dependencies * [SAMRAI] (https://github.com/LLNL/SAMRAI) -* [CPODES] (https://simtk.org/projects/cpodes) - * [Sundials] (https://github.com/LLNL/sundials) * [HDF5] (https://support.hdfgroup.org/HDF5) @@ -27,12 +25,6 @@ Dependencies * [BOOST] (http://www.boost.org) -Since CPODES is currently not distributed with Sundials, and SAMRAI -does not supports the lastest SUNDIALS release, modifications to -these libraries had to be made. The modified libraries are distributed -with this code under the 'base' directory and need to be built before -building the main code. - References ---------- diff --git a/scripts/condo_cmake_cvode_2d b/scripts/condo_cmake_cvode_2d index 99790dd5..e642c257 100755 --- a/scripts/condo_cmake_cvode_2d +++ b/scripts/condo_cmake_cvode_2d @@ -11,7 +11,7 @@ cmake -DCMAKE_CXX_COMPILER=mpiCC -DCMAKE_C_COMPILER=mpicc \ -DCMAKE_Fortran_COMPILER=mpif77 \ -DSAMRAI_DIR=$HOME/SAMRAI/SAMRAI-v4.0.0 \ -DHYPRE_DIR=$HOME/hypre/hypre-2.11.1/src/hypre \ - -DSUNDIALS_DIR=$HOME/sundials/sundials-5.3.0/install \ + -DSUNDIALS_DIR=$HOME/sundials/sundials-5.4.0/install \ -DNetCDF_CXX_ROOT=$NETCDF_DIR \ -DCMAKE_BUILD_TYPE=Release -DBUILD_SHARED_LIBS=OFF \ -DAMPE_WITH_CLANG_FORMAT=ON \ diff --git a/scripts/condo_cmake_cvode_3d b/scripts/condo_cmake_cvode_3d new file mode 100755 index 00000000..ca4dce59 --- /dev/null +++ b/scripts/condo_cmake_cvode_3d @@ -0,0 +1,22 @@ +#!/bin/sh +source ../scripts/modules.condo + +rm -rf CMakeCache.txt +rm -rf CMakeFiles/ +rm cmake_install.cmake +rm Makefile +rm -f ../source/fortran/2d/*.f + +cmake -DCMAKE_CXX_COMPILER=mpiCC -DCMAKE_C_COMPILER=mpicc \ + -DCMAKE_Fortran_COMPILER=mpif77 \ + -DSAMRAI_DIR=$HOME/SAMRAI/SAMRAI-v4.0.0 \ + -DHYPRE_DIR=$HOME/hypre/hypre-2.11.1/src/hypre \ + -DSUNDIALS_DIR=$HOME/sundials/sundials-5.4.0/install \ + -DNetCDF_CXX_ROOT=$NETCDF_DIR \ + -DCMAKE_BUILD_TYPE=Release -DBUILD_SHARED_LIBS=OFF \ + -DAMPE_WITH_CLANG_FORMAT=ON \ + -DCMAKE_PREFIX_PATH=${HOME}/bin \ + -DNDIM="3" \ + -DUSE_CVODE=ON \ + .. + diff --git a/source/CVODEAbstractFunctions.h b/source/CVODEAbstractFunctions.h index 480a94f8..e1391744 100644 --- a/source/CVODEAbstractFunctions.h +++ b/source/CVODEAbstractFunctions.h @@ -59,8 +59,7 @@ class CVODEAbstractFunctions * IMPORTANT: This function must not modify the vector y. */ virtual int evaluateRHSFunction(double t, SundialsAbstractVector* y, - SundialsAbstractVector* y_dot, - int fd_flag) = 0; + SundialsAbstractVector* y_dot) = 0; /** * User-supplied function for setting up the preconditioner @@ -81,6 +80,12 @@ class CVODEAbstractFunctions SundialsAbstractVector* r, SundialsAbstractVector* z, double gamma, double delta, int lr) = 0; + virtual int applyProjection(double t, SundialsAbstractVector* y, + SundialsAbstractVector* corr, double epsProj, + SundialsAbstractVector* err) = 0; + + virtual int evaluateJTimesRHSFunction(double t, SundialsAbstractVector* y, + SundialsAbstractVector* y_dot) = 0; }; #endif diff --git a/source/CVODESolver.C b/source/CVODESolver.C index d59569a6..1eb08ee1 100644 --- a/source/CVODESolver.C +++ b/source/CVODESolver.C @@ -84,6 +84,8 @@ CVODESolver::CVODESolver(const std::string& object_name, setCVSpgmrToleranceScaleFactor(0); d_CVODE_needs_initialization = true; + d_uses_projectionfn = false; + d_uses_jtimesrhsfn = false; } CVODESolver::~CVODESolver() @@ -184,6 +186,8 @@ void CVODESolver::initializeCVODE() ierr = CVSpilsSetLinearSolver(d_cvode_mem, d_linear_solver); CVODE_SAMRAI_ERROR(ierr); + CVodeSetLSetupFrequency(d_cvode_mem, d_max_precond_steps); + if (!(d_max_order < 1)) { ierr = CVodeSetMaxOrd(d_cvode_mem, d_max_order); CVODE_SAMRAI_ERROR(ierr); @@ -226,8 +230,22 @@ void CVODESolver::initializeCVODE() CVODE_SAMRAI_ERROR(ierr); } + if (d_uses_projectionfn) { + CVProjFn proj_fn = CVODESolver::CVODEProjEval; + ierr = CVodeSetProjFn(d_cvode_mem, proj_fn); + } + + if (d_uses_jtimesrhsfn) { + CVRhsFn jtimesrhs_fn = CVODESolver::CVODEJTimesRHSFuncEval; + ierr = CVodeSetJacTimesRhsFn(d_cvode_mem, jtimesrhs_fn); + } + + ierr = CVodeSetLSNormFactor(d_cvode_mem, -1.0); + CVODE_SAMRAI_ERROR(ierr); + } // if no need to initialize CVODE, function does nothing + d_CVODE_needs_initialization = false; } diff --git a/source/CVODESolver.h b/source/CVODESolver.h index c436aeca..ff685b5b 100644 --- a/source/CVODESolver.h +++ b/source/CVODESolver.h @@ -649,6 +649,16 @@ class CVODESolver d_CVODE_needs_initialization = true; } + void setProjectionFunction(bool uses_projectionfn) + { + d_uses_projectionfn = uses_projectionfn; + } + + void setJTimesRhsFunction(bool uses_jtimesrhsfn) + { + d_uses_jtimesrhsfn = uses_jtimesrhsfn; + } + /** * Get solution vector. */ @@ -792,6 +802,8 @@ class CVODESolver */ void printCVODEStatistics(std::ostream& os) const; + void setMaxPrecondSteps(int max_steps) { d_max_precond_steps = max_steps; } + // CVODE optional return values. /** @@ -1127,14 +1139,9 @@ class CVODESolver static int CVODERHSFuncEval(realtype t, N_Vector y, N_Vector y_dot, void* my_solver) { - // evaluateRHSFunction requires an argument "fd_flag" - // we set it to 0 always since this interface does not let - // us provide one - int fd_flag = 0; return ((CVODESolver*)my_solver) ->getCVODEFunctions() - ->evaluateRHSFunction(t, SABSVEC_CAST(y), SABSVEC_CAST(y_dot), - fd_flag); + ->evaluateRHSFunction(t, SABSVEC_CAST(y), SABSVEC_CAST(y_dot)); } /* @@ -1165,6 +1172,27 @@ class CVODESolver return success; } + static int CVODEProjEval(realtype t, N_Vector y, N_Vector corr, + realtype epsProj, N_Vector err, void* my_solver) + { + int success = + ((CVODESolver*)my_solver) + ->getCVODEFunctions() + ->applyProjection(t, SABSVEC_CAST(y), SABSVEC_CAST(corr), epsProj, + SABSVEC_CAST(err)); + return success; + } + + static int CVODEJTimesRHSFuncEval(realtype t, N_Vector y, N_Vector y_dot, + void* my_solver) + { + int success = ((CVODESolver*)my_solver) + ->getCVODEFunctions() + ->evaluateJTimesRHSFunction(t, SABSVEC_CAST(y), + SABSVEC_CAST(y_dot)); + return success; + } + /* * Open CVODE log file, allocate main memory for CVODE and initialize * CVODE memory record. CVODE is initialized based on current state @@ -1264,6 +1292,25 @@ class CVODESolver * CVODEAbstractFunctions. */ bool d_uses_preconditioner; + + /* + * Boolean flag indicating whether a user-supplied projection + * routine is provided in the concrete subclass of + * CVODEAbstractFunctions. + */ + bool d_uses_projectionfn; + + /* + * Boolean flag indicating whether a different RHS + * routine for Jacobian -vector products is provided + * in the concrete subclass of CVODEAbstractFunctions. + */ + bool d_uses_jtimesrhsfn; + + /* + * Maximum number of steps between Jacobian evaluations + */ + int d_max_precond_steps; }; #endif diff --git a/source/QuatIntegrator.C b/source/QuatIntegrator.C index 282d7d92..f24594b4 100644 --- a/source/QuatIntegrator.C +++ b/source/QuatIntegrator.C @@ -71,14 +71,10 @@ #ifdef USE_CPODE #define BDF CP_BDF #define ONE_STEP CP_ONE_STEP -#define NEWTON CP_NEWTON -#define SUNDIALS_PREFIX CP #include "cpodes/cpodes_spils.h" #else #define BDF CV_BDF #define ONE_STEP CV_ONE_STEP -#define NEWTON CV_NEWTON -#define SUNDIALS_PREFIX CV #include "cvode/cvode_spils.h" #endif @@ -87,9 +83,6 @@ #ifndef TRUE #define TRUE (1) #endif -#ifndef FALSE -#define FALSE (0) -#endif //----------------------------------------------------------------------- @@ -1624,13 +1617,13 @@ void QuatIntegrator::setSundialsOptions() #ifdef USE_CPODE // if ( d_show_solver_stats ) d_sundials_solver->printDiagnostics( true ); d_sundials_solver->setMaxKrylovDimension(d_max_krylov_dimension); - d_sundials_solver->setMaxPrecondSteps(d_max_precond_steps); // d_sundials_solver->setCPSpgmrToleranceScaleFactor(0.005); // d_sundials_solver->setCPNonlinConvCoef(0.1); - d_sundials_solver->setIterationType(NEWTON); + d_sundials_solver->setIterationType(CP_NEWTON); #else // d_sundials_solver->setCVSpgmrToleranceScaleFactor(0.005); #endif + d_sundials_solver->setMaxPrecondSteps(d_max_precond_steps); d_sundials_solver->setLinearMultistepMethod(BDF); d_sundials_solver->setSteppingMethod(ONE_STEP); d_sundials_solver->setInitialStepSize(d_previous_timestep); @@ -1649,6 +1642,10 @@ void QuatIntegrator::setSundialsOptions() bool needs_initialization = true; d_sundials_solver->setFinalValueOfIndependentVariable(d_end_time, needs_initialization); +#ifndef USE_CPODE + d_sundials_solver->setProjectionFunction(true); + d_sundials_solver->setJTimesRhsFunction(true); +#endif } //----------------------------------------------------------------------- @@ -4526,8 +4523,6 @@ int QuatIntegrator::applyConcentrationPreconditioner( return retcode; } -#ifdef USE_CPODE - //----------------------------------------------------------------------- // Virtual function from CPODESAbstractFunction @@ -4563,8 +4558,6 @@ int QuatIntegrator::applyProjection(double time, SundialsAbstractVector* y, return 0; // Always successful } -#endif - //======================================================================= void QuatIntegrator::correctRhsForSymmetry( diff --git a/source/QuatIntegrator.h b/source/QuatIntegrator.h index 1db7fae1..2e215d03 100644 --- a/source/QuatIntegrator.h +++ b/source/QuatIntegrator.h @@ -259,6 +259,19 @@ class QuatIntegrator : public mesh::StandardTagAndInitStrategy // // Methods inherited from CVODEAbstractFunctions // + int applyProjection(double time, SundialsAbstractVector* y, + SundialsAbstractVector* corr, double epsProj, + SundialsAbstractVector* err); + int evaluateJTimesRHSFunction(double t, SundialsAbstractVector* y, + SundialsAbstractVector* y_dot) + { + evaluateRHSFunction(t, y, y_dot, 1); + } + int evaluateRHSFunction(double t, SundialsAbstractVector* y, + SundialsAbstractVector* y_dot) + { + evaluateRHSFunction(t, y, y_dot, 0); + } int CVSpgmrPrecondSet(double t, SundialsAbstractVector* y, SundialsAbstractVector* fy, int jok, int* jcurPtr, diff --git a/source/Sundials_SAMRAIVector.C b/source/Sundials_SAMRAIVector.C index a36216e7..0d3820a4 100644 --- a/source/Sundials_SAMRAIVector.C +++ b/source/Sundials_SAMRAIVector.C @@ -199,6 +199,8 @@ int Sundials_SAMRAIVector::testReciprocal(const SundialsAbstractVector* x) } #ifndef USE_CPODE +// note: this function has been modified to include the missing +// "allreduce" in the SAMRAI getLenght function sunindextype Sundials_SAMRAIVector::getLength() const { int64_t len = d_samrai_vector->getLength(); diff --git a/tests/AlCu/test3d.py b/tests/AlCu/test3d.py index b13cbc1c..c1695c06 100644 --- a/tests/AlCu/test3d.py +++ b/tests/AlCu/test3d.py @@ -25,7 +25,7 @@ lines=output.split(b'\n') end_reached = False -end_time = 0.00049 +end_time = 0.000499 first_concentration=-1. target_solid_fraction = 0.0614 for line in lines: @@ -57,7 +57,7 @@ if end_reached: words=line.split() sfraction=eval(words[6]) - if abs(sfraction-target_solid_fraction)>1.e-4: + if abs(sfraction-target_solid_fraction)>2.e-4: print("Wrong solid fraction") print("{} instead of {}".format(sfraction,target_solid_fraction)) sys.exit(1) diff --git a/tests/GaussianT/test3d.py b/tests/GaussianT/test3d.py index a20f1995..21e9d3e5 100644 --- a/tests/GaussianT/test3d.py +++ b/tests/GaussianT/test3d.py @@ -65,7 +65,7 @@ words=line.split() maxT = eval(words[3]) -tol = 1.e-1 +tol = 2.e-1 expected_max_T = 1423. + final_time * 500. if abs(maxT-expected_max_T)>tol: print("max. T = {}".format(maxT)) diff --git a/tests/TwoBilayers/test2d.py b/tests/TwoBilayers/test2d.py index b1447f33..4831a002 100644 --- a/tests/TwoBilayers/test2d.py +++ b/tests/TwoBilayers/test2d.py @@ -58,7 +58,7 @@ words=line.split() oe = eval(words[4]) print("Orient enegy = {}".format(oe)) - if abs(oe-0.0049)>1.e-4: + if abs(oe-0.0049)>3.e-4: print("Wrong orient energy") sys.exit(1) diff --git a/tests/TwoBilayers/test3d.py b/tests/TwoBilayers/test3d.py index 830912f8..da521254 100644 --- a/tests/TwoBilayers/test3d.py +++ b/tests/TwoBilayers/test3d.py @@ -58,7 +58,7 @@ words=line.split() oe = eval(words[4]) print("Orient enegy = {}".format(oe)) - if abs(oe-0.00025)>1.e-5: + if abs(oe-0.00026)>1.e-5: print("Wrong orient energy") sys.exit(1)