diff --git a/CMakeLists.txt b/CMakeLists.txt index 6bd5cc39fe..24f4baa097 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -30,13 +30,13 @@ project(sundials C) # Set some variables with info on the SUNDIALS project set(PACKAGE_BUGREPORT "woodward6@llnl.gov") set(PACKAGE_NAME "SUNDIALS") -set(PACKAGE_STRING "SUNDIALS 5.2.0") +set(PACKAGE_STRING "SUNDIALS 5.3.0") set(PACKAGE_TARNAME "sundials") # set SUNDIALS version numbers # (use "" for the version label if none is needed) set(PACKAGE_VERSION_MAJOR "5") -set(PACKAGE_VERSION_MINOR "2") +set(PACKAGE_VERSION_MINOR "3") set(PACKAGE_VERSION_PATCH "0") set(PACKAGE_VERSION_LABEL "") @@ -63,37 +63,37 @@ mark_as_advanced(CLEAR # Specify the VERSION and SOVERSION for shared libraries -set(arkodelib_VERSION "4.2.0") +set(arkodelib_VERSION "4.3.0") set(arkodelib_SOVERSION "4") -set(cvodelib_VERSION "5.2.0") +set(cvodelib_VERSION "5.3.0") set(cvodelib_SOVERSION "5") -set(cvodeslib_VERSION "5.2.0") +set(cvodeslib_VERSION "5.3.0") set(cvodeslib_SOVERSION "5") -set(idalib_VERSION "5.2.0") +set(idalib_VERSION "5.3.0") set(idalib_SOVERSION "5") -set(idaslib_VERSION "4.2.0") +set(idaslib_VERSION "4.3.0") set(idaslib_SOVERSION "4") -set(kinsollib_VERSION "5.2.0") +set(kinsollib_VERSION "5.3.0") set(kinsollib_SOVERSION "5") set(cpodeslib_VERSION "0.0.0") set(cpodeslib_SOVERSION "0") -set(nveclib_VERSION "5.2.0") +set(nveclib_VERSION "5.3.0") set(nveclib_SOVERSION "5") -set(sunmatrixlib_VERSION "3.2.0") +set(sunmatrixlib_VERSION "3.3.0") set(sunmatrixlib_SOVERSION "3") -set(sunlinsollib_VERSION "3.2.0") +set(sunlinsollib_VERSION "3.3.0") set(sunlinsollib_SOVERSION "3") -set(sunnonlinsollib_VERSION "2.2.0") +set(sunnonlinsollib_VERSION "2.3.0") set(sunnonlinsollib_SOVERSION "2") # Specify the location of additional CMAKE modules @@ -232,6 +232,13 @@ show_variable(SUNDIALS_INDEX_TYPE STRING "${DOCSTR}" "") mark_as_advanced(SUNDIALS_INDEX_TYPE) include(SundialsIndexSize) +# --------------------------------------------------------------- +# Option to specify monitoring +# --------------------------------------------------------------- + +set(DOCSTR "Build with simulation monitoring capabilities enabled") +sundials_option(SUNDIALS_BUILD_WITH_MONITORING BOOL ${DOCSTR} OFF) + # --------------------------------------------------------------- # Enable Fortran interface? # --------------------------------------------------------------- @@ -858,6 +865,17 @@ endif(CUDA_ENABLE) # Now that all languages are setup, we can configure them more. # --------------------------------------------------------------- +# --------------------------------------------------------------- +# Option to use specialized fused kernels in the packages. +# Currently only available in CVODE. +# --------------------------------------------------------------- + +if(CUDA_ENABLE AND CMAKE_CUDA_COMPILER AND BUILD_CVODE) + set(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS FALSE CACHE BOOL "Build specialized fused CUDA kernels") +else() + set(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS FALSE CACHE BOOL "Build specialized fused CUDA kernels" FORCE) +endif() + # --------------------------------------------------------------- # Decide how to compile MPI codes. We must check for MPI if # MPI is enabled or if Trilinos is enabled because the Trilinos diff --git a/INSTALL_GUIDE.pdf b/INSTALL_GUIDE.pdf index 82f9cd87e2..dfcb21916e 100644 Binary files a/INSTALL_GUIDE.pdf and b/INSTALL_GUIDE.pdf differ diff --git a/README.md b/README.md index 781dc4d75a..18aa60f264 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,5 @@ # SUNDIALS: SUite of Nonlinear and DIfferential/ALgebraic equation Solvers # -### Version 5.2.0 (Mar 2020) ### +### Version 5.3.0 (May 2020) ### **Center for Applied Scientific Computing, Lawrence Livermore National Laboratory** diff --git a/doc/arkode/ark_examples.pdf b/doc/arkode/ark_examples.pdf index 35351fd2ef..1303ff3d0d 100644 Binary files a/doc/arkode/ark_examples.pdf and b/doc/arkode/ark_examples.pdf differ diff --git a/doc/arkode/ark_guide.pdf b/doc/arkode/ark_guide.pdf index 89382838d4..d3d78e4b5e 100644 Binary files a/doc/arkode/ark_guide.pdf and b/doc/arkode/ark_guide.pdf differ diff --git a/doc/cvode/cv_examples.pdf b/doc/cvode/cv_examples.pdf index 23c84191dd..2e25f21e0a 100644 Binary files a/doc/cvode/cv_examples.pdf and b/doc/cvode/cv_examples.pdf differ diff --git a/doc/cvode/cv_guide.pdf b/doc/cvode/cv_guide.pdf index ea9d2c5cec..8c7eecac22 100644 Binary files a/doc/cvode/cv_guide.pdf and b/doc/cvode/cv_guide.pdf differ diff --git a/doc/cvodes/cvs_examples.pdf b/doc/cvodes/cvs_examples.pdf index daa44e4bf2..0f7ceec7a8 100644 Binary files a/doc/cvodes/cvs_examples.pdf and b/doc/cvodes/cvs_examples.pdf differ diff --git a/doc/cvodes/cvs_guide.pdf b/doc/cvodes/cvs_guide.pdf index 0d8bb9680a..793e1fcb1d 100644 Binary files a/doc/cvodes/cvs_guide.pdf and b/doc/cvodes/cvs_guide.pdf differ diff --git a/doc/ida/ida_examples.pdf b/doc/ida/ida_examples.pdf index cdbee2c72d..2bd00076b2 100644 Binary files a/doc/ida/ida_examples.pdf and b/doc/ida/ida_examples.pdf differ diff --git a/doc/ida/ida_guide.pdf b/doc/ida/ida_guide.pdf index 48e708e395..371d860486 100644 Binary files a/doc/ida/ida_guide.pdf and b/doc/ida/ida_guide.pdf differ diff --git a/doc/idas/idas_examples.pdf b/doc/idas/idas_examples.pdf index ba0d8d4347..e601ff8716 100644 Binary files a/doc/idas/idas_examples.pdf and b/doc/idas/idas_examples.pdf differ diff --git a/doc/idas/idas_guide.pdf b/doc/idas/idas_guide.pdf index 4b72003d96..2cf26f9558 100644 Binary files a/doc/idas/idas_guide.pdf and b/doc/idas/idas_guide.pdf differ diff --git a/doc/kinsol/kin_examples.pdf b/doc/kinsol/kin_examples.pdf index 47d6d670be..97133def05 100644 Binary files a/doc/kinsol/kin_examples.pdf and b/doc/kinsol/kin_examples.pdf differ diff --git a/doc/kinsol/kin_guide.pdf b/doc/kinsol/kin_guide.pdf index 0d164150cd..25dbeb2938 100644 Binary files a/doc/kinsol/kin_guide.pdf and b/doc/kinsol/kin_guide.pdf differ diff --git a/examples/arkode/C_serial/CMakeLists.txt b/examples/arkode/C_serial/CMakeLists.txt index 88a0a585e8..0b7cbfcf23 100644 --- a/examples/arkode/C_serial/CMakeLists.txt +++ b/examples/arkode/C_serial/CMakeLists.txt @@ -15,29 +15,33 @@ # CMakeLists.txt file for ARKODE serial examples # --------------------------------------------------------------- -# Example lists are tuples "name\;type" where the type is +# Example lists are tuples "name\;args\;type" where the type is # 'develop' for examples excluded from 'make test' in releases # Examples using SUNDIALS linear solvers set(ARKODE_examples - "ark_analytic\;" - "ark_analytic_nonlin\;develop" - "ark_brusselator\;develop" - "ark_brusselator_fp\;develop" - "ark_brusselator1D\;develop" - "ark_heat1D\;develop" - "ark_heat1D_adapt\;develop" - "ark_KrylovDemo_prec\;develop" - "ark_robertson\;develop" - "ark_robertson_constraints\;develop" - "ark_robertson_root\;develop" - "ark_brusselator_mri\;develop" - "ark_onewaycouple_mri\;develop" - "ark_twowaycouple_mri\;develop" - "ark_reaction_diffusion_mri\;develop" - "ark_brusselator_1D_mri\;develop" + "ark_analytic\;\;" + "ark_analytic_nonlin\;\;develop" + "ark_brusselator\;\;develop" + "ark_brusselator_fp\;\;develop" + "ark_brusselator1D\;\;develop" + "ark_heat1D\;\;develop" + "ark_heat1D_adapt\;\;develop" + "ark_KrylovDemo_prec\;\;develop" + "ark_robertson\;\;develop" + "ark_robertson_constraints\;\;develop" + "ark_robertson_root\;\;develop" + "ark_brusselator_mri\;\;develop" + "ark_onewaycouple_mri\;\;develop" + "ark_twowaycouple_mri\;\;develop" + "ark_reaction_diffusion_mri\;\;develop" + "ark_brusselator_1D_mri\;\;develop" ) +if(SUNDIALS_BUILD_WITH_MONITORING) + list(APPEND ARKODE_examples "ark_brusselator_fp\;1\;develop") +endif() + # Examples using LAPACK linear solvers set(ARKODE_examples_BL ) @@ -81,22 +85,34 @@ foreach(example_tuple ${ARKODE_examples}) # parse the example tuple list(GET example_tuple 0 example) - list(GET example_tuple 1 example_type) + list(GET example_tuple 1 example_args) + list(GET example_tuple 2 example_type) + + if (NOT TARGET ${example}) + # example source files + add_executable(${example} ${example}.c) + + # folder for IDEs + set_target_properties(${example} PROPERTIES FOLDER "Examples") - # example source files - add_executable(${example} ${example}.c) + # libraries to link against + target_link_libraries(${example} ${SUNDIALS_LIBS}) + endif() - set_target_properties(${example} PROPERTIES FOLDER "Examples") + # check if example args are provided and set the test name + if("${example_args}" STREQUAL "") + set(test_name ${example}) + else() + string(REGEX REPLACE " " "_" test_name ${example}_${example_args}) + endif() # add example to regression tests - sundials_add_test(${example} ${example} + sundials_add_test(${test_name} ${example} + TEST_ARGS ${example_args} ANSWER_DIR ${CMAKE_CURRENT_SOURCE_DIR} - ANSWER_FILE ${example}.out + ANSWER_FILE ${test_name}.out EXAMPLE_TYPE ${example_type}) - # libraries to link against - target_link_libraries(${example} ${SUNDIALS_LIBS}) - # install example source and out files if(EXAMPLES_INSTALL) install(FILES ${example}.c ${example}.out diff --git a/examples/arkode/C_serial/ark_brusselator_fp.c b/examples/arkode/C_serial/ark_brusselator_fp.c index 150def1652..52e516a3c8 100644 --- a/examples/arkode/C_serial/ark_brusselator_fp.c +++ b/examples/arkode/C_serial/ark_brusselator_fp.c @@ -78,7 +78,7 @@ static int fi(realtype t, N_Vector y, N_Vector ydot, void *user_data); static int check_flag(void *flagvalue, const char *funcname, int opt); /* Main Program */ -int main() +int main(int argc, char *argv[]) { /* general problem parameters */ realtype T0 = RCONST(0.0); /* initial time */ @@ -93,17 +93,23 @@ int main() int maxcor = 10; /* maximum # of nonlinear iterations/step */ realtype a, b, ep, u0, v0, w0; realtype rdata[3]; + int monitor = 0; /* turn on/off monitoring */ /* general problem variables */ int flag; /* reusable error-checking flag */ N_Vector y = NULL; /* empty vector for storing solution */ SUNNonlinearSolver NLS = NULL; /* empty nonlinear solver object */ void *arkode_mem = NULL; /* empty ARKode memory structure */ - FILE *UFID; + FILE *UFID, *INFOFID; realtype t, tout; int iout; long int nst, nst_a, nfe, nfi, nni, ncfn, netf; + /* read inputs */ + if (argc == 2) { + monitor = atoi(argv[1]); + } + /* set up the test problem according to the desired test */ if (test == 1) { u0 = RCONST(3.9); @@ -134,6 +140,9 @@ int main() printf(" problem parameters: a = %"GSYM", b = %"GSYM", ep = %"GSYM"\n",a,b,ep); printf(" reltol = %.1"ESYM", abstol = %.1"ESYM"\n\n",reltol,abstol); + /* Open up info output file */ + if (monitor) INFOFID = fopen("ark_brusselator_fp-info.txt","w"); + /* Initialize data structures */ rdata[0] = a; /* set user data */ rdata[1] = b; @@ -153,6 +162,12 @@ int main() /* Initialize fixed-point nonlinear solver and attach to ARKStep */ NLS = SUNNonlinSol_FixedPoint(y, fp_m); if (check_flag((void *)NLS, "SUNNonlinSol_FixedPoint", 0)) return 1; + if (monitor) { + flag = SUNNonlinSolSetPrintLevel_FixedPoint(NLS, 1); + if (check_flag(&flag, "SUNNonlinSolSetPrintLevel_Newton", 1)) return(1); + flag = SUNNonlinSolSetInfoFile_FixedPoint(NLS, INFOFID); + if (check_flag(&flag, "SUNNonlinSolSetPrintLevel_Newton", 1)) return(1); + } flag = ARKStepSetNonlinearSolver(arkode_mem, NLS); if (check_flag(&flag, "ARKStepSetNonlinearSolver", 1)) return 1; @@ -196,6 +211,7 @@ int main() } printf(" ----------------------------------------------\n"); fclose(UFID); + if (monitor) fclose(INFOFID); /* Print some final statistics */ flag = ARKStepGetNumSteps(arkode_mem, &nst); diff --git a/examples/arkode/C_serial/ark_brusselator_fp_1.out b/examples/arkode/C_serial/ark_brusselator_fp_1.out new file mode 100644 index 0000000000..041ebe3d2a --- /dev/null +++ b/examples/arkode/C_serial/ark_brusselator_fp_1.out @@ -0,0 +1,27 @@ + +Brusselator ODE test problem, fixed-point solver: + initial conditions: u0 = 3, v0 = 3, w0 = 3.5 + problem parameters: a = 0.5, b = 3, ep = 0.0005 + reltol = 1.0e-06, abstol = 1.0e-10 + + t u v w + ---------------------------------------------- + 1.000000 1.897255 1.274939 2.997155 + 2.000000 0.346125 2.366448 2.999481 + 3.000000 0.147442 2.862061 2.999781 + 4.000000 0.140733 3.226731 2.999788 + 5.000000 0.142659 3.583206 2.999788 + 6.000000 0.145095 3.936910 2.999782 + 7.000000 0.147720 4.287893 2.999780 + 8.000000 0.150542 4.635957 2.999775 + 9.000000 0.153590 4.980863 2.999768 + 10.000000 0.156901 5.322330 2.999763 + ---------------------------------------------- + +Final Solver Statistics: + Internal solver steps = 729 (attempted = 730) + Total RHS evals: Fe = 4383, Fi = 18793 + Total number of fixed-point iterations = 14410 + Total number of nonlinear solver convergence failures = 0 + Total number of error test failures = 1 + diff --git a/examples/cvode/C_mpimanyvector/CMakeLists.txt b/examples/cvode/C_mpimanyvector/CMakeLists.txt index 7f3aa5dac2..e85a799677 100644 --- a/examples/cvode/C_mpimanyvector/CMakeLists.txt +++ b/examples/cvode/C_mpimanyvector/CMakeLists.txt @@ -38,6 +38,11 @@ else() set(NVECP_LIB sundials_nvecmpimanyvector_shared sundials_nvecparallel_shared) endif() +if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + list(APPEND CVODE_LIB + sundials_cvode_fused_stubs_${LINK_LIBRARY_TYPE}) +endif() + # Set-up linker flags and link libraries set(SUNDIALS_LIBS ${CVODE_LIB} ${NVECP_LIB} ${EXTRA_LINK_LIBS}) @@ -87,6 +92,9 @@ if(EXAMPLES_INSTALL) # Prepare substitution variables for Makefile and/or CMakeLists templates set(SOLVER "CVODE") set(SOLVER_LIB "sundials_cvode") + if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + set(LIBS "-lsundials_cvode_fused_stubs ${LIBS}") + endif() examples2string(CVODE_examples EXAMPLES) diff --git a/examples/cvode/C_openmp/CMakeLists.txt b/examples/cvode/C_openmp/CMakeLists.txt index 92939c1ebe..c6910b49ca 100644 --- a/examples/cvode/C_openmp/CMakeLists.txt +++ b/examples/cvode/C_openmp/CMakeLists.txt @@ -33,6 +33,11 @@ else() set(NVECOMP_LIB sundials_nvecopenmp_shared) endif() +if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + list(APPEND CVODE_LIB + sundials_cvode_fused_stubs_${LINK_LIBRARY_TYPE}) +endif() + # Set-up linker flags and link libraries set(SUNDIALS_LIBS ${CVODE_LIB} ${NVECOMP_LIB} ${EXTRA_LINK_LIBS}) @@ -81,6 +86,9 @@ if(EXAMPLES_INSTALL) # Prepare substitution variables for Makefile and/or CMakeLists templates set(SOLVER "CVODE") set(SOLVER_LIB "sundials_cvode") + if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + set(LIBS "-lsundials_cvode_fused_stubs ${LIBS}") + endif() examples2string(CVODE_examples EXAMPLES) diff --git a/examples/cvode/C_openmpdev/CMakeLists.txt b/examples/cvode/C_openmpdev/CMakeLists.txt index 2e5d496c6b..f3bc1b0a08 100644 --- a/examples/cvode/C_openmpdev/CMakeLists.txt +++ b/examples/cvode/C_openmpdev/CMakeLists.txt @@ -32,6 +32,11 @@ else() set(NVECOMP_LIB sundials_nvecopenmpdev_shared) endif() +if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + list(APPEND CVODE_LIB + sundials_cvode_fused_stubs_${LINK_LIBRARY_TYPE}) +endif() + # Set-up linker flags and link libraries set(SUNDIALS_LIBS ${CVODE_LIB} ${NVECOMP_LIB} ${EXTRA_LINK_LIBS}) @@ -81,6 +86,9 @@ if(EXAMPLES_INSTALL) # Prepare substitution variables for Makefile and/or CMakeLists templates set(SOLVER "CVODE") set(SOLVER_LIB "sundials_cvode") + if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + set(LIBS "-lsundials_cvode_fused_stubs ${LIBS}") + endif() examples2string(CVODE_examples EXAMPLES) diff --git a/examples/cvode/F2003_serial/CMakeLists.txt b/examples/cvode/F2003_serial/CMakeLists.txt index acca42ae3f..4b694d76b2 100644 --- a/examples/cvode/F2003_serial/CMakeLists.txt +++ b/examples/cvode/F2003_serial/CMakeLists.txt @@ -46,6 +46,12 @@ else() set(CVODE_LIB sundials_fcvode_mod_shared) endif() +if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + list(APPEND CVODE_LIB + sundials_cvode_${LINK_LIBRARY_TYPE} + sundials_cvode_fused_stubs_${LINK_LIBRARY_TYPE}) +endif() + # Set-up linker flags and link libraries set(SUNDIALS_LIBS ${CVODE_LIB} ${EXTRA_LINK_LIBS}) @@ -139,6 +145,9 @@ if(EXAMPLES_INSTALL) set(SOLVER_FLIB "sundials_fcvode_mod") set(NVEC_LIB "sundials_nvecserial") set(NVEC_FLIB "sundials_fnvecserial_mod") + if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + set(LIBS "-lsundials_cvode_fused_stubs ${LIBS}") + endif() examples2string(FCVODE_examples EXAMPLES) diff --git a/examples/cvode/cuda/CMakeLists.txt b/examples/cvode/cuda/CMakeLists.txt index cc0dac3b1a..2ee41ea0bb 100644 --- a/examples/cvode/cuda/CMakeLists.txt +++ b/examples/cvode/cuda/CMakeLists.txt @@ -21,13 +21,16 @@ # Examples using SUNDIALS linear solvers set(CVODE_examples - "cvAdvDiff_kry_cuda\;develop" - "cvAdvDiff_kry_cuda_managed\;develop" + "cvAdvDiff_kry_cuda\;\;develop" + "cvAdvDiff_kry_cuda_managed\;\;develop" + "cvAdvDiff_diag_cuda\;0 0\;develop" + "cvAdvDiff_diag_cuda\;0 1\;develop" + "cvAdvDiff_diag_cuda\;1 1\;develop" ) # Examples using cuSolverSP linear solvers set(CVODE_examples_cusolver - "cvRoberts_block_cusolversp_batchqr\;develop" + "cvRoberts_block_cusolversp_batchqr\;\;develop" ) if(SUNDIALS_INDEX_SIZE MATCHES "32") @@ -56,6 +59,11 @@ else() ${EXTRA_LINK_LIBS}) endif() +if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + set(SUNDIALS_LIBS ${SUNDIALS_LIBS} + sundials_cvode_fused_cuda_${LINK_LIBRARY_TYPE}) +endif() + # Add source directory to include directories include_directories(.) @@ -64,30 +72,41 @@ foreach(example_tuple ${all_examples}) # parse the example tuple list(GET example_tuple 0 example) - list(GET example_tuple 1 example_type) + list(GET example_tuple 1 example_args) + list(GET example_tuple 2 example_type) set_source_files_properties(${example}.cu PROPERTIES CUDA_SOURCE_PROPERTY_FORMAT OBJ) - # example source files - add_executable(${example} ${example}.cu) + if (NOT TARGET ${example}) + # example source files + add_executable(${example} ${example}.cu) - set_target_properties(${example} PROPERTIES FOLDER "Examples") + # folder for IDEs + set_target_properties(${example} PROPERTIES FOLDER "Examples") + + # libraries to link against + target_link_libraries(${example} PRIVATE ${SUNDIALS_LIBS}) + endif() + + # check if example args are provided and set the test name + if("${example_args}" STREQUAL "") + set(test_name ${example}) + else() + string(REGEX REPLACE " " "_" test_name ${example}_${example_args}) + endif() # add example to regression tests - sundials_add_test(${example} ${example} + sundials_add_test(${test_name} ${example} + TEST_ARGS ${example_args} ANSWER_DIR ${CMAKE_CURRENT_SOURCE_DIR} - ANSWER_FILE ${example}.out + ANSWER_FILE ${test_name}.out EXAMPLE_TYPE ${example_type}) - # libraries to link against - target_link_libraries(${example} PRIVATE ${SUNDIALS_LIBS}) - # install example source and out files if(EXAMPLES_INSTALL) - install(FILES ${example}.cu ${example}.out + install(FILES ${example}.cu ${test_name}.out DESTINATION ${EXAMPLES_INSTALL_PATH}/cvode/cuda) endif() - endforeach(example_tuple ${CVODE_examples}) @@ -101,11 +120,15 @@ if(EXAMPLES_INSTALL) set(SOLVER "CVODE") set(SOLVER_LIB "sundials_cvode") set(NVECTOR_LIB "sundials_nveccuda") + if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + set(LIBS "-lsundials_cvode_fused_cuda ${LIBS}") + endif() examples2string(CVODE_examples EXAMPLES) - if(SUNDIALS_INDEX_SIZE MATCHES "32" AND cuda_arch_ok) + if(SUNDIALS_INDEX_SIZE MATCHES "32") set(SUNLS_LIB "sundials_sunlinsolcusolversp") + set(SUNMAT_LIB "sundials_sunmatrixcusparse") examples2string(CVODE_examples_cusolver EXAMPLES_CUSOLVER) endif() diff --git a/examples/cvode/cuda/cvAdvDiff_diag_cuda.cu b/examples/cvode/cuda/cvAdvDiff_diag_cuda.cu new file mode 100644 index 0000000000..c9eba6ce72 --- /dev/null +++ b/examples/cvode/cuda/cvAdvDiff_diag_cuda.cu @@ -0,0 +1,364 @@ +/* + * ----------------------------------------------------------------- + * Programmer(s): Cody J. Balos @ LLNL + * ----------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2002-2020, Lawrence Livermore National Security + * and Southern Methodist University. + * All rights reserved. + * + * See the top-level LICENSE and NOTICE files for details. + * + * SPDX-License-Identifier: BSD-3-Clause + * SUNDIALS Copyright End + * ----------------------------------------------------------------- + * Example problem: + * + * The following is a simple example problem, with the program for + * its solution by CVODE. The problem is the semi-discrete + * form of the advection-diffusion equation in 1-D: + * du/dt = d^2 u / dx^2 + .5 du/dx + * on the interval 0 <= x <= 2, and the time interval 0 <= t <= 5. + * Homogeneous Dirichlet boundary conditions are posed, and the + * initial condition is the following: + * u(x,t=0) = x(2-x)exp(2x) . + * The PDE is discretized on a uniform grid of size MX+2 with + * central differencing, and with boundary values eliminated, + * leaving an ODE system of size NEQ = MX. + * This program solves the problem with the ADAMS integration method, + * and with Newton iteration using diagonal approximate Jacobians. + * It can use scalar (default) relative and absolute tolerances or a + * vector of absolute tolerances (controlled by a runtime argument). + * The constraint u_i >= 0 is posed for all components. + * Output is printed at t = .5, 1.0, ..., 5. + * Run statistics (optional outputs) are printed at the end. + * + * ./cvAdvDiff_diag_cuda [0 (scalar atol) | 1 (vector atol)] + * [0 (unfused) | 1 (fused)] + * ----------------------------------------------------------------- + */ + +#include +#include +#include + +#include + +#include /* prototypes for CVODE fcts., consts. */ +#include /* prototypes for CVODE diagonal solver */ +#include /* access to cuda N_Vector */ +#include /* definition of type realtype */ + +/* Problem Constants */ + +#define ZERO RCONST(0.0) + +#define XMAX RCONST(2.0) /* domain boundary */ +#define MX 10 /* mesh dimension */ +#define NEQ MX /* number of equations */ +#define ATOL RCONST(1e-10) /* scalar absolute tolerance */ +#define T0 ZERO /* initial time */ +#define T1 RCONST(0.5) /* first output time */ +#define DTOUT RCONST(0.5) /* output time increment */ +#define NOUT 10 /* number of output times */ + +/* Type : UserData + contains mesh spacing and problem parameters. */ + +typedef struct { + realtype dx; + realtype hdcoef; + realtype hacoef; +} *UserData; + +/* Private Helper Functions */ + +static void SetIC(N_Vector u, realtype dx); + +static void PrintIntro(int toltype, int usefused); + +static void PrintData(realtype t, realtype umax, long int nst); + +static void PrintFinalStats(void *cvode_mem); + +/* Functions Called by the Solver */ + +static int f(realtype t, N_Vector u, N_Vector udot, void *user_data); + +/* Private function to check function return values */ + +static int check_retval(void *returnvalue, const char *funcname, int opt); + +/***************************** Main Program ******************************/ + +int main(int argc, char *argv[]) +{ + realtype dx, reltol, abstol, t, tout, umax; + N_Vector u; + UserData data; + void *cvode_mem; + int iout, retval, toltype, usefused; + long int nst; + + u = NULL; + data = NULL; + cvode_mem = NULL; + toltype = 0; + usefused = 0; + + if (argc >= 2) { + /* use vector or scalar atol? */ + toltype = atoi(argv[1]); + /* use fused operations? */ + if (argc == 3) + usefused = atoi(argv[2]); + } + + data = (UserData) malloc(sizeof *data); /* Allocate data memory */ + if(check_retval((void *)data, "malloc", 2)) return 1; + + u = N_VNew_Cuda(NEQ); /* Allocate u vector */ + if(check_retval((void *)u, "N_VNew", 0)) return 1; + + reltol = ZERO; /* Set the tolerances */ + abstol = ATOL; + + dx = data->dx = XMAX/((realtype)(MX+1)); /* Set grid coefficients in data */ + data->hdcoef = RCONST(1.0)/(dx*dx); + data->hacoef = RCONST(0.5)/(RCONST(2.0)*dx); + + SetIC(u, dx); /* Initialize u vector */ + + /* Call CVodeCreate to create the solver memory and specify the + * Adams-Moulton LMM */ + cvode_mem = CVodeCreate(CV_ADAMS); + if(check_retval((void *)cvode_mem, "CVodeCreate", 0)) return 1; + + retval = CVodeSetUserData(cvode_mem, data); + if(check_retval(&retval, "CVodeSetUserData", 1)) return 1; + + /* Call CVodeInit to initialize the integrator memory and specify the + * user's right hand side function in u'=f(t,u), the inital time T0, and + * the initial dependent variable vector u. */ + retval = CVodeInit(cvode_mem, f, T0, u); + if(check_retval(&retval, "CVodeInit", 1)) return(1); + + /* Call CVodeSStolerances to specify the scalar relative tolerance + * and scalar absolute tolerances */ + + if (toltype == 0) { + retval = CVodeSStolerances(cvode_mem, reltol, abstol); + if (check_retval(&retval, "CVodeSStolerances", 1)) return(1); + } else { + N_Vector vabstol = N_VClone_Cuda(u); + if (check_retval(&vabstol, "N_VClone_Cuda", 0)) return(1); + N_VConst(abstol, vabstol); + retval = CVodeSVtolerances(cvode_mem, reltol, vabstol); + if (check_retval(&retval, "CVodeSVtolerances", 1)) return(1); + N_VDestroy(vabstol); + } + + /* Call CVDiag to create and attach CVODE-specific diagonal linear solver */ + retval = CVDiag(cvode_mem); + if(check_retval(&retval, "CVDiag", 1)) return(1); + + /* Tell CVode to use fused kernels if they are available. */ + retval = CVodeSetUseIntegratorFusedKernels(cvode_mem, usefused); + check_retval(&retval, "CVodeSetUseIntegratorFusedKernels", 1); + + PrintIntro(toltype, usefused); + + umax = N_VMaxNorm(u); + + t = T0; + PrintData(t, umax, 0); + + /* In loop over output points, call CVode, print results, test for error */ + + for (iout=1, tout=T1; iout <= NOUT; iout++, tout += DTOUT) { + retval = CVode(cvode_mem, tout, u, &t, CV_NORMAL); + if(check_retval(&retval, "CVode", 1)) break; + umax = N_VMaxNorm(u); + retval = CVodeGetNumSteps(cvode_mem, &nst); + check_retval(&retval, "CVodeGetNumSteps", 1); + PrintData(t, umax, nst); + } + + PrintFinalStats(cvode_mem); /* Print some final statistics */ + + N_VDestroy(u); /* Free the u vector */ + CVodeFree(&cvode_mem); /* Free the integrator memory */ + free(data); /* Free user data */ + + return(0); +} + +/************************ Private Helper Functions ***********************/ + +/* Set initial conditions in u vector */ + +static void SetIC(N_Vector u, realtype dx) +{ + int i; + sunindextype N; + realtype x; + realtype *udata; + + /* Set pointer to data array and get local length of u. */ + udata = N_VGetHostArrayPointer_Cuda(u); + N = N_VGetLength(u); + + /* Load initial profile into u vector */ + for (i=1; i<=N; i++) { + x = i*dx; + udata[i-1] = x*(XMAX - x)*exp(RCONST(2.0)*x); + } + N_VCopyToDevice_Cuda(u); +} + +/* Print problem introduction */ + +static void PrintIntro(int toltype, int usefused) +{ + printf("\n 1-D advection-diffusion equation, mesh size =%3d \n", MX); + printf("\n Diagonal linear solver CVDiag \n"); + if (usefused) + printf(" Using fused CVODE kernels \n"); + if (toltype == 0) + printf(" Using scalar ATOL\n"); + else + printf(" Using vector ATOL\n"); + printf("\n"); + + return; +} + +/* Print data */ + +static void PrintData(realtype t, realtype umax, long int nst) +{ + +#if defined(SUNDIALS_EXTENDED_PRECISION) + printf("At t = %4.2Lf max.norm(u) =%14.6Le nst =%4ld \n", t, umax, nst); +#elif defined(SUNDIALS_DOUBLE_PRECISION) + printf("At t = %4.2f max.norm(u) =%14.6e nst =%4ld \n", t, umax, nst); +#else + printf("At t = %4.2f max.norm(u) =%14.6e nst =%4ld \n", t, umax, nst); +#endif + + return; +} + +/* Print some final statistics located in the iopt array */ + +static void PrintFinalStats(void *cvode_mem) +{ + long int nst, nfe, nni, ncfn, netf; + int retval; + + retval = CVodeGetNumSteps(cvode_mem, &nst); + check_retval(&retval, "CVodeGetNumSteps", 1); + retval = CVodeGetNumRhsEvals(cvode_mem, &nfe); + check_retval(&retval, "CVodeGetNumRhsEvals", 1); + retval = CVodeGetNumErrTestFails(cvode_mem, &netf); + check_retval(&retval, "CVodeGetNumErrTestFails", 1); + retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni); + check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1); + retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn); + check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1); + + printf("\nFinal Statistics: \n\n"); + printf("nst = %-6ld nfe = %-6ld ", nst, nfe); + printf("nni = %-6ld ncfn = %-6ld netf = %ld\n \n", nni, ncfn, netf); +} + + /***************** Function Called by the Solver ***********************/ + + /* f routine. Compute f(t,u). */ + +__global__ +static void f_kernel(sunindextype N, + realtype hordc, realtype horac, + const realtype* u, realtype* udot) +{ + sunindextype i = blockDim.x*blockIdx.x + threadIdx.x; + realtype ui, ult, urt, hdiff, hadv; + + if (i < N) { + /* Extract u at x_i and two neighboring points */ + ui = u[i]; + ult = (i == 0) ? ZERO : u[i-1]; + urt = (i == N-1) ? ZERO : u[i+1]; + + /* Set diffusion and advection terms and load into udot */ + hdiff = hordc*(ult - RCONST(2.0)*ui + urt); + hadv = horac*(urt - ult); + udot[i] = hdiff + hadv; + } +} + +static int f(realtype t, N_Vector u, N_Vector udot, void *user_data) +{ + realtype hordc, horac; + realtype *udata, *dudata; + sunindextype N; + size_t grid, block; + UserData data; + cudaError_t cuerr; + + udata = N_VGetDeviceArrayPointer_Cuda(u); + dudata = N_VGetDeviceArrayPointer_Cuda(udot); + + /* Extract needed problem constants from data */ + data = (UserData) user_data; + hordc = data->hdcoef; + horac = data->hacoef; + + /* Extract parameters for parallel computation. */ + N = N_VGetLength(u); /* Number of elements of u. */ + + block = 64; + grid = (block + N - 1)/block; + f_kernel<<>>(N, hordc, horac, udata, dudata); + + cudaDeviceSynchronize(); + cuerr = cudaGetLastError(); + if (cuerr != cudaSuccess) { + fprintf(stderr, "ERROR in f: f_kernel --> %s\n", cudaGetErrorString(cuerr)); + return(-1); + } + + return(0); +} + + /* Check function return value... + opt == 0 means SUNDIALS function allocates memory so check if + returned NULL pointer + opt == 1 means SUNDIALS function returns an integer value so check if + retval < 0 + opt == 2 means function allocates memory so check if returned + NULL pointer */ + +static int check_retval(void *returnvalue, const char *funcname, int opt) +{ + int *retval; + + /* Check if SUNDIALS function returned NULL pointer - no memory allocated */ + if (opt == 0 && returnvalue == NULL) { + fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n", funcname); + return(1); } + + /* Check if retval < 0 */ + else if (opt == 1) { + retval = (int *) returnvalue; + if (*retval < 0) { + fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n", funcname, *retval); + return(1); }} + + /* Check if function returned NULL pointer - no memory allocated */ + else if (opt == 2 && returnvalue == NULL) { + fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n", funcname); + return(1); } + + return(0); +} \ No newline at end of file diff --git a/examples/cvode/cuda/cvAdvDiff_diag_cuda_0_0.out b/examples/cvode/cuda/cvAdvDiff_diag_cuda_0_0.out new file mode 100644 index 0000000000..211ab1796e --- /dev/null +++ b/examples/cvode/cuda/cvAdvDiff_diag_cuda_0_0.out @@ -0,0 +1,22 @@ + + 1-D advection-diffusion equation, mesh size = 10 + + Diagonal linear solver CVDiag + Using scalar ATOL + +At t = 0.00 max.norm(u) = 1.569909e+01 nst = 0 +At t = 0.50 max.norm(u) = 3.052879e+00 nst = 444 +At t = 1.00 max.norm(u) = 8.753297e-01 nst = 632 +At t = 1.50 max.norm(u) = 2.494935e-01 nst = 787 +At t = 2.00 max.norm(u) = 7.110094e-02 nst = 883 +At t = 2.50 max.norm(u) = 2.026233e-02 nst = 970 +At t = 3.00 max.norm(u) = 5.774352e-03 nst =1093 +At t = 3.50 max.norm(u) = 1.645572e-03 nst =1179 +At t = 4.00 max.norm(u) = 4.689201e-04 nst =1304 +At t = 4.50 max.norm(u) = 1.336304e-04 nst =1385 +At t = 5.00 max.norm(u) = 3.808141e-05 nst =1448 + +Final Statistics: + +nst = 1448 nfe = 2264 nni = 2261 ncfn = 34 netf = 56 + diff --git a/examples/cvode/cuda/cvAdvDiff_diag_cuda_0_1.out b/examples/cvode/cuda/cvAdvDiff_diag_cuda_0_1.out new file mode 100644 index 0000000000..1d8f94eeff --- /dev/null +++ b/examples/cvode/cuda/cvAdvDiff_diag_cuda_0_1.out @@ -0,0 +1,23 @@ + + 1-D advection-diffusion equation, mesh size = 10 + + Diagonal linear solver CVDiag + Using fused CVODE kernels + Using scalar ATOL + +At t = 0.00 max.norm(u) = 1.569909e+01 nst = 0 +At t = 0.50 max.norm(u) = 3.052879e+00 nst = 460 +At t = 1.00 max.norm(u) = 8.753297e-01 nst = 638 +At t = 1.50 max.norm(u) = 2.494935e-01 nst = 799 +At t = 2.00 max.norm(u) = 7.110094e-02 nst = 956 +At t = 2.50 max.norm(u) = 2.026233e-02 nst =1072 +At t = 3.00 max.norm(u) = 5.774352e-03 nst =1160 +At t = 3.50 max.norm(u) = 1.645568e-03 nst =1255 +At t = 4.00 max.norm(u) = 4.689520e-04 nst =1339 +At t = 4.50 max.norm(u) = 1.336373e-04 nst =1399 +At t = 5.00 max.norm(u) = 3.808158e-05 nst =1465 + +Final Statistics: + +nst = 1465 nfe = 2345 nni = 2342 ncfn = 21 netf = 70 + diff --git a/examples/cvode/cuda/cvAdvDiff_diag_cuda_1_1.out b/examples/cvode/cuda/cvAdvDiff_diag_cuda_1_1.out new file mode 100644 index 0000000000..9d602136cb --- /dev/null +++ b/examples/cvode/cuda/cvAdvDiff_diag_cuda_1_1.out @@ -0,0 +1,23 @@ + + 1-D advection-diffusion equation, mesh size = 10 + + Diagonal linear solver CVDiag + Using fused CVODE kernels + Using vector ATOL + +At t = 0.00 max.norm(u) = 1.569909e+01 nst = 0 +At t = 0.50 max.norm(u) = 3.052879e+00 nst = 460 +At t = 1.00 max.norm(u) = 8.753297e-01 nst = 638 +At t = 1.50 max.norm(u) = 2.494935e-01 nst = 799 +At t = 2.00 max.norm(u) = 7.110094e-02 nst = 956 +At t = 2.50 max.norm(u) = 2.026233e-02 nst =1072 +At t = 3.00 max.norm(u) = 5.774352e-03 nst =1160 +At t = 3.50 max.norm(u) = 1.645568e-03 nst =1255 +At t = 4.00 max.norm(u) = 4.689520e-04 nst =1339 +At t = 4.50 max.norm(u) = 1.336373e-04 nst =1399 +At t = 5.00 max.norm(u) = 3.808158e-05 nst =1465 + +Final Statistics: + +nst = 1465 nfe = 2345 nni = 2342 ncfn = 21 netf = 70 + diff --git a/examples/cvode/cuda/cvRoberts_block_cusolversp_batchqr.out b/examples/cvode/cuda/cvRoberts_block_cusolversp_batchqr.out index 45522fcff6..0f6d30f883 100644 --- a/examples/cvode/cuda/cvRoberts_block_cusolversp_batchqr.out +++ b/examples/cvode/cuda/cvRoberts_block_cusolversp_batchqr.out @@ -43,90 +43,90 @@ group 60: At t = 4.0000e+02 y = 4.505332e-01 3.223118e-06 5.494636e- group 70: At t = 4.0000e+02 y = 4.505332e-01 3.223118e-06 5.494636e-01 group 80: At t = 4.0000e+02 y = 4.505332e-01 3.223118e-06 5.494636e-01 group 90: At t = 4.0000e+02 y = 4.505332e-01 3.223118e-06 5.494636e-01 -group 0: At t = 4.0000e+03 y = 1.831798e-01 8.940768e-07 8.168193e-01 -group 10: At t = 4.0000e+03 y = 1.831798e-01 8.940768e-07 8.168193e-01 -group 20: At t = 4.0000e+03 y = 1.831798e-01 8.940768e-07 8.168193e-01 -group 30: At t = 4.0000e+03 y = 1.831798e-01 8.940768e-07 8.168193e-01 -group 40: At t = 4.0000e+03 y = 1.831798e-01 8.940768e-07 8.168193e-01 -group 50: At t = 4.0000e+03 y = 1.831798e-01 8.940768e-07 8.168193e-01 -group 60: At t = 4.0000e+03 y = 1.831798e-01 8.940768e-07 8.168193e-01 -group 70: At t = 4.0000e+03 y = 1.831798e-01 8.940768e-07 8.168193e-01 -group 80: At t = 4.0000e+03 y = 1.831798e-01 8.940768e-07 8.168193e-01 -group 90: At t = 4.0000e+03 y = 1.831798e-01 8.940768e-07 8.168193e-01 -group 0: At t = 4.0000e+04 y = 3.897912e-02 1.621575e-07 9.610207e-01 -group 10: At t = 4.0000e+04 y = 3.897912e-02 1.621575e-07 9.610207e-01 -group 20: At t = 4.0000e+04 y = 3.897912e-02 1.621575e-07 9.610207e-01 -group 30: At t = 4.0000e+04 y = 3.897912e-02 1.621575e-07 9.610207e-01 -group 40: At t = 4.0000e+04 y = 3.897912e-02 1.621575e-07 9.610207e-01 -group 50: At t = 4.0000e+04 y = 3.897912e-02 1.621575e-07 9.610207e-01 -group 60: At t = 4.0000e+04 y = 3.897912e-02 1.621575e-07 9.610207e-01 -group 70: At t = 4.0000e+04 y = 3.897912e-02 1.621575e-07 9.610207e-01 -group 80: At t = 4.0000e+04 y = 3.897912e-02 1.621575e-07 9.610207e-01 -group 90: At t = 4.0000e+04 y = 3.897912e-02 1.621575e-07 9.610207e-01 -group 0: At t = 4.0000e+05 y = 4.938822e-03 1.985206e-08 9.950612e-01 -group 10: At t = 4.0000e+05 y = 4.938822e-03 1.985206e-08 9.950612e-01 -group 20: At t = 4.0000e+05 y = 4.938822e-03 1.985206e-08 9.950612e-01 -group 30: At t = 4.0000e+05 y = 4.938822e-03 1.985206e-08 9.950612e-01 -group 40: At t = 4.0000e+05 y = 4.938822e-03 1.985206e-08 9.950612e-01 -group 50: At t = 4.0000e+05 y = 4.938822e-03 1.985206e-08 9.950612e-01 -group 60: At t = 4.0000e+05 y = 4.938822e-03 1.985206e-08 9.950612e-01 -group 70: At t = 4.0000e+05 y = 4.938822e-03 1.985206e-08 9.950612e-01 -group 80: At t = 4.0000e+05 y = 4.938822e-03 1.985206e-08 9.950612e-01 -group 90: At t = 4.0000e+05 y = 4.938822e-03 1.985206e-08 9.950612e-01 -group 0: At t = 4.0000e+06 y = 5.168467e-04 2.068440e-09 9.994832e-01 -group 10: At t = 4.0000e+06 y = 5.168467e-04 2.068440e-09 9.994832e-01 -group 20: At t = 4.0000e+06 y = 5.168467e-04 2.068440e-09 9.994832e-01 -group 30: At t = 4.0000e+06 y = 5.168467e-04 2.068440e-09 9.994832e-01 -group 40: At t = 4.0000e+06 y = 5.168467e-04 2.068440e-09 9.994832e-01 -group 50: At t = 4.0000e+06 y = 5.168467e-04 2.068440e-09 9.994832e-01 -group 60: At t = 4.0000e+06 y = 5.168467e-04 2.068440e-09 9.994832e-01 -group 70: At t = 4.0000e+06 y = 5.168467e-04 2.068440e-09 9.994832e-01 -group 80: At t = 4.0000e+06 y = 5.168467e-04 2.068440e-09 9.994832e-01 -group 90: At t = 4.0000e+06 y = 5.168467e-04 2.068440e-09 9.994832e-01 -group 0: At t = 4.0000e+07 y = 5.202426e-05 2.081078e-10 9.999480e-01 -group 10: At t = 4.0000e+07 y = 5.202426e-05 2.081078e-10 9.999480e-01 -group 20: At t = 4.0000e+07 y = 5.202426e-05 2.081078e-10 9.999480e-01 -group 30: At t = 4.0000e+07 y = 5.202426e-05 2.081078e-10 9.999480e-01 -group 40: At t = 4.0000e+07 y = 5.202426e-05 2.081078e-10 9.999480e-01 -group 50: At t = 4.0000e+07 y = 5.202426e-05 2.081078e-10 9.999480e-01 -group 60: At t = 4.0000e+07 y = 5.202426e-05 2.081078e-10 9.999480e-01 -group 70: At t = 4.0000e+07 y = 5.202426e-05 2.081078e-10 9.999480e-01 -group 80: At t = 4.0000e+07 y = 5.202426e-05 2.081078e-10 9.999480e-01 -group 90: At t = 4.0000e+07 y = 5.202426e-05 2.081078e-10 9.999480e-01 -group 0: At t = 4.0000e+08 y = 5.212615e-06 2.085057e-11 9.999948e-01 -group 10: At t = 4.0000e+08 y = 5.212615e-06 2.085057e-11 9.999948e-01 -group 20: At t = 4.0000e+08 y = 5.212615e-06 2.085057e-11 9.999948e-01 -group 30: At t = 4.0000e+08 y = 5.212615e-06 2.085057e-11 9.999948e-01 -group 40: At t = 4.0000e+08 y = 5.212615e-06 2.085057e-11 9.999948e-01 -group 50: At t = 4.0000e+08 y = 5.212615e-06 2.085057e-11 9.999948e-01 -group 60: At t = 4.0000e+08 y = 5.212615e-06 2.085057e-11 9.999948e-01 -group 70: At t = 4.0000e+08 y = 5.212615e-06 2.085057e-11 9.999948e-01 -group 80: At t = 4.0000e+08 y = 5.212615e-06 2.085057e-11 9.999948e-01 -group 90: At t = 4.0000e+08 y = 5.212615e-06 2.085057e-11 9.999948e-01 -group 0: At t = 4.0000e+09 y = 5.146550e-07 2.058622e-12 9.999995e-01 -group 10: At t = 4.0000e+09 y = 5.146550e-07 2.058622e-12 9.999995e-01 -group 20: At t = 4.0000e+09 y = 5.146550e-07 2.058622e-12 9.999995e-01 -group 30: At t = 4.0000e+09 y = 5.146550e-07 2.058622e-12 9.999995e-01 -group 40: At t = 4.0000e+09 y = 5.146550e-07 2.058622e-12 9.999995e-01 -group 50: At t = 4.0000e+09 y = 5.146550e-07 2.058622e-12 9.999995e-01 -group 60: At t = 4.0000e+09 y = 5.146550e-07 2.058622e-12 9.999995e-01 -group 70: At t = 4.0000e+09 y = 5.146550e-07 2.058622e-12 9.999995e-01 -group 80: At t = 4.0000e+09 y = 5.146550e-07 2.058622e-12 9.999995e-01 -group 90: At t = 4.0000e+09 y = 5.146550e-07 2.058622e-12 9.999995e-01 -group 0: At t = 4.0000e+10 y = 5.147690e-08 2.059076e-13 9.999999e-01 -group 10: At t = 4.0000e+10 y = 5.147690e-08 2.059076e-13 9.999999e-01 -group 20: At t = 4.0000e+10 y = 5.147690e-08 2.059076e-13 9.999999e-01 -group 30: At t = 4.0000e+10 y = 5.147690e-08 2.059076e-13 9.999999e-01 -group 40: At t = 4.0000e+10 y = 5.147690e-08 2.059076e-13 9.999999e-01 -group 50: At t = 4.0000e+10 y = 5.147690e-08 2.059076e-13 9.999999e-01 -group 60: At t = 4.0000e+10 y = 5.147690e-08 2.059076e-13 9.999999e-01 -group 70: At t = 4.0000e+10 y = 5.147690e-08 2.059076e-13 9.999999e-01 -group 80: At t = 4.0000e+10 y = 5.147690e-08 2.059076e-13 9.999999e-01 -group 90: At t = 4.0000e+10 y = 5.147690e-08 2.059076e-13 9.999999e-01 +group 0: At t = 4.0000e+03 y = 1.831797e-01 8.940765e-07 8.168194e-01 +group 10: At t = 4.0000e+03 y = 1.831797e-01 8.940765e-07 8.168194e-01 +group 20: At t = 4.0000e+03 y = 1.831797e-01 8.940765e-07 8.168194e-01 +group 30: At t = 4.0000e+03 y = 1.831797e-01 8.940765e-07 8.168194e-01 +group 40: At t = 4.0000e+03 y = 1.831797e-01 8.940765e-07 8.168194e-01 +group 50: At t = 4.0000e+03 y = 1.831797e-01 8.940765e-07 8.168194e-01 +group 60: At t = 4.0000e+03 y = 1.831797e-01 8.940765e-07 8.168194e-01 +group 70: At t = 4.0000e+03 y = 1.831797e-01 8.940765e-07 8.168194e-01 +group 80: At t = 4.0000e+03 y = 1.831797e-01 8.940765e-07 8.168194e-01 +group 90: At t = 4.0000e+03 y = 1.831797e-01 8.940765e-07 8.168194e-01 +group 0: At t = 4.0000e+04 y = 3.897911e-02 1.621575e-07 9.610207e-01 +group 10: At t = 4.0000e+04 y = 3.897911e-02 1.621575e-07 9.610207e-01 +group 20: At t = 4.0000e+04 y = 3.897911e-02 1.621575e-07 9.610207e-01 +group 30: At t = 4.0000e+04 y = 3.897911e-02 1.621575e-07 9.610207e-01 +group 40: At t = 4.0000e+04 y = 3.897911e-02 1.621575e-07 9.610207e-01 +group 50: At t = 4.0000e+04 y = 3.897911e-02 1.621575e-07 9.610207e-01 +group 60: At t = 4.0000e+04 y = 3.897911e-02 1.621575e-07 9.610207e-01 +group 70: At t = 4.0000e+04 y = 3.897911e-02 1.621575e-07 9.610207e-01 +group 80: At t = 4.0000e+04 y = 3.897911e-02 1.621575e-07 9.610207e-01 +group 90: At t = 4.0000e+04 y = 3.897911e-02 1.621575e-07 9.610207e-01 +group 0: At t = 4.0000e+05 y = 4.938738e-03 1.985178e-08 9.950612e-01 +group 10: At t = 4.0000e+05 y = 4.938738e-03 1.985178e-08 9.950612e-01 +group 20: At t = 4.0000e+05 y = 4.938738e-03 1.985178e-08 9.950612e-01 +group 30: At t = 4.0000e+05 y = 4.938738e-03 1.985178e-08 9.950612e-01 +group 40: At t = 4.0000e+05 y = 4.938738e-03 1.985178e-08 9.950612e-01 +group 50: At t = 4.0000e+05 y = 4.938738e-03 1.985178e-08 9.950612e-01 +group 60: At t = 4.0000e+05 y = 4.938738e-03 1.985178e-08 9.950612e-01 +group 70: At t = 4.0000e+05 y = 4.938738e-03 1.985178e-08 9.950612e-01 +group 80: At t = 4.0000e+05 y = 4.938738e-03 1.985178e-08 9.950612e-01 +group 90: At t = 4.0000e+05 y = 4.938738e-03 1.985178e-08 9.950612e-01 +group 0: At t = 4.0000e+06 y = 5.166270e-04 2.067563e-09 9.994834e-01 +group 10: At t = 4.0000e+06 y = 5.166270e-04 2.067563e-09 9.994834e-01 +group 20: At t = 4.0000e+06 y = 5.166270e-04 2.067563e-09 9.994834e-01 +group 30: At t = 4.0000e+06 y = 5.166270e-04 2.067563e-09 9.994834e-01 +group 40: At t = 4.0000e+06 y = 5.166270e-04 2.067563e-09 9.994834e-01 +group 50: At t = 4.0000e+06 y = 5.166270e-04 2.067563e-09 9.994834e-01 +group 60: At t = 4.0000e+06 y = 5.166270e-04 2.067563e-09 9.994834e-01 +group 70: At t = 4.0000e+06 y = 5.166270e-04 2.067563e-09 9.994834e-01 +group 80: At t = 4.0000e+06 y = 5.166270e-04 2.067563e-09 9.994834e-01 +group 90: At t = 4.0000e+06 y = 5.166270e-04 2.067563e-09 9.994834e-01 +group 0: At t = 4.0000e+07 y = 5.202720e-05 2.081194e-10 9.999480e-01 +group 10: At t = 4.0000e+07 y = 5.202720e-05 2.081194e-10 9.999480e-01 +group 20: At t = 4.0000e+07 y = 5.202720e-05 2.081194e-10 9.999480e-01 +group 30: At t = 4.0000e+07 y = 5.202720e-05 2.081194e-10 9.999480e-01 +group 40: At t = 4.0000e+07 y = 5.202720e-05 2.081194e-10 9.999480e-01 +group 50: At t = 4.0000e+07 y = 5.202720e-05 2.081194e-10 9.999480e-01 +group 60: At t = 4.0000e+07 y = 5.202720e-05 2.081194e-10 9.999480e-01 +group 70: At t = 4.0000e+07 y = 5.202720e-05 2.081194e-10 9.999480e-01 +group 80: At t = 4.0000e+07 y = 5.202720e-05 2.081194e-10 9.999480e-01 +group 90: At t = 4.0000e+07 y = 5.202720e-05 2.081194e-10 9.999480e-01 +group 0: At t = 4.0000e+08 y = 5.214602e-06 2.085852e-11 9.999948e-01 +group 10: At t = 4.0000e+08 y = 5.214602e-06 2.085852e-11 9.999948e-01 +group 20: At t = 4.0000e+08 y = 5.214602e-06 2.085852e-11 9.999948e-01 +group 30: At t = 4.0000e+08 y = 5.214602e-06 2.085852e-11 9.999948e-01 +group 40: At t = 4.0000e+08 y = 5.214602e-06 2.085852e-11 9.999948e-01 +group 50: At t = 4.0000e+08 y = 5.214602e-06 2.085852e-11 9.999948e-01 +group 60: At t = 4.0000e+08 y = 5.214602e-06 2.085852e-11 9.999948e-01 +group 70: At t = 4.0000e+08 y = 5.214602e-06 2.085852e-11 9.999948e-01 +group 80: At t = 4.0000e+08 y = 5.214602e-06 2.085852e-11 9.999948e-01 +group 90: At t = 4.0000e+08 y = 5.214602e-06 2.085852e-11 9.999948e-01 +group 0: At t = 4.0000e+09 y = 5.160314e-07 2.064127e-12 9.999995e-01 +group 10: At t = 4.0000e+09 y = 5.160314e-07 2.064127e-12 9.999995e-01 +group 20: At t = 4.0000e+09 y = 5.160314e-07 2.064127e-12 9.999995e-01 +group 30: At t = 4.0000e+09 y = 5.160314e-07 2.064127e-12 9.999995e-01 +group 40: At t = 4.0000e+09 y = 5.160314e-07 2.064127e-12 9.999995e-01 +group 50: At t = 4.0000e+09 y = 5.160314e-07 2.064127e-12 9.999995e-01 +group 60: At t = 4.0000e+09 y = 5.160314e-07 2.064127e-12 9.999995e-01 +group 70: At t = 4.0000e+09 y = 5.160314e-07 2.064127e-12 9.999995e-01 +group 80: At t = 4.0000e+09 y = 5.160314e-07 2.064127e-12 9.999995e-01 +group 90: At t = 4.0000e+09 y = 5.160314e-07 2.064127e-12 9.999995e-01 +group 0: At t = 4.0000e+10 y = 5.791263e-08 2.316505e-13 9.999999e-01 +group 10: At t = 4.0000e+10 y = 5.791263e-08 2.316505e-13 9.999999e-01 +group 20: At t = 4.0000e+10 y = 5.791263e-08 2.316505e-13 9.999999e-01 +group 30: At t = 4.0000e+10 y = 5.791263e-08 2.316505e-13 9.999999e-01 +group 40: At t = 4.0000e+10 y = 5.791263e-08 2.316505e-13 9.999999e-01 +group 50: At t = 4.0000e+10 y = 5.791263e-08 2.316505e-13 9.999999e-01 +group 60: At t = 4.0000e+10 y = 5.791263e-08 2.316505e-13 9.999999e-01 +group 70: At t = 4.0000e+10 y = 5.791263e-08 2.316505e-13 9.999999e-01 +group 80: At t = 4.0000e+10 y = 5.791263e-08 2.316505e-13 9.999999e-01 +group 90: At t = 4.0000e+10 y = 5.791263e-08 2.316505e-13 9.999999e-01 Final Statistics: -nst = 604 nfe = 859 nsetups = 121 nje = 12 -nni = 856 ncfn = 0 netf = 30 nge = 0 +nst = 547 nfe = 812 nsetups = 117 nje = 13 +nni = 809 ncfn = 0 netf = 29 nge = 0 cuSolverSp numerical factorization workspace size (in bytes) = 73472 -cuSolverSp internal Q, R buffer size (in bytes) = 9600 +cuSolverSp internal Q, R buffer size (in bytes) = 9600 \ No newline at end of file diff --git a/examples/cvode/fcmix_parallel/CMakeLists.txt b/examples/cvode/fcmix_parallel/CMakeLists.txt index 5bfa7e0bae..61a692396d 100644 --- a/examples/cvode/fcmix_parallel/CMakeLists.txt +++ b/examples/cvode/fcmix_parallel/CMakeLists.txt @@ -44,6 +44,11 @@ else() set(FNVECP_LIB sundials_fnvecparallel_shared) endif() +if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + list(APPEND CVODE_LIB + sundials_cvode_fused_stubs_${LINK_LIBRARY_TYPE}) +endif() + # Only static FCMIX libraries are available set(FCVODE_LIB sundials_fcvode_static) @@ -98,6 +103,9 @@ if(EXAMPLES_INSTALL) set(SOLVER "CVODE") set(SOLVER_LIB "sundials_cvode") set(SOLVER_FLIB "sundials_fcvode") + if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + set(LIBS "-lsundials_cvode_fused_stubs ${LIBS}") + endif() examples2string(FCVODE_examples EXAMPLES) diff --git a/examples/cvode/fcmix_serial/CMakeLists.txt b/examples/cvode/fcmix_serial/CMakeLists.txt index d0b0a9d5d7..04d7cb9314 100644 --- a/examples/cvode/fcmix_serial/CMakeLists.txt +++ b/examples/cvode/fcmix_serial/CMakeLists.txt @@ -60,6 +60,11 @@ else() set(FNVECS_LIB sundials_fnvecserial_shared) endif() +if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + list(APPEND CVODE_LIB + sundials_cvode_fused_stubs_${LINK_LIBRARY_TYPE}) +endif() + # Only static FCMIX libraries are available set(FCVODE_LIB sundials_fcvode_static) @@ -264,6 +269,9 @@ if(EXAMPLES_INSTALL) set(SOLVER "CVODE") set(SOLVER_LIB "sundials_cvode") set(SOLVER_FLIB "sundials_fcvode") + if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + set(LIBS "-lsundials_cvode_fused_stubs ${LIBS}") + endif() examples2string(FCVODE_examples EXAMPLES) diff --git a/examples/cvode/parallel/CMakeLists.txt b/examples/cvode/parallel/CMakeLists.txt index c4b4223c55..7cd7aad65a 100644 --- a/examples/cvode/parallel/CMakeLists.txt +++ b/examples/cvode/parallel/CMakeLists.txt @@ -41,6 +41,11 @@ else() set(NVECP_LIB sundials_nvecparallel_shared) endif() +if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + list(APPEND CVODE_LIB + sundials_cvode_fused_stubs_${LINK_LIBRARY_TYPE}) +endif() + # Set-up linker flags and link libraries set(SUNDIALS_LIBS ${CVODE_LIB} ${NVECP_LIB} ${EXTRA_LINK_LIBS}) @@ -90,6 +95,9 @@ if(EXAMPLES_INSTALL) # Prepare substitution variables for Makefile and/or CMakeLists templates set(SOLVER "CVODE") set(SOLVER_LIB "sundials_cvode") + if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + set(LIBS "-lsundials_cvode_fused_stubs ${LIBS}") + endif() examples2string(CVODE_examples EXAMPLES) diff --git a/examples/cvode/parhyp/CMakeLists.txt b/examples/cvode/parhyp/CMakeLists.txt index 61f7b62882..dea5f716ba 100644 --- a/examples/cvode/parhyp/CMakeLists.txt +++ b/examples/cvode/parhyp/CMakeLists.txt @@ -40,10 +40,14 @@ else() set(NVECP_LIB sundials_nvecparhyp_shared) endif() +if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + list(APPEND CVODE_LIB + sundials_cvode_fused_stubs_${LINK_LIBRARY_TYPE}) +endif() + # Set-up linker flags and link libraries set(SUNDIALS_LIBS ${CVODE_LIB} ${NVECP_LIB} ${EXTRA_LINK_LIBS}) - # Add the build and install targets for each example foreach(example_tuple ${CVODE_examples}) @@ -90,6 +94,9 @@ if(EXAMPLES_INSTALL) # Prepare substitution variables for Makefile and/or CMakeLists templates set(SOLVER "CVODE") set(SOLVER_LIB "sundials_cvode") + if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + set(LIBS "-lsundials_cvode_fused_stubs ${LIBS}") + endif() examples2string(CVODE_examples EXAMPLES) diff --git a/examples/cvode/petsc/CMakeLists.txt b/examples/cvode/petsc/CMakeLists.txt index 621d3d57da..3f368319ae 100644 --- a/examples/cvode/petsc/CMakeLists.txt +++ b/examples/cvode/petsc/CMakeLists.txt @@ -43,6 +43,11 @@ else() set(SUNNLS_LIB sundials_sunnonlinsolpetscsnes_shared) endif() +if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + list(APPEND CVODE_LIB + sundials_cvode_fused_stubs_${LINK_LIBRARY_TYPE}) +endif() + # Set-up linker flags and link libraries set(SUNDIALS_LIBS ${CVODE_LIB} ${NVECP_LIB} ${SUNNLS_LIB} ${EXTRA_LINK_LIBS}) @@ -90,6 +95,9 @@ if(EXAMPLES_INSTALL) # Prepare substitution variables for Makefile and/or CMakeLists templates set(SOLVER "CVODE") set(SOLVER_LIB "sundials_cvode") + if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + set(LIBS "-lsundials_cvode_fused_stubs ${LIBS}") + endif() examples2string(CVODE_examples EXAMPLES) diff --git a/examples/cvode/raja/CMakeLists.txt b/examples/cvode/raja/CMakeLists.txt index 32f49729fc..9c759e3f3c 100644 --- a/examples/cvode/raja/CMakeLists.txt +++ b/examples/cvode/raja/CMakeLists.txt @@ -38,6 +38,11 @@ else() set(NVECS_LIB sundials_nveccudaraja_shared) endif() +if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + list(APPEND CVODE_LIB + sundials_cvode_fused_stubs_${LINK_LIBRARY_TYPE}) +endif() + set(SUNDIALS_LIBS ${CVODE_LIB} ${NVECS_LIB} ${EXTRA_LINK_LIBS}) # Add the build and install targets for each CVODE example @@ -80,6 +85,9 @@ if(EXAMPLES_INSTALL) set(SOLVER "CVODE") set(SOLVER_LIB "sundials_cvode") set(NVECTOR_LIB "sundials_nveccudaraja") + if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + set(LIBS "-lsundials_cvode_fused_stubs ${LIBS}") + endif() examples2string(CVODE_examples EXAMPLES) diff --git a/examples/cvode/serial/CMakeLists.txt b/examples/cvode/serial/CMakeLists.txt index ab17d6a821..aae4fbbff8 100644 --- a/examples/cvode/serial/CMakeLists.txt +++ b/examples/cvode/serial/CMakeLists.txt @@ -15,40 +15,52 @@ # CMakeLists.txt file for CVODE serial examples # --------------------------------------------------------------- -# Example lists are tuples "name\;type" where the type is +# Example lists are tuples "name\;args\;type" where the type is # 'develop' for examples excluded from 'make test' in releases # Examples using SUNDIALS linear solvers set(CVODE_examples - "cvRoberts_dns\;" - "cvRoberts_dns_uw\;develop" - "cvRoberts_dns_negsol\;develop" - "cvRoberts_dns_constraints\;develop" - "cvAdvDiff_bnd\;develop" - "cvDirectDemo_ls\;develop" - "cvDiurnal_kry_bp\;develop" - "cvDiurnal_kry\;develop" - "cvDisc_dns\;develop" - "cvKrylovDemo_ls\;develop" - "cvKrylovDemo_prec\;develop" + "cvAdvDiff_bnd\;\;develop" + "cvDirectDemo_ls\;\;develop" + "cvDiurnal_kry_bp\;\;develop" + "cvDiurnal_kry\;\;develop" + "cvDisc_dns\;\;develop" + "cvKrylovDemo_ls\;\;develop" + "cvKrylovDemo_prec\;\;develop" + "cvParticle_dns\;\;develop" + "cvPendulum_dns\;\;develop" + "cvRoberts_dns\;\;" + "cvRoberts_dns_uw\;\;develop" + "cvRoberts_dns_negsol\;\;develop" + "cvRoberts_dns_constraints\;\;develop" #cvAdvDiffReac_kry\;develop" # not released ) +if(SUNDIALS_BUILD_WITH_MONITORING) + list(APPEND CVODE_examples "cvKrylovDemo_ls\;1\;develop") +endif() + # Examples using LAPACK linear solvers set(CVODE_examples_BL - "cvAdvDiff_bndL\;develop" - "cvRoberts_dnsL\;develop" + "cvAdvDiff_bndL\;\;develop" + "cvRoberts_dnsL\;\;develop" ) # Examples using KLU linear solver set(CVODE_examples_KLU - "cvRoberts_klu\;develop" - "cvRoberts_block_klu\;develop" + "cvRoberts_klu\;\;develop" + "cvRoberts_block_klu\;\;develop" ) # Examples using SuperLU_MT linear solver set(CVODE_examples_SUPERLUMT - "cvRoberts_sps\;develop" + "cvRoberts_sps\;\;develop" + ) + +# Auxiliary files to install +set(CVODE_extras + plot_cvParticle.py + plot_cvPendulum.py ) # Specify libraries to link against (through the target that was used to @@ -61,34 +73,55 @@ else() set(NVECS_LIB sundials_nvecserial_shared) endif() +if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + list(APPEND CVODE_LIB + sundials_cvode_fused_stubs_${LINK_LIBRARY_TYPE}) +endif() + # Set-up linker flags and link libraries set(SUNDIALS_LIBS ${CVODE_LIB} ${NVECS_LIB} ${EXTRA_LINK_LIBS}) - # Add the build and install targets for each example foreach(example_tuple ${CVODE_examples}) # parse the example tuple list(GET example_tuple 0 example) - list(GET example_tuple 1 example_type) + list(GET example_tuple 1 example_args) + list(GET example_tuple 2 example_type) + + # check if this example has already been added, only need to add + # example source files once for testing with different inputs + if(NOT TARGET ${example}) + # example source files + add_executable(${example} ${example}.c) + + # folder to organize targets in an IDE + set_target_properties(${example} PROPERTIES FOLDER "Examples") - # example source files - add_executable(${example} ${example}.c) + # libraries to link against + target_link_libraries(${example} ${SUNDIALS_LIBS}) + endif() - set_target_properties(${example} PROPERTIES FOLDER "Examples") + # check if example args are provided and set the test name + if("${example_args}" STREQUAL "") + set(test_name ${example}) + else() + string(REGEX REPLACE " " "_" test_name ${example}_${example_args}) + endif() # add example to regression tests - sundials_add_test(${example} ${example} + sundials_add_test(${test_name} ${example} + TEST_ARGS ${example_args} ANSWER_DIR ${CMAKE_CURRENT_SOURCE_DIR} - ANSWER_FILE ${example}.out + ANSWER_FILE ${test_name}.out EXAMPLE_TYPE ${example_type}) - # libraries to link against - target_link_libraries(${example} ${SUNDIALS_LIBS}) + # find all .out files for this example + file(GLOB example_out ${example}*.out) - # install example source and out files + # install example source and .out files if(EXAMPLES_INSTALL) - install(FILES ${example}.c ${example}.out + install(FILES ${example}.c ${example_out} DESTINATION ${EXAMPLES_INSTALL_PATH}/cvode/serial) endif() @@ -121,25 +154,42 @@ if(LAPACK_FOUND) # parse the example tuple list(GET example_tuple 0 example) - list(GET example_tuple 1 example_type) + list(GET example_tuple 1 example_args) + list(GET example_tuple 2 example_type) - # example source files - add_executable(${example} ${example}.c) + # check if this example has already been added, only need to add + # example source files once for testing with different inputs + if(NOT TARGET ${example}) + # example source files + add_executable(${example} ${example}.c) - set_target_properties(${example} PROPERTIES FOLDER "Examples") + # folder to organize targets in an IDE + set_target_properties(${example} PROPERTIES FOLDER "Examples") + + # libraries to link against + target_link_libraries(${example} ${SUNDIALS_LIBS} ${SUNLINSOLLAPACK_LIBS}) + endif() + + # check if example args are provided and set the test name + if("${example_args}" STREQUAL "") + set(test_name ${example}) + else() + string(REGEX REPLACE " " "_" test_name ${example}_${example_args}) + endif() # add example to regression tests - sundials_add_test(${example} ${example} + sundials_add_test(${test_name} ${example} + TEST_ARGS ${example_args} ANSWER_DIR ${CMAKE_CURRENT_SOURCE_DIR} - ANSWER_FILE ${example}.out + ANSWER_FILE ${test_name}.out EXAMPLE_TYPE ${example_type}) - # libraries to link against - target_link_libraries(${example} ${SUNDIALS_LIBS} ${SUNLINSOLLAPACK_LIBS}) + # find all .out files for this example + file(GLOB example_out ${example}*.out) - # install example source and out files + # install example source and .out files if(EXAMPLES_INSTALL) - install(FILES ${example}.c ${example}.out + install(FILES ${example}.c ${example_out} DESTINATION ${EXAMPLES_INSTALL_PATH}/cvode/serial) endif() @@ -165,25 +215,42 @@ if(KLU_FOUND) # parse the example tuple list(GET example_tuple 0 example) - list(GET example_tuple 1 example_type) + list(GET example_tuple 1 example_args) + list(GET example_tuple 2 example_type) - # example source files - add_executable(${example} ${example}.c) + # check if this example has already been added, only need to add + # example source files once for testing with different inputs + if(NOT TARGET ${example}) + # add example source files + add_executable(${example} ${example}.c) - set_target_properties(${example} PROPERTIES FOLDER "Examples") + # folder to organize targets in an IDE + set_target_properties(${example} PROPERTIES FOLDER "Examples") + + # libraries to link against + target_link_libraries(${example} ${SUNDIALS_LIBS} ${SUNLINSOLKLU_LIBS}) + endif() + + # check if example args are provided and set the test name + if("${example_args}" STREQUAL "") + set(test_name ${example}) + else() + string(REGEX REPLACE " " "_" test_name ${example}_${example_args}) + endif() # add example to regression tests - sundials_add_test(${example} ${example} + sundials_add_test(${test_name} ${example} + TEST_ARGS ${example_args} ANSWER_DIR ${CMAKE_CURRENT_SOURCE_DIR} - ANSWER_FILE ${example}.out + ANSWER_FILE ${test_name}.out EXAMPLE_TYPE ${example_type}) - # libraries to link against - target_link_libraries(${example} ${SUNDIALS_LIBS} ${SUNLINSOLKLU_LIBS}) + # find all .out files for this example + file(GLOB example_out ${example}*.out) - # install example source and out files + # install example source and .out files if(EXAMPLES_INSTALL) - install(FILES ${example}.c ${example}.out + install(FILES ${example}.c ${example_out} DESTINATION ${EXAMPLES_INSTALL_PATH}/cvode/serial) endif() @@ -214,25 +281,35 @@ if(SUPERLUMT_FOUND) # parse the example tuple list(GET example_tuple 0 example) - list(GET example_tuple 1 example_type) + list(GET example_tuple 1 example_args) + list(GET example_tuple 2 example_type) - # example source files - add_executable(${example} ${example}.c) + # check if this example has already been added, only need to add + # example source files once for testing with different inputs + if(NOT TARGET ${example}) + # add example source files + add_executable(${example} ${example}.c) - set_target_properties(${example} PROPERTIES FOLDER "Examples") + # folder to organize targets in an IDE + set_target_properties(${example} PROPERTIES FOLDER "Examples") - # add example to regression tests - sundials_add_test(${example} ${example} - ANSWER_DIR ${CMAKE_CURRENT_SOURCE_DIR} - ANSWER_FILE ${example}.out - EXAMPLE_TYPE ${example_type}) + # libraries to link against + target_link_libraries(${example} ${SUNDIALS_LIBS} ${SUNLINSOLSLUMT_LIBS}) + endif() - # libraries to link against - target_link_libraries(${example} ${SUNDIALS_LIBS} ${SUNLINSOLSLUMT_LIBS}) + # check if example args are provided and set the test name + if("${example_args}" STREQUAL "") + set(test_name ${example}) + else() + string(REGEX REPLACE " " "_" test_name ${example}_${example_args}) + endif() + + # find all .out files for this example + file(GLOB example_out ${example}*.out) - # install example source and out files + # install example source and .out files if(EXAMPLES_INSTALL) - install(FILES ${example}.c ${example}.out + install(FILES ${example}.c ${example_out} DESTINATION ${EXAMPLES_INSTALL_PATH}/cvode/serial) endif() @@ -247,9 +324,17 @@ if(EXAMPLES_INSTALL) # Install the README file install(FILES README DESTINATION ${EXAMPLES_INSTALL_PATH}/cvode/serial) + # Install the extra files + foreach(extrafile ${CVODE_extras}) + install(FILES ${extrafile} DESTINATION ${EXAMPLES_INSTALL_PATH}/cvode/serial) + endforeach() + # Prepare substitution variables for Makefile and/or CMakeLists templates set(SOLVER "CVODE") set(SOLVER_LIB "sundials_cvode") + if(SUNDIALS_BUILD_PACKAGE_FUSED_KERNELS) + set(LIBS "-lsundials_cvode_fused_stubs ${LIBS}") + endif() examples2string(CVODE_examples EXAMPLES) diff --git a/examples/cvode/serial/cvKrylovDemo_ls.c b/examples/cvode/serial/cvKrylovDemo_ls.c index cfa96d388e..8cf582eee1 100644 --- a/examples/cvode/serial/cvKrylovDemo_ls.c +++ b/examples/cvode/serial/cvKrylovDemo_ls.c @@ -1,6 +1,6 @@ /* ----------------------------------------------------------------- * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh and - * Radu Serban @ LLNL + * Radu Serban, Cody J. Balos @ LLNL * ----------------------------------------------------------------- * SUNDIALS Copyright Start * Copyright (c) 2002-2020, Lawrence Livermore National Security @@ -13,7 +13,7 @@ * SUNDIALS Copyright End * ----------------------------------------------------------------- * This example loops through the available iterative linear solvers: - * SPGMR, SPBCG and SPTFQMR. + * SPGMR, SPFGMR, SPBCG and SPTFQMR. * * Example problem: * @@ -32,10 +32,10 @@ * 0 <= t <= 86400 sec (1 day). * The PDE system is treated by central differences on a uniform * 10 x 10 mesh, with simple polynomial initial profiles. - * The problem is solved with CVODE, with the BDF/GMRES, + * The problem is solved with CVODE, with the BDF/GMRES, BDF/FGMRES * BDF/Bi-CGStab, and BDF/TFQMR methods (i.e. using the SUNLinSol_SPGMR, - * SUNLinSol_SPBCGS and SUNLinSol_SPTFQMR linear solvers) and the - * block-diagonal part of the Newton matrix as a left preconditioner. + * SUNLinSol_SPFGMR, SUNLinSol_SPBCGS, and SUNLinSol_SPTFQMR linear solvers) + * and the block-diagonal part of the Newton matrix as a left preconditioner. * A copy of the block-diagonal part of the Jacobian is saved and * conditionally reused within the Precond routine. * -----------------------------------------------------------------*/ @@ -44,13 +44,15 @@ #include #include -#include /* main integrator header file */ -#include /* access to SPGMR SUNLinearSolver */ -#include /* access to SPBCGS SUNLinearSolver */ -#include /* access to SPTFQMR SUNLinearSolver */ -#include /* serial N_Vector types, fct. and macros */ -#include /* use generic DENSE solver in preconditioning */ -#include /* definition of realtype */ +#include /* main integrator header file */ +#include /* access to SPGMR SUNLinearSolver */ +#include /* access to SPFGMR SUNLinearSolver */ +#include /* access to SPBCGS SUNLinearSolver */ +#include /* access to SPTFQMR SUNLinearSolver */ +#include /* access to Newton SUNNonlinearSolver */ +#include /* serial N_Vector types, fct. and macros */ +#include /* use generic DENSE solver in preconditioning */ +#include /* definition of realtype */ /* helpful macros */ @@ -105,8 +107,9 @@ /* Linear Solver Loop Constants */ #define USE_SPGMR 0 -#define USE_SPBCG 1 -#define USE_SPTFQMR 2 +#define USE_SPFGMR 1 +#define USE_SPBCG 2 +#define USE_SPTFQMR 3 /* User-defined vector and matrix accessor macros: IJKth, IJth */ @@ -132,22 +135,25 @@ #define IJth(a,i,j) (a[j-1][i-1]) /* Type : UserData - contains preconditioner blocks, pivot arrays, and problem constants */ + contains preconditioner blocks, pivot arrays, and problem constants, + solution vector, and linsolver type */ typedef struct { realtype **P[MX][MY], **Jbd[MX][MY]; sunindextype *pivot[MX][MY]; realtype q4, om, dx, dy, hdco, haco, vdco; + N_Vector u; + int linsolver; } *UserData; /* Private Helper Functions */ static UserData AllocUserData(void); -static void InitUserData(UserData data); +static void InitUserData(UserData data, N_Vector u); static void FreeUserData(UserData data); static void SetInitialProfiles(N_Vector u, realtype dx, realtype dy); static void PrintOutput(void *cvode_mem, N_Vector u, realtype t); -static void PrintFinalStats(void *cvode_mem, int linsolver); +static void PrintStats(void *cvode_mem, int linsolver, int stats); static int check_retval(void *returnvalue, const char *funcname, int opt); /* Functions Called by the Solver */ @@ -162,6 +168,7 @@ static int PSolve(realtype tn, N_Vector u, N_Vector fu, realtype gamma, realtype delta, int lr, void *user_data); +static int myMonitorFunction(void *cvode_mem, void *user_data); /* *------------------------------- @@ -169,26 +176,40 @@ static int PSolve(realtype tn, N_Vector u, N_Vector fu, *------------------------------- */ -int main(void) +int main(int argc, char* argv[]) { realtype abstol, reltol, t, tout; N_Vector u; UserData data; SUNLinearSolver LS; + SUNNonlinearSolver NLS; void *cvode_mem; int linsolver, iout, retval; + FILE* infofp; + int monitor; u = NULL; data = NULL; LS = NULL; cvode_mem = NULL; + monitor = 0; + + if (argc == 2) { + monitor = atoi(argv[1]); + } + + /* Open info file if monitoring is turned on */ + if (monitor) { + infofp = fopen("cvKrylovDemo_ls-info.txt", "w+"); + if (check_retval((void *)infofp, "fopen", 0)) return(1); + } /* Allocate memory, and set problem data, initial values, tolerances */ u = N_VNew_Serial(NEQ); - if(check_retval((void *)u, "N_VNew_Serial", 0)) return(1); + if (check_retval((void *)u, "N_VNew_Serial", 0)) return(1); data = AllocUserData(); - if(check_retval((void *)data, "AllocUserData", 2)) return(1); - InitUserData(data); + if (check_retval((void *)data, "AllocUserData", 2)) return(1); + InitUserData(data, u); SetInitialProfiles(u, data->dx, data->dy); abstol=ATOL; reltol=RTOL; @@ -196,30 +217,55 @@ int main(void) /* Call CVodeCreate to create the solver memory and specify the * Backward Differentiation Formula */ cvode_mem = CVodeCreate(CV_BDF); - if(check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1); + if (check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1); /* Set the pointer to user-defined data */ retval = CVodeSetUserData(cvode_mem, data); - if(check_retval(&retval, "CVodeSetUserData", 1)) return(1); + if (check_retval(&retval, "CVodeSetUserData", 1)) return(1); /* Call CVodeInit to initialize the integrator memory and specify the * user's right hand side function in u'=f(t,u), the inital time T0, and * the initial dependent variable vector u. */ retval = CVodeInit(cvode_mem, f, T0, u); - if(check_retval(&retval, "CVodeInit", 1)) return(1); + if (check_retval(&retval, "CVodeInit", 1)) return(1); /* Call CVodeSStolerances to specify the scalar relative tolerance * and scalar absolute tolerances */ retval = CVodeSStolerances(cvode_mem, reltol, abstol); if (check_retval(&retval, "CVodeSStolerances", 1)) return(1); - /* START: Loop through SPGMR, SPBCG and SPTFQMR linear solver modules */ - for (linsolver = 0; linsolver < 3; ++linsolver) { + /* Set a function that CVode will call every 50 successful time steps. + * This will be used to monitor the solution and integrator statistics. */ + if (monitor) { + retval = CVodeSetMonitorFn(cvode_mem, myMonitorFunction); + if (check_retval(&retval, "CVodeSetMonitorFn", 1)) return(1); + retval = CVodeSetMonitorFrequency(cvode_mem, 50); + if (check_retval(&retval, "CVodeSetMonitorFrequency", 1)) return(1); + } + + /* Create the SUNNonlinearSolver */ + NLS = SUNNonlinSol_Newton(u); + if (check_retval(&retval, "SUNNonlinSol_Newton", 0)) return(1); + if (monitor) { + /* Set the print level set to 1, so that the nonlinear residual + is printed every newton iteration. */ + retval = SUNNonlinSolSetPrintLevel_Newton(NLS, 1); + if (check_retval(&retval, "SUNNonlinSolSetPrintLevel_Newton", 1)) return(1); + retval = SUNNonlinSolSetInfoFile_Newton(NLS, infofp); + if (check_retval(&retval, "SUNNonlinSolSetInfoFile_Newton", 1)) return(1); + } + + /* Call CVodeSetNonlinearSolver to attach the nonlinear solver to CVode */ + retval = CVodeSetNonlinearSolver(cvode_mem, NLS); + if (check_retval(&retval, "CVodeSetNonlinearSolver", 1)) return(1); + + /* START: Loop through SPGMR, SPFGMR, SPBCG and SPTFQMR linear solver modules */ + for (linsolver = 0; linsolver < 4; ++linsolver) { if (linsolver != 0) { /* Re-initialize user data */ - InitUserData(data); + InitUserData(data, u); SetInitialProfiles(u, data->dx, data->dy); /* Re-initialize CVode for the solution of the same problem, but @@ -232,6 +278,9 @@ int main(void) /* Free previous linear solver and attach a new linear solver module */ SUNLinSolFree(LS); + /* Set the linear sovler type in user data */ + data->linsolver = linsolver; + switch(linsolver) { /* (a) SPGMR */ @@ -241,74 +290,132 @@ int main(void) printf(" -------"); printf(" \n| SPGMR |\n"); printf(" -------\n"); + if (monitor) { + fprintf(infofp, " ---------"); + fprintf(infofp, " \n| SPGMR |\n"); + fprintf(infofp, " ---------\n"); + } /* Call SUNLinSol_SPGMR to specify the linear solver SPGMR with left preconditioning and the default maximum Krylov dimension */ LS = SUNLinSol_SPGMR(u, PREC_LEFT, 0); - if(check_retval((void *)LS, "SUNLinSol_SPGMR", 0)) return(1); + if (check_retval((void *)LS, "SUNLinSol_SPGMR", 0)) return(1); + if (monitor) { + retval = SUNLinSolSetPrintLevel_SPGMR(LS, 1); + if (check_retval(&retval, "SUNLinSolSetPrintLevel_SPGMR", 1)) return(1); + retval = SUNLinSolSetInfoFile_SPGMR(LS, infofp); + if (check_retval(&retval, "SUNLinSolSetInfoFile_SPGMR", 1)) return(1); + } + retval = CVodeSetLinearSolver(cvode_mem, LS, NULL); + if (check_retval(&retval, "CVodeSetLinearSolver", 1)) return 1; + + break; + + /* (b) SPFGMR */ + case(USE_SPFGMR): + + /* Print header */ + printf(" ---------"); + printf(" \n| SPFGMR |\n"); + printf(" ---------\n"); + if (monitor) { + fprintf(infofp, " ---------"); + fprintf(infofp, " \n| SPFGMR |\n"); + fprintf(infofp, " ---------\n"); + } + /* Call SUNLinSol_SPFGMR to specify the linear solver SPFGMR with + left preconditioning and the default maximum Krylov dimension */ + LS = SUNLinSol_SPFGMR(u, PREC_LEFT, 0); + if (check_retval((void *)LS, "SUNLinSol_SPFGMR", 0)) return(1); + if (monitor) { + retval = SUNLinSolSetPrintLevel_SPFGMR(LS, 1); + if (check_retval(&retval, "SUNLinSolSetPrintLevel_SPFGMR", 1)) return(1); + retval = SUNLinSolSetInfoFile_SPFGMR(LS, infofp); + if (check_retval(&retval, "SUNLinSolSetInfoFile_SPFGMR", 1)) return(1); + } retval = CVodeSetLinearSolver(cvode_mem, LS, NULL); - if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return 1; + if (check_retval(&retval, "CVodeSetLinearSolver", 1)) return 1; break; - /* (b) SPBCG */ + /* (c) SPBCG */ case(USE_SPBCG): /* Print header */ printf(" -------"); printf(" \n| SPBCGS |\n"); printf(" -------\n"); + if (monitor) { + fprintf(infofp, " ---------"); + fprintf(infofp, " \n| SPBCGS |\n"); + fprintf(infofp, " ---------\n"); + } /* Call SUNLinSol_SPBCGS to specify the linear solver SPBCGS with left preconditioning and the default maximum Krylov dimension */ LS = SUNLinSol_SPBCGS(u, PREC_LEFT, 0); - if(check_retval((void *)LS, "SUNLinSol_SPBCGS", 0)) return(1); - + if (check_retval((void *)LS, "SUNLinSol_SPBCGS", 0)) return(1); + if (monitor) { + retval = SUNLinSolSetPrintLevel_SPBCGS(LS, 1); + if (check_retval(&retval, "SUNLinSolSetPrintLevel_SPBCGS", 1)) return(1); + retval = SUNLinSolSetInfoFile_SPBCGS(LS, infofp); + if (check_retval(&retval, "SUNLinSolSetInfoFile_SPBCGS", 1)) return(1); + } retval = CVodeSetLinearSolver(cvode_mem, LS, NULL); - if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return 1; + if (check_retval(&retval, "CVodeSetLinearSolver", 1)) return 1; break; - /* (c) SPTFQMR */ + /* (d) SPTFQMR */ case(USE_SPTFQMR): /* Print header */ printf(" ---------"); printf(" \n| SPTFQMR |\n"); printf(" ---------\n"); + if (monitor) { + fprintf(infofp, " ---------"); + fprintf(infofp, " \n| SPTFQMR |\n"); + fprintf(infofp, " ---------\n"); + } /* Call SUNLinSol_SPTFQMR to specify the linear solver SPTFQMR with left preconditioning and the default maximum Krylov dimension */ LS = SUNLinSol_SPTFQMR(u, PREC_LEFT, 0); - if(check_retval((void *)LS, "SUNLinSol_SPTFQMR", 0)) return(1); - + if (check_retval((void *)LS, "SUNLinSol_SPTFQMR", 0)) return(1); + if (monitor) { + retval = SUNLinSolSetPrintLevel_SPTFQMR(LS, 1); + if (check_retval(&retval, "SUNLinSolSetPrintLevel_SPTFQMR", 1)) return(1); + retval = SUNLinSolSetInfoFile_SPTFQMR(LS, infofp); + if (check_retval(&retval, "SUNLinSolSetInfoFile_SPTFQMR", 1)) return(1); + } retval = CVodeSetLinearSolver(cvode_mem, LS, NULL); - if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return 1; + if (check_retval(&retval, "CVodeSetLinearSolver", 1)) return 1; break; - } /* Set preconditioner setup and solve routines Precond and PSolve, and the pointer to the user-defined block data */ retval = CVodeSetPreconditioner(cvode_mem, Precond, PSolve); - if(check_retval(&retval, "CVodeSetPreconditioner", 1)) return(1); + if (check_retval(&retval, "CVodeSetPreconditioner", 1)) return(1); - /* In loop over output points, call CVode, print results, test for error */ + /* In loop over output points, call CVode, print results, and test for error */ printf(" \n2-species diurnal advection-diffusion problem\n\n"); for (iout=1, tout = TWOHR; iout <= NOUT; iout++, tout += TWOHR) { retval = CVode(cvode_mem, tout, u, &t, CV_NORMAL); - PrintOutput(cvode_mem, u, t); - if(check_retval(&retval, "CVode", 1)) break; + if (!monitor) PrintOutput(cvode_mem, u, t); + if (check_retval(&retval, "CVode", 1)) break; } - - PrintFinalStats(cvode_mem, linsolver); + if (monitor) PrintOutput(cvode_mem, u, t); + PrintStats(cvode_mem, linsolver, 1); } /* END: Loop through SPGMR, SPBCG and SPTFQMR linear solver modules */ /* Free memory */ + if (monitor) fclose(infofp); N_VDestroy(u); FreeUserData(data); CVodeFree(&cvode_mem); @@ -345,8 +452,9 @@ static UserData AllocUserData(void) /* Load problem constants in data */ -static void InitUserData(UserData data) +static void InitUserData(UserData data, N_Vector u) { + data->u = u; data->om = PI/HALFDAY; data->dx = (XMAX-XMIN)/(MX-1); data->dy = (YMAX-YMIN)/(MY-1); @@ -444,12 +552,12 @@ static void PrintOutput(void *cvode_mem, N_Vector u, realtype t) /* Get and print final statistics */ -static void PrintFinalStats(void *cvode_mem, int linsolver) +static void PrintStats(void *cvode_mem, int linsolver, int final) { long int lenrw, leniw ; long int lenrwLS, leniwLS; long int nst, nfe, nsetups, nni, ncfn, netf; - long int nli, npe, nps, ncfl, nfeLS; + long int nje, nli, npe, nps, ncfl, nfeLS, njts, njte; int retval; retval = CVodeGetWorkSpace(cvode_mem, &lenrw, &leniw); @@ -469,18 +577,12 @@ static void PrintFinalStats(void *cvode_mem, int linsolver) retval = CVodeGetLinWorkSpace(cvode_mem, &lenrwLS, &leniwLS); check_retval(&retval, "CVodeGetLinWorkSpace", 1); - retval = CVodeGetNumLinIters(cvode_mem, &nli); - check_retval(&retval, "CVodeGetNumLinIters", 1); - retval = CVodeGetNumPrecEvals(cvode_mem, &npe); - check_retval(&retval, "CVodeGetNumPrecEvals", 1); - retval = CVodeGetNumPrecSolves(cvode_mem, &nps); - check_retval(&retval, "CVodeGetNumPrecSolves", 1); - retval = CVodeGetNumLinConvFails(cvode_mem, &ncfl); - check_retval(&retval, "CVodeGetNumLinConvFails", 1); - retval = CVodeGetNumLinRhsEvals(cvode_mem, &nfeLS); - check_retval(&retval, "CVodeGetNumLinRhsEvals", 1); - - printf("\nFinal Statistics.. \n\n"); + CVodeGetLinSolveStats(cvode_mem, &nje, &nfeLS, &nli, &ncfl, &npe, + &nps, &njts, &njte); + check_retval(&retval, "CVodeGetLinWorkSpace", 1); + + if (final) printf("\nFinal Statistics.. \n\n"); + else printf("\nIntermediate Statistics.. \n\n"); printf("lenrw = %5ld leniw = %5ld\n" , lenrw, leniw); printf("lenrwLS = %5ld leniwLS = %5ld\n" , lenrwLS, leniwLS); printf("nst = %5ld\n" , nst); @@ -748,3 +850,18 @@ static int PSolve(realtype tn, N_Vector u, N_Vector fu, return(0); } + + +/* Function that is called at some step interval by CVODE */ + +static int myMonitorFunction(void* cvode_mem, void* user_data) +{ + UserData data = (UserData) user_data; + realtype t = 0; + + CVodeGetCurrentTime(cvode_mem, &t); + PrintOutput(cvode_mem, data->u, t); + PrintStats(cvode_mem, data->linsolver, 0); + + return(0); +} diff --git a/examples/cvode/serial/cvKrylovDemo_ls.out b/examples/cvode/serial/cvKrylovDemo_ls.out index 522a87d9ad..6f1196281d 100644 --- a/examples/cvode/serial/cvKrylovDemo_ls.out +++ b/examples/cvode/serial/cvKrylovDemo_ls.out @@ -64,6 +64,74 @@ nsetups = 78 netf = 27 npe = 8 nps = 1176 ncfn = 0 ncfl = 0 +====================================================================== + + --------- +| SPFGMR | + --------- + +2-species diurnal advection-diffusion problem + +t = 7.20e+03 no. steps = 192 order = 5 stepsize = 1.34e+02 +c1 (bot.left/middle/top rt.) = 1.047e+04 2.964e+04 1.119e+04 +c2 (bot.left/middle/top rt.) = 2.527e+11 7.154e+11 2.700e+11 + +t = 1.44e+04 no. steps = 225 order = 5 stepsize = 3.12e+02 +c1 (bot.left/middle/top rt.) = 6.659e+06 5.316e+06 7.301e+06 +c2 (bot.left/middle/top rt.) = 2.582e+11 2.057e+11 2.833e+11 + +t = 2.16e+04 no. steps = 248 order = 5 stepsize = 3.12e+02 +c1 (bot.left/middle/top rt.) = 2.665e+07 1.036e+07 2.931e+07 +c2 (bot.left/middle/top rt.) = 2.993e+11 1.028e+11 3.313e+11 + +t = 2.88e+04 no. steps = 293 order = 3 stepsize = 9.18e+01 +c1 (bot.left/middle/top rt.) = 8.702e+06 1.292e+07 9.650e+06 +c2 (bot.left/middle/top rt.) = 3.380e+11 5.029e+11 3.751e+11 + +t = 3.60e+04 no. steps = 324 order = 5 stepsize = 1.15e+02 +c1 (bot.left/middle/top rt.) = 1.404e+04 2.029e+04 1.561e+04 +c2 (bot.left/middle/top rt.) = 3.387e+11 4.895e+11 3.765e+11 + +t = 4.32e+04 no. steps = 378 order = 5 stepsize = 5.86e+02 +c1 (bot.left/middle/top rt.) = 1.180e-09 7.956e-08 1.490e-09 +c2 (bot.left/middle/top rt.) = 3.382e+11 1.355e+11 3.804e+11 + +t = 5.04e+04 no. steps = 391 order = 5 stepsize = 4.27e+02 +c1 (bot.left/middle/top rt.) = 6.580e-11 2.809e-08 1.507e-10 +c2 (bot.left/middle/top rt.) = 3.358e+11 4.930e+11 3.864e+11 + +t = 5.76e+04 no. steps = 402 order = 5 stepsize = 4.96e+02 +c1 (bot.left/middle/top rt.) = 2.069e-10 1.559e-08 2.264e-10 +c2 (bot.left/middle/top rt.) = 3.320e+11 9.650e+11 3.909e+11 + +t = 6.48e+04 no. steps = 412 order = 5 stepsize = 7.86e+02 +c1 (bot.left/middle/top rt.) = 3.735e-11 -1.076e-09 -3.688e-11 +c2 (bot.left/middle/top rt.) = 3.313e+11 8.922e+11 3.963e+11 + +t = 7.20e+04 no. steps = 421 order = 5 stepsize = 7.86e+02 +c1 (bot.left/middle/top rt.) = -1.410e-11 -1.148e-09 9.417e-11 +c2 (bot.left/middle/top rt.) = 3.330e+11 6.186e+11 4.039e+11 + +t = 7.92e+04 no. steps = 430 order = 5 stepsize = 7.86e+02 +c1 (bot.left/middle/top rt.) = 1.671e-12 -2.171e-10 -4.471e-12 +c2 (bot.left/middle/top rt.) = 3.334e+11 6.669e+11 4.120e+11 + +t = 8.64e+04 no. steps = 439 order = 5 stepsize = 7.86e+02 +c1 (bot.left/middle/top rt.) = 5.687e-13 -9.924e-11 -1.149e-12 +c2 (bot.left/middle/top rt.) = 3.352e+11 9.106e+11 4.162e+11 + + +Final Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 3254 leniwLS = 46 +nst = 439 +nfe = 565 nfeLS = 858 +nni = 562 nli = 858 +nsetups = 77 netf = 26 +npe = 8 nps = 858 +ncfn = 0 ncfl = 0 + ====================================================================== ------- @@ -132,8 +200,6 @@ nsetups = 72 netf = 25 npe = 8 nps = 1453 ncfn = 0 ncfl = 0 -====================================================================== - --------- | SPTFQMR | --------- diff --git a/examples/cvode/serial/cvKrylovDemo_ls_1.out b/examples/cvode/serial/cvKrylovDemo_ls_1.out new file mode 100644 index 0000000000..2a61824c44 --- /dev/null +++ b/examples/cvode/serial/cvKrylovDemo_ls_1.out @@ -0,0 +1,654 @@ + ------- +| SPGMR | + ------- + +2-species diurnal advection-diffusion problem + +t = 7.88e-01 no. steps = 50 order = 5 stepsize = 2.51e-02 +c1 (bot.left/middle/top rt.) = 2.151e+03 8.392e+03 2.151e+03 +c2 (bot.left/middle/top rt.) = 2.500e+11 9.756e+11 2.500e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2454 leniwLS = 42 +nst = 50 +nfe = 80 nfeLS = 55 +nni = 77 nli = 55 +nsetups = 16 netf = 6 +npe = 1 nps = 110 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 2.21e+00 no. steps = 100 order = 5 stepsize = 3.77e-02 +c1 (bot.left/middle/top rt.) = 4.149e-01 1.618e+00 4.149e-01 +c2 (bot.left/middle/top rt.) = 2.500e+11 9.757e+11 2.500e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2454 leniwLS = 42 +nst = 100 +nfe = 131 nfeLS = 105 +nni = 128 nli = 105 +nsetups = 19 netf = 6 +npe = 2 nps = 210 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 3.67e+03 no. steps = 150 order = 4 stepsize = 3.74e+01 +c1 (bot.left/middle/top rt.) = 1.236e-02 4.596e-02 1.288e-02 +c2 (bot.left/middle/top rt.) = 2.509e+11 9.328e+11 2.615e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2454 leniwLS = 42 +nst = 150 +nfe = 197 nfeLS = 162 +nni = 194 nli = 162 +nsetups = 34 netf = 8 +npe = 3 nps = 321 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 5.77e+03 no. steps = 200 order = 3 stepsize = 3.60e+01 +c1 (bot.left/middle/top rt.) = 3.285e+02 1.066e+03 3.480e+02 +c2 (bot.left/middle/top rt.) = 2.518e+11 8.171e+11 2.667e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2454 leniwLS = 42 +nst = 200 +nfe = 259 nfeLS = 214 +nni = 256 nli = 214 +nsetups = 38 netf = 10 +npe = 4 nps = 425 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 1.43e+04 no. steps = 250 order = 5 stepsize = 3.77e+02 +c1 (bot.left/middle/top rt.) = 6.387e+06 5.244e+06 7.000e+06 +c2 (bot.left/middle/top rt.) = 2.581e+11 2.116e+11 2.831e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2454 leniwLS = 42 +nst = 250 +nfe = 320 nfeLS = 265 +nni = 317 nli = 265 +nsetups = 45 netf = 11 +npe = 5 nps = 535 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 2.79e+04 no. steps = 300 order = 4 stepsize = 4.41e+01 +c1 (bot.left/middle/top rt.) = 1.189e+07 1.588e+07 1.317e+07 +c2 (bot.left/middle/top rt.) = 3.375e+11 4.520e+11 3.743e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2454 leniwLS = 42 +nst = 300 +nfe = 391 nfeLS = 337 +nni = 388 nli = 337 +nsetups = 54 netf = 17 +npe = 6 nps = 676 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 3.74e+04 no. steps = 350 order = 5 stepsize = 1.02e+02 +c1 (bot.left/middle/top rt.) = 4.321e+02 5.131e+02 4.808e+02 +c2 (bot.left/middle/top rt.) = 3.387e+11 4.022e+11 3.769e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2454 leniwLS = 42 +nst = 350 +nfe = 461 nfeLS = 428 +nni = 458 nli = 428 +nsetups = 62 netf = 22 +npe = 6 nps = 837 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 4.98e+04 no. steps = 400 order = 4 stepsize = 3.15e+02 +c1 (bot.left/middle/top rt.) = -3.654e-08 -4.295e-06 -6.514e-08 +c2 (bot.left/middle/top rt.) = 3.361e+11 4.418e+11 3.860e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2454 leniwLS = 42 +nst = 400 +nfe = 518 nfeLS = 521 +nni = 515 nli = 521 +nsetups = 68 netf = 23 +npe = 7 nps = 985 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 7.88e+04 no. steps = 450 order = 5 stepsize = 7.01e+02 +c1 (bot.left/middle/top rt.) = 6.689e-13 3.345e-11 1.269e-12 +c2 (bot.left/middle/top rt.) = 3.334e+11 6.542e+11 4.117e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2454 leniwLS = 42 +nst = 450 +nfe = 586 nfeLS = 614 +nni = 583 nli = 614 +nsetups = 77 netf = 27 +npe = 8 nps = 1143 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 8.64e+04 no. steps = 461 order = 5 stepsize = 7.01e+02 +c1 (bot.left/middle/top rt.) = 4.144e-14 2.072e-12 7.862e-14 +c2 (bot.left/middle/top rt.) = 3.352e+11 9.107e+11 4.163e+11 + + +Final Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2454 leniwLS = 42 +nst = 461 +nfe = 597 nfeLS = 636 +nni = 594 nli = 636 +nsetups = 78 netf = 27 +npe = 8 nps = 1176 +ncfn = 0 ncfl = 0 + +====================================================================== + + --------- +| SPFGMR | + --------- + +2-species diurnal advection-diffusion problem + +t = 7.88e-01 no. steps = 50 order = 5 stepsize = 2.51e-02 +c1 (bot.left/middle/top rt.) = 2.151e+03 8.392e+03 2.151e+03 +c2 (bot.left/middle/top rt.) = 2.500e+11 9.756e+11 2.500e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 3254 leniwLS = 46 +nst = 50 +nfe = 80 nfeLS = 55 +nni = 77 nli = 55 +nsetups = 16 netf = 6 +npe = 1 nps = 55 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 2.21e+00 no. steps = 100 order = 5 stepsize = 3.77e-02 +c1 (bot.left/middle/top rt.) = 4.149e-01 1.618e+00 4.149e-01 +c2 (bot.left/middle/top rt.) = 2.500e+11 9.757e+11 2.500e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 3254 leniwLS = 46 +nst = 100 +nfe = 131 nfeLS = 105 +nni = 128 nli = 105 +nsetups = 19 netf = 6 +npe = 2 nps = 105 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 3.94e+03 no. steps = 150 order = 4 stepsize = 9.97e+01 +c1 (bot.left/middle/top rt.) = 8.547e-02 3.136e-01 8.929e-02 +c2 (bot.left/middle/top rt.) = 2.510e+11 9.210e+11 2.622e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 3254 leniwLS = 46 +nst = 150 +nfe = 200 nfeLS = 171 +nni = 197 nli = 171 +nsetups = 37 netf = 9 +npe = 3 nps = 171 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 8.36e+03 no. steps = 200 order = 5 stepsize = 1.34e+02 +c1 (bot.left/middle/top rt.) = 6.996e+04 1.728e+05 7.519e+04 +c2 (bot.left/middle/top rt.) = 2.535e+11 6.262e+11 2.724e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 3254 leniwLS = 46 +nst = 200 +nfe = 261 nfeLS = 225 +nni = 258 nli = 225 +nsetups = 43 netf = 12 +npe = 4 nps = 225 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 2.23e+04 no. steps = 250 order = 5 stepsize = 3.12e+02 +c1 (bot.left/middle/top rt.) = 2.698e+07 1.250e+07 2.967e+07 +c2 (bot.left/middle/top rt.) = 3.075e+11 1.307e+11 3.403e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 3254 leniwLS = 46 +nst = 250 +nfe = 315 nfeLS = 359 +nni = 312 nli = 359 +nsetups = 46 netf = 12 +npe = 5 nps = 359 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 3.01e+04 no. steps = 300 order = 4 stepsize = 1.75e+02 +c1 (bot.left/middle/top rt.) = 5.017e+06 8.298e+06 5.570e+06 +c2 (bot.left/middle/top rt.) = 3.383e+11 5.600e+11 3.757e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 3254 leniwLS = 46 +nst = 300 +nfe = 396 nfeLS = 484 +nni = 393 nli = 484 +nsetups = 62 netf = 20 +npe = 6 nps = 484 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 3.86e+04 no. steps = 350 order = 5 stepsize = 8.38e+01 +c1 (bot.left/middle/top rt.) = 4.314e+00 4.141e+00 4.806e+00 +c2 (bot.left/middle/top rt.) = 3.387e+11 3.251e+11 3.774e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 3254 leniwLS = 46 +nst = 350 +nfe = 462 nfeLS = 656 +nni = 459 nli = 656 +nsetups = 67 netf = 24 +npe = 6 nps = 656 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 5.68e+04 no. steps = 400 order = 5 stepsize = 7.44e+02 +c1 (bot.left/middle/top rt.) = -1.055e-12 -8.927e-09 2.668e-11 +c2 (bot.left/middle/top rt.) = 3.324e+11 9.373e+11 3.904e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 3254 leniwLS = 46 +nst = 400 +nfe = 519 nfeLS = 758 +nni = 516 nli = 758 +nsetups = 74 netf = 25 +npe = 7 nps = 758 +ncfn = 0 ncfl = 0 + +====================================================================== + +t = 8.64e+04 no. steps = 439 order = 5 stepsize = 7.86e+02 +c1 (bot.left/middle/top rt.) = 5.687e-13 -9.924e-11 -1.149e-12 +c2 (bot.left/middle/top rt.) = 3.352e+11 9.106e+11 4.162e+11 + + +Final Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 3254 leniwLS = 46 +nst = 439 +nfe = 565 nfeLS = 858 +nni = 562 nli = 858 +nsetups = 77 netf = 26 +npe = 8 nps = 858 +ncfn = 0 ncfl = 0 + +====================================================================== + + ------- +| SPBCGS | + ------- + +2-species diurnal advection-diffusion problem + +t = 7.88e-01 no. steps = 50 order = 5 stepsize = 2.51e-02 +c1 (bot.left/middle/top rt.) = 2.151e+03 8.392e+03 2.151e+03 +c2 (bot.left/middle/top rt.) = 2.500e+11 9.756e+11 2.500e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2202 leniwLS = 41 +nst = 50 +nfe = 80 nfeLS = 110 +nni = 77 nli = 55 +nsetups = 16 netf = 6 +npe = 1 nps = 165 +ncfn = 0 ncfl = 0 + +t = 2.21e+00 no. steps = 100 order = 5 stepsize = 3.77e-02 +c1 (bot.left/middle/top rt.) = 4.149e-01 1.618e+00 4.149e-01 +c2 (bot.left/middle/top rt.) = 2.500e+11 9.757e+11 2.500e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2202 leniwLS = 41 +nst = 100 +nfe = 131 nfeLS = 210 +nni = 128 nli = 105 +nsetups = 19 netf = 6 +npe = 2 nps = 315 +ncfn = 0 ncfl = 0 + +t = 3.88e+03 no. steps = 150 order = 4 stepsize = 1.45e+02 +c1 (bot.left/middle/top rt.) = 5.801e-02 2.135e-01 6.057e-02 +c2 (bot.left/middle/top rt.) = 2.509e+11 9.235e+11 2.620e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2202 leniwLS = 41 +nst = 150 +nfe = 197 nfeLS = 312 +nni = 194 nli = 156 +nsetups = 34 netf = 8 +npe = 3 nps = 470 +ncfn = 0 ncfl = 0 + +t = 8.84e+03 no. steps = 200 order = 5 stepsize = 1.58e+02 +c1 (bot.left/middle/top rt.) = 1.311e+05 3.042e+05 1.413e+05 +c2 (bot.left/middle/top rt.) = 2.538e+11 5.888e+11 2.734e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2202 leniwLS = 41 +nst = 200 +nfe = 261 nfeLS = 418 +nni = 258 nli = 209 +nsetups = 41 netf = 11 +npe = 4 nps = 629 +ncfn = 0 ncfl = 0 + +t = 2.36e+04 no. steps = 250 order = 5 stepsize = 4.34e+02 +c1 (bot.left/middle/top rt.) = 2.596e+07 1.592e+07 2.856e+07 +c2 (bot.left/middle/top rt.) = 3.194e+11 1.885e+11 3.533e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2202 leniwLS = 41 +nst = 250 +nfe = 318 nfeLS = 520 +nni = 315 nli = 260 +nsetups = 46 netf = 12 +npe = 5 nps = 788 +ncfn = 0 ncfl = 0 + +t = 3.34e+04 no. steps = 300 order = 5 stepsize = 2.67e+02 +c1 (bot.left/middle/top rt.) = 5.232e+05 9.140e+05 5.813e+05 +c2 (bot.left/middle/top rt.) = 3.385e+11 5.914e+11 3.761e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2202 leniwLS = 41 +nst = 300 +nfe = 395 nfeLS = 642 +nni = 392 nli = 321 +nsetups = 58 netf = 20 +npe = 5 nps = 975 +ncfn = 0 ncfl = 0 + +t = 3.92e+04 no. steps = 350 order = 5 stepsize = 9.48e+01 +c1 (bot.left/middle/top rt.) = 1.804e-01 1.542e-01 2.011e-01 +c2 (bot.left/middle/top rt.) = 3.387e+11 2.896e+11 3.776e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2202 leniwLS = 41 +nst = 350 +nfe = 460 nfeLS = 752 +nni = 457 nli = 376 +nsetups = 63 netf = 23 +npe = 6 nps = 1140 +ncfn = 0 ncfl = 0 + +t = 5.86e+04 no. steps = 400 order = 5 stepsize = 5.96e+02 +c1 (bot.left/middle/top rt.) = 7.894e-14 -5.680e-12 1.260e-13 +c2 (bot.left/middle/top rt.) = 3.317e+11 9.872e+11 3.915e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2202 leniwLS = 41 +nst = 400 +nfe = 522 nfeLS = 870 +nni = 519 nli = 435 +nsetups = 70 netf = 25 +npe = 7 nps = 1312 +ncfn = 0 ncfl = 0 + +t = 8.64e+04 no. steps = 447 order = 5 stepsize = 5.96e+02 +c1 (bot.left/middle/top rt.) = -9.999e-27 -1.866e-22 -1.581e-27 +c2 (bot.left/middle/top rt.) = 3.352e+11 9.107e+11 4.163e+11 + + +Final Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2202 leniwLS = 41 +nst = 447 +nfe = 569 nfeLS = 964 +nni = 566 nli = 482 +nsetups = 72 netf = 25 +npe = 8 nps = 1453 +ncfn = 0 ncfl = 0 + + --------- +| SPTFQMR | + --------- + +2-species diurnal advection-diffusion problem + +t = 7.88e-01 no. steps = 50 order = 5 stepsize = 2.51e-02 +c1 (bot.left/middle/top rt.) = 2.151e+03 8.392e+03 2.151e+03 +c2 (bot.left/middle/top rt.) = 2.500e+11 9.756e+11 2.500e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2602 leniwLS = 43 +nst = 50 +nfe = 80 nfeLS = 110 +nni = 77 nli = 55 +nsetups = 16 netf = 6 +npe = 1 nps = 165 +ncfn = 0 ncfl = 0 + +t = 2.21e+00 no. steps = 100 order = 5 stepsize = 3.77e-02 +c1 (bot.left/middle/top rt.) = 4.149e-01 1.618e+00 4.149e-01 +c2 (bot.left/middle/top rt.) = 2.500e+11 9.757e+11 2.500e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2602 leniwLS = 43 +nst = 100 +nfe = 131 nfeLS = 210 +nni = 128 nli = 105 +nsetups = 19 netf = 6 +npe = 2 nps = 315 +ncfn = 0 ncfl = 0 + +t = 3.65e+03 no. steps = 150 order = 4 stepsize = 3.82e+01 +c1 (bot.left/middle/top rt.) = 1.067e-02 3.969e-02 1.111e-02 +c2 (bot.left/middle/top rt.) = 2.509e+11 9.336e+11 2.614e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2602 leniwLS = 43 +nst = 150 +nfe = 197 nfeLS = 318 +nni = 194 nli = 156 +nsetups = 34 netf = 8 +npe = 3 nps = 483 +ncfn = 0 ncfl = 0 + +t = 5.70e+03 no. steps = 200 order = 3 stepsize = 3.36e+01 +c1 (bot.left/middle/top rt.) = 2.634e+02 8.596e+02 2.788e+02 +c2 (bot.left/middle/top rt.) = 2.518e+11 8.218e+11 2.666e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2602 leniwLS = 43 +nst = 200 +nfe = 259 nfeLS = 422 +nni = 256 nli = 208 +nsetups = 38 netf = 10 +npe = 4 nps = 639 +ncfn = 0 ncfl = 0 + +t = 1.45e+04 no. steps = 250 order = 5 stepsize = 3.27e+02 +c1 (bot.left/middle/top rt.) = 6.824e+06 5.357e+06 7.483e+06 +c2 (bot.left/middle/top rt.) = 2.583e+11 2.022e+11 2.834e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2602 leniwLS = 43 +nst = 250 +nfe = 318 nfeLS = 525 +nni = 315 nli = 259 +nsetups = 44 netf = 11 +npe = 5 nps = 798 +ncfn = 0 ncfl = 0 + +t = 2.94e+04 no. steps = 300 order = 5 stepsize = 3.39e+02 +c1 (bot.left/middle/top rt.) = 6.796e+06 1.069e+07 7.541e+06 +c2 (bot.left/middle/top rt.) = 3.382e+11 5.329e+11 3.754e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2602 leniwLS = 43 +nst = 300 +nfe = 382 nfeLS = 671 +nni = 379 nli = 313 +nsetups = 51 netf = 14 +npe = 6 nps = 1040 +ncfn = 0 ncfl = 0 + +t = 3.87e+04 no. steps = 350 order = 5 stepsize = 9.42e+01 +c1 (bot.left/middle/top rt.) = 3.613e+00 3.443e+00 4.025e+00 +c2 (bot.left/middle/top rt.) = 3.387e+11 3.228e+11 3.774e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2602 leniwLS = 43 +nst = 350 +nfe = 448 nfeLS = 840 +nni = 445 nli = 370 +nsetups = 55 netf = 18 +npe = 6 nps = 1324 +ncfn = 0 ncfl = 0 + +t = 5.67e+04 no. steps = 400 order = 5 stepsize = 6.79e+02 +c1 (bot.left/middle/top rt.) = 4.357e-09 5.640e-07 2.278e-09 +c2 (bot.left/middle/top rt.) = 3.324e+11 9.339e+11 3.904e+11 + + +Intermediate Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2602 leniwLS = 43 +nst = 400 +nfe = 505 nfeLS = 987 +nni = 502 nli = 424 +nsetups = 62 netf = 19 +npe = 7 nps = 1554 +ncfn = 0 ncfl = 0 + +t = 8.64e+04 no. steps = 445 order = 5 stepsize = 6.89e+02 +c1 (bot.left/middle/top rt.) = -1.673e-20 -1.128e-15 -3.081e-21 +c2 (bot.left/middle/top rt.) = 3.352e+11 9.106e+11 4.163e+11 + + +Final Statistics.. + +lenrw = 2689 leniw = 53 +lenrwLS = 2602 leniwLS = 43 +nst = 445 +nfe = 558 nfeLS = 1143 +nni = 555 nli = 477 +nsetups = 67 netf = 21 +npe = 8 nps = 1805 +ncfn = 0 ncfl = 0 + diff --git a/examples/cvode/serial/cvParticle_dns.c b/examples/cvode/serial/cvParticle_dns.c new file mode 100644 index 0000000000..a5d8a29659 --- /dev/null +++ b/examples/cvode/serial/cvParticle_dns.c @@ -0,0 +1,620 @@ +/* ----------------------------------------------------------------------------- + * Programmer(s): David J. Gardner @ LLNL + * ----------------------------------------------------------------------------- + * Based on an example from Jean-Luc Fattebert @ ORNL + * ----------------------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2002-2020, Lawrence Livermore National Security + * and Southern Methodist University. + * All rights reserved. + * + * See the top-level LICENSE and NOTICE files for details. + * + * SPDX-License-Identifier: BSD-3-Clause + * SUNDIALS Copyright End + * ----------------------------------------------------------------------------- + * This example solves the equation for a particle moving conterclockwise with + * velocity alpha on the unit circle in the xy-plane. The ODE system is given by + * + * x' = -alpha * y + * y' = alpha * x + * + * where x and y are subject to the constraint + * + * x^2 + y^2 - 1 = 0 + * + * with initial condition x = 1 and y = 0 at t = 0. The system has the analytic + * solution + * + * x(t) = cos(alpha * t) + * y(t) = sin(alpha * t) + * + * For a description of the command line options for this example run the + * program with the --help flag. + * ---------------------------------------------------------------------------*/ + +#include +#include +#include +#include + +#include /* access to CVODE */ +#include /* access to serial N_Vector */ +#include /* access to dense SUNMatrix */ +#include /* access to dense SUNLinearSolver */ + +/* Precision specific formatting macros */ +#if defined(SUNDIALS_EXTENDED_PRECISION) +#define GSYM "Lg" +#define ESYM "Le" +#define FSYM "Lf" +#else +#define GSYM "g" +#define ESYM "e" +#define FSYM "f" +#endif + +/* Precision specific math function macros */ +#if defined(SUNDIALS_DOUBLE_PRECISION) +#define SIN(x) (sin((x))) +#define COS(x) (cos((x))) +#define SQRT(x) (sqrt((x))) +#elif defined(SUNDIALS_SINGLE_PRECISION) +#define SIN(x) (sinf((x))) +#define COS(x) (cosf((x))) +#define SQRT(x) (sqrtf((x))) +#elif defined(SUNDIALS_EXTENDED_PRECISION) +#define SIN(x) (sinl((x))) +#define COS(x) (cosl((x))) +#define SQRT(x) (sqrtl((x))) +#endif + +/* Problem Constants */ +#define PI RCONST(3.141592653589793238462643383279502884197169) +#define ZERO RCONST(0.0) +#define ONE RCONST(1.0) +#define TWO RCONST(2.0) + +/* User-defined data structure */ +typedef struct UserData_ +{ + realtype alpha; /* particle velocity */ + + int orbits; /* number of orbits */ + realtype torbit; /* orbit time */ + + realtype rtol; /* integration tolerances */ + realtype atol; + + int proj; /* enable/disable solution projection */ + int projerr; /* enable/disable error projection */ + + int tstop; /* use tstop mode */ + int nout; /* number of outputs per orbit */ + +} *UserData; + +/* Functions provided to CVODE */ +static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data); +static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J, + void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); +static int Proj(realtype t, N_Vector ycur, N_Vector corr, realtype epsProj, + N_Vector err, void *user_data); + +/* Utility functions */ +static int InitUserData(int *argc, char ***argv, UserData udata); +static int PrintUserData(UserData udata); +static void InputHelp(); +static int ComputeSolution(realtype t, N_Vector y, UserData udata); +static int ComputeError(realtype t, N_Vector y, N_Vector e, realtype *ec, + UserData udata); +static int WriteOutput(realtype t, N_Vector y, N_Vector e, realtype ec, + int screenfile, FILE *YFID, FILE *EFID); +static int PrintStats(void *cvode_mem); +static int check_retval(void *returnvalue, const char *funcname, int opt); + + +/* ----------------------------------------------------------------------------- + * Main Program + * ---------------------------------------------------------------------------*/ + +int main(int argc, char* argv[]) +{ + int retval; /* reusable return flag */ + int out = 0; /* output counter */ + int totalout = 0; /* output counter */ + realtype t = ZERO; /* current integration time */ + realtype dtout = ZERO; /* output spacing */ + realtype tout = ZERO; /* next output time */ + realtype ec = ZERO; /* constraint error */ + UserData udata = NULL; /* user data structure */ + + void *cvode_mem = NULL; /* CVODE memory */ + N_Vector y = NULL; /* solution vector */ + realtype *ydata = NULL; /* solution vector data */ + N_Vector e = NULL; /* error vector */ + SUNMatrix A = NULL; /* Jacobian matrix */ + SUNLinearSolver LS = NULL; /* linear solver */ + + FILE *YFID = NULL; /* solution output file */ + FILE *EFID = NULL; /* error output file */ + + /* Allocate and initialize user data structure */ + udata = (UserData) malloc(sizeof *udata); + if (check_retval((void *)udata, "malloc", 0)) return(1); + + retval = InitUserData(&argc, &argv, udata); + if (check_retval(&retval, "InitUserData", 1)) return(1); + + /* Create serial vector to store the solution */ + y = N_VNew_Serial(2); + if (check_retval((void *)y, "N_VNew_Serial", 0)) return(1); + + /* Set initial contion */ + ydata = N_VGetArrayPointer(y); + ydata[0] = ONE; + ydata[1] = ZERO; + + /* Create serial vector to store the solution error */ + e = N_VClone(y); + if (check_retval((void *)y, "N_VClone", 0)) return(1); + + /* Set initial error */ + N_VConst(ZERO, e); + + /* Create CVODE memory */ + cvode_mem = CVodeCreate(CV_BDF); + if (check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1); + + /* Initialize CVODE */ + retval = CVodeInit(cvode_mem, f, t, y); + if (check_retval(&retval, "CVodeInit", 1)) return(1); + + /* Attach user-defined data structure to CVODE */ + retval = CVodeSetUserData(cvode_mem, udata); + if(check_retval(&retval, "CVodeSetUserData", 1)) return(1); + + /* Set integration tolerances */ + retval = CVodeSStolerances(cvode_mem, udata->rtol, udata->atol); + if (check_retval(&retval, "CVodeSStolerances", 1)) return(1); + + /* Create dense SUNMatrix for use in linear solves */ + A = SUNDenseMatrix(2, 2); + if(check_retval((void *)A, "SUNDenseMatrix", 0)) return(1); + + /* Create dense SUNLinearSolver object */ + LS = SUNLinSol_Dense(y, A); + if(check_retval((void *)LS, "SUNLinSol_Dense", 0)) return(1); + + /* Attach the matrix and linear solver to CVODE */ + retval = CVodeSetLinearSolver(cvode_mem, LS, A); + if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1); + + /* Set a user-supplied Jacobian function */ + retval = CVodeSetJacFn(cvode_mem, Jac); + if(check_retval(&retval, "CVodeSetJacFn", 1)) return(1); + + /* Set a user-supplied projection function */ + if (udata->proj) + { + retval = CVodeSetProjFn(cvode_mem, Proj); + if(check_retval(&retval, "CVodeSetProjFn", 1)) return(1); + + retval = CVodeSetProjErrEst(cvode_mem, udata->projerr); + if(check_retval(&retval, "CVodeSetProjErrEst", 1)) return(1); + } + + /* Set max steps between outputs */ + retval = CVodeSetMaxNumSteps(cvode_mem, 100000); + if (check_retval(&retval, "CVodeSetMaxNumSteps", 1)) return(1); + + /* Output problem setup */ + retval = PrintUserData(udata); + if(check_retval(&retval, "PrintUserData", 1)) return(1); + + /* Output initial condition */ + printf("\n t x y"); + printf(" err x err y err constr\n"); + WriteOutput(t, y, e, ec, 0, NULL, NULL); + + if (udata->nout > 0) + { + YFID = fopen("cvParticle_solution.txt","w"); + EFID = fopen("cvParticle_error.txt","w"); + WriteOutput(t, y, e, ec, 1, YFID, EFID); + } + + /* Integrate in time and periodically output the solution and error */ + if (udata->nout > 0) + { + totalout = udata->orbits * udata->nout; + dtout = udata->torbit / udata->nout; + } + else + { + totalout = 1; + dtout = udata->torbit * udata->orbits; + } + tout = dtout; + + for (out = 0; out < totalout; out++) + { + /* Stop at output time (do not interpolate output) */ + if (udata->tstop || udata->nout == 0) + { + retval = CVodeSetStopTime(cvode_mem, tout); + if (check_retval(&retval, "CVodeSetStopTime", 1)) return(1); + } + + /* Advance in time */ + retval = CVode(cvode_mem, tout, y, &t, CV_NORMAL); + if (check_retval(&retval, "CVode", 1)) break; + + /* Output solution and error */ + if (udata->nout > 0) + { + retval = ComputeError(t, y, e, &ec, udata); + if (check_retval(&retval, "ComputeError", 1)) break; + + WriteOutput(t, y, e, ec, 1, YFID, EFID); + if (check_retval(&retval, "WriteOutput", 1)) break; + } + + /* Update output time */ + if (out < totalout - 1) + { + tout += dtout; + } + else + { + tout = udata->torbit * udata->orbits; + } + } + + /* Close output files */ + if (udata->nout > 0) + { + fclose(YFID); + fclose(EFID); + } + + /* Output final solution and error to screen */ + ComputeError(t, y, e, &ec, udata); + if (check_retval(&retval, "ComputeError", 1)) return(1); + + WriteOutput(t, y, e, ec, 0, NULL, NULL); + if (check_retval(&retval, "WriteOutput", 1)) return(1); + + /* Print some final statistics */ + PrintStats(cvode_mem); + + /* Free memory */ + N_VDestroy(y); + SUNMatDestroy(A); + SUNLinSolFree(LS); + CVodeFree(&cvode_mem); + + return(0); +} + + +/* ----------------------------------------------------------------------------- + * Functions provided to CVODE + * ---------------------------------------------------------------------------*/ + + +/* Compute the right-hand side function, y' = f(t,y) */ +static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data) +{ + UserData udata = (UserData) user_data; + realtype *ydata = N_VGetArrayPointer(y); + realtype *fdata = N_VGetArrayPointer(ydot); + + fdata[0] = -(udata->alpha) * ydata[1]; + fdata[1] = (udata->alpha) * ydata[0]; + + return(0); +} + + +/* Compute the Jacobian of the right-hand side function, J(t,y) = df/dy */ +static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J, + void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) +{ + UserData udata = (UserData) user_data; + realtype *Jdata = SUNDenseMatrix_Data(J); + + Jdata[0] = ZERO; + Jdata[1] = -(udata->alpha); + Jdata[2] = (udata->alpha); + Jdata[3] = ZERO; + + return(0); +} + +/* Project the solution onto the constraint manifold */ +static int Proj(realtype t, N_Vector ycur, N_Vector corr, realtype epsProj, + N_Vector err, void *user_data) +{ + realtype *ydata = N_VGetArrayPointer(ycur); + realtype *cdata = N_VGetArrayPointer(corr); + realtype *edata = NULL; + realtype x = ydata[0]; + realtype y = ydata[1]; + realtype xp, yp, r; + realtype errxp, erryp; + + /* project onto the unit circle */ + r = SQRT(x * x + y * y); + + xp = x / r; + yp = y / r; + + /* correction to the unprojected solution */ + cdata[0] = xp - x; + cdata[1] = yp - y; + + /* project the error */ + if (err != NULL) + { + edata = N_VGetArrayPointer(err); + + errxp = edata[0] * yp * yp - edata[1] * xp * yp; + erryp = -edata[0] * xp * yp + edata[1] * xp * xp; + + edata[0] = errxp; + edata[1] = erryp; + } + + return(0); +} + + +/* ----------------------------------------------------------------------------- + * Private helper functions + * ---------------------------------------------------------------------------*/ + +static int InitUserData(int *argc, char ***argv, UserData udata) +{ + int arg_idx = 1; + + /* set default values */ + udata->alpha = ONE; + + udata->orbits = 100; + udata->torbit = (TWO * PI) / udata->alpha; + + udata->rtol = RCONST(1.0e-4); + udata->atol = RCONST(1.0e-9); + + udata->proj = 1; + udata->projerr = 0; + + udata->tstop = 0; + udata->nout = 0; + + /* check for input args */ + while (arg_idx < (*argc)) + { + if (strcmp((*argv)[arg_idx],"--alpha") == 0) + { + arg_idx++; + udata->alpha = atof((*argv)[arg_idx++]); + udata->torbit = (TWO * PI) / udata->alpha; + } + else if (strcmp((*argv)[arg_idx],"--orbits") == 0) + { + arg_idx++; + udata->orbits = atoi((*argv)[arg_idx++]); + } + else if (strcmp((*argv)[arg_idx],"--rtol") == 0) + { + arg_idx++; + udata->rtol = atof((*argv)[arg_idx++]); + } + else if (strcmp((*argv)[arg_idx],"--atol") == 0) + { + arg_idx++; + udata->atol = atof((*argv)[arg_idx++]); + } + else if (strcmp((*argv)[arg_idx],"--proj") == 0) + { + arg_idx++; + udata->proj = atoi((*argv)[arg_idx++]); + } + else if (strcmp((*argv)[arg_idx],"--projerr") == 0) + { + arg_idx++; + udata->projerr = atoi((*argv)[arg_idx++]); + } + else if (strcmp((*argv)[arg_idx],"--nout") == 0) + { + arg_idx++; + udata->nout = atoi((*argv)[arg_idx++]); + } + else if (strcmp((*argv)[arg_idx],"--tstop") == 0) + { + arg_idx++; + udata->tstop = 1; + } + else if (strcmp((*argv)[arg_idx],"--help") == 0 ) + { + InputHelp(); + return(-1); + } + else + { + fprintf(stderr, "ERROR: Invalid input %s",(*argv)[arg_idx]); + InputHelp(); + return(-1); + } + } + + /* If projection is disabled then disable error projection */ + if (!(udata->proj)) udata->projerr = 0; + + return(0); +} + +static int PrintUserData(UserData udata) +{ + if (udata == NULL) return(-1); + + printf("\nParticle traveling on the unit circle example\n"); + printf("---------------------------------------------\n"); + printf("alpha = %0.4" ESYM"\n", udata->alpha); + printf("num orbits = %d\n", udata->orbits); + printf("---------------------------------------------\n"); + printf("rtol = %" GSYM"\n", udata->rtol); + printf("atol = %" GSYM"\n", udata->atol); + printf("proj sol = %d\n", udata->proj); + printf("proj err = %d\n", udata->projerr); + printf("nout = %d\n", udata->nout); + printf("tstop = %d\n", udata->tstop); + printf("---------------------------------------------\n"); + + return(0); +} + + +/* Print command line options */ +static void InputHelp() +{ + printf("\nCommand line options:\n"); + printf(" --alpha : particle velocity\n"); + printf(" --orbits : number of orbits to perform\n"); + printf(" --rtol : relative tolerance\n"); + printf(" --atol : absoltue tolerance\n"); + printf(" --proj <1 or 0> : enable (1) / disable (0) projection\n"); + printf(" --projerr <1 or 0> : enable (1) / disable (0) error projection\n"); + printf(" --nout : outputs per period\n"); + printf(" --tstop : stop at output time (do not interpolate)\n"); + return; +} + + +/* Compute the analytical solution */ +static int ComputeSolution(realtype t, N_Vector y, UserData udata) +{ + realtype *ydata = N_VGetArrayPointer(y); + + ydata[0] = COS((udata->alpha) * t); + ydata[1] = SIN((udata->alpha) * t); + + return(0); +} + + +/* Compute the error in the solution and constraint */ +static int ComputeError(realtype t, N_Vector y, N_Vector e, realtype *ec, + UserData udata) +{ + realtype *ydata = N_VGetArrayPointer(y); + int retval; + + /* solution error */ + retval = ComputeSolution(t, e, udata); + if (check_retval(&retval, "ComputeSolution", 1)) return(1); + N_VLinearSum(ONE, y, -ONE, e, e); + + /* constraint error */ + *ec = ydata[0] * ydata[0] + ydata[1] * ydata[1] - ONE; + + return(0); +} + +/* Output the solution to the screen or disk */ +static int WriteOutput(realtype t, N_Vector y, N_Vector e, realtype ec, + int screenfile, FILE* YFID, FILE* EFID) +{ + realtype *ydata = N_VGetArrayPointer(y); + realtype *edata = N_VGetArrayPointer(e); + + if (screenfile == 0) + { + /* output solution and error to screen */ + printf("%0.4" ESYM" %14.6" ESYM" %14.6" ESYM" %14.6" ESYM" %14.6" ESYM" %14.6" ESYM"\n", + t, ydata[0], ydata[1], edata[0], edata[1], ec); + } + else + { + /* check file pointers */ + if (YFID == NULL || EFID == NULL) return(1); + + /* output solution to disk */ + fprintf(YFID, "%24.16" ESYM" %24.16" ESYM" %24.16"ESYM"\n", + t, ydata[0], ydata[1]); + + /* output error to disk */ + fprintf(EFID, + "%24.16" ESYM" %24.16" ESYM" %24.16"ESYM" %24.16"ESYM"\n", + t, edata[0], edata[1], ec); + } + + return(0); +} + + +/* Print final statistics */ +static int PrintStats(void *cvode_mem) +{ + int retval; + long int nst, nfe, nsetups, nje, nni, ncfn, netf; + + retval = CVodeGetNumSteps(cvode_mem, &nst); + check_retval(&retval, "CVodeGetNumSteps", 1); + retval = CVodeGetNumRhsEvals(cvode_mem, &nfe); + check_retval(&retval, "CVodeGetNumRhsEvals", 1); + retval = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups); + check_retval(&retval, "CVodeGetNumLinSolvSetups", 1); + retval = CVodeGetNumErrTestFails(cvode_mem, &netf); + check_retval(&retval, "CVodeGetNumErrTestFails", 1); + retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni); + check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1); + retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn); + check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1); + + retval = CVodeGetNumJacEvals(cvode_mem, &nje); + check_retval(&retval, "CVodeGetNumJacEvals", 1); + + printf("\nIntegration Statistics:\n"); + + printf("Number of steps taken = %-6ld\n", nst); + printf("Number of function evaluations = %-6ld\n", nfe); + + printf("Number of linear solver setups = %-6ld\n", nsetups); + printf("Number of Jacobian evaluations = %-6ld\n", nje); + + printf("Number of nonlinear solver iterations = %-6ld\n", nni); + printf("Number of convergence failures = %-6ld\n", ncfn); + printf("Number of error test failures = %-6ld\n", netf); + + return(0); +} + +/* Check function return value */ +static int check_retval(void *returnvalue, const char *funcname, int opt) +{ + int *retval; + + /* Opt 0: Check if a NULL pointer was returned - no memory allocated */ + if (opt == 0 && returnvalue == NULL) + { + fprintf(stderr, "\nERROR: %s() returned a NULL pointer\n\n", + funcname); + return(1); + } + /* Opt 1: Check if retval < 0 */ + else if (opt == 1) + { + retval = (int *) returnvalue; + if (*retval < 0) + { + fprintf(stderr, "\nERROR: %s() returned = %d\n\n", + funcname, *retval); + return(1); + } + } + + return(0); +} diff --git a/examples/cvode/serial/cvParticle_dns.out b/examples/cvode/serial/cvParticle_dns.out new file mode 100644 index 0000000000..a2639038fa --- /dev/null +++ b/examples/cvode/serial/cvParticle_dns.out @@ -0,0 +1,26 @@ + +Particle traveling on the unit circle example +--------------------------------------------- +alpha = 1.0000e+00 +num orbits = 100 +--------------------------------------------- +rtol = 0.0001 +atol = 1e-09 +proj sol = 1 +proj err = 0 +nout = 0 +tstop = 0 +--------------------------------------------- + + t x y err x err y err constr +0.0000e+00 1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 +6.2832e+02 9.977240e-01 6.742947e-02 -2.275957e-03 6.742947e-02 2.220446e-16 + +Integration Statistics: +Number of steps taken = 4044 +Number of function evaluations = 5831 +Number of linear solver setups = 806 +Number of Jacobian evaluations = 103 +Number of nonlinear solver iterations = 5828 +Number of convergence failures = 20 +Number of error test failures = 329 diff --git a/examples/cvode/serial/cvPendulum_dns.c b/examples/cvode/serial/cvPendulum_dns.c new file mode 100644 index 0000000000..800fede1b2 --- /dev/null +++ b/examples/cvode/serial/cvPendulum_dns.c @@ -0,0 +1,771 @@ +/* ----------------------------------------------------------------------------- + * Programmer(s): Radu Serban and David J. Gardner @ LLNL + * ----------------------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2002-2020, Lawrence Livermore National Security + * and Southern Methodist University. + * All rights reserved. + * + * See the top-level LICENSE and NOTICE files for details. + * + * SPDX-License-Identifier: BSD-3-Clause + * SUNDIALS Copyright End + * ----------------------------------------------------------------------------- + * This example solves a simple pendulum equation in Cartesian coordinates where + * the pendulum bob has mass 1 and is suspended from the origin with a rod of + * length 1. The governing equations are + * + * x' = vx + * y' = vy + * vx' = -x * T + * vy' = -y * T - g + * + * with the constraints + * + * x^2 + y^2 - 1 = 0 + * x * vx + y * vy = 0 + * + * where x and y are the pendulum bob position, vx and vy are the bob velocity + * in the x and y directions respectively, T is the tension in the rod, and + * g is acceleration due to gravity chosen such that the pendulum has period 2. + * The initial condition at t = 0 is x = 1, y = 0, vx = 0, and vy = 0. + * + * A reference solution is computed using the pendulum equation in terms of the + * angle between the x-axis and the pendulum rod i.e., theta in [0, -pi]. The + * governing equations are + * + * theta' = vtheta + * vtheta' = -g * cos(theta) + * + * where theta is the angle from the x-axis, vtheta is the angular velocity, and + * g the same acceleration due to gravity from above. The initial condition at + * t = 0 is theta = 0 and vtheta = 0. + * + * The Cartesian formulation is run to a final time tf (default 30) with and + * without projection for various integration tolerances. The error in the + * position and velocity at tf compared to the reference solution, the error in + * the position constraint equation, and various integrator statistics are + * printed to the screen for each run. + * + * When projection is enabled a user-supplied function is used to project the + * position, velocity, and error to the constraint manifold. + * + * Optional command line inputs may be used to change the final simulation time + * (default 30), the initial tolerance (default 1e-5), the number of outputs + * (default 1), or disable error projection. Use the option --help for a list + * of the command line flags. + * ---------------------------------------------------------------------------*/ + +#include +#include +#include +#include + +#include /* access to CVODE */ +#include /* access to serial N_Vector */ +#include /* access to dense SUNmatrix */ +#include /* access to dense SUNLinearSolver */ + +/* Precision specific formatting macros */ +#if defined(SUNDIALS_EXTENDED_PRECISION) +#define GSYM "Lg" +#define ESYM "Le" +#define FSYM "Lf" +#else +#define GSYM "g" +#define ESYM "e" +#define FSYM "f" +#endif + +/* Precision specific math function macros */ +#if defined(SUNDIALS_DOUBLE_PRECISION) +#define SIN(x) (sin((x))) +#define COS(x) (cos((x))) +#define SQRT(x) (sqrt((x))) +#define ABS(x) (fabs((x))) +#elif defined(SUNDIALS_SINGLE_PRECISION) +#define SIN(x) (sinf((x))) +#define COS(x) (cosf((x))) +#define SQRT(x) (sqrtf((x))) +#define ABS(x) (fabsf((x))) +#elif defined(SUNDIALS_EXTENDED_PRECISION) +#define SIN(x) (sinl((x))) +#define COS(x) (cosl((x))) +#define SQRT(x) (sqrtl((x))) +#define ABS(x) (fabsl((x))) +#endif + +/* Problem Constants */ +#define ZERO RCONST(0.0) +#define ONE RCONST(1.0) +#define GRAV RCONST(13.750371636040745654980191559621114395801712) + +/* Functions provided to CVODE */ +static int fref(realtype t, N_Vector yy, N_Vector fy, void *f_data); + +static int f(realtype t, N_Vector yy, N_Vector fy, void *f_data); +static int proj(realtype t, N_Vector yy, N_Vector corr, + realtype epsProj, N_Vector err, void *pdata); + +/* Functions to integrate the Cartesian and reference solutions */ +int GetSol(void *cvode_mem, N_Vector yy0, realtype tol, realtype tf, + int nout, booleantype proj, booleantype projerr, N_Vector yref); + +int RefSol(realtype tf, N_Vector yref, int nout); + +/* Utility functions */ +static int ReadInputs(int *argc, char ***argv, realtype *tol, realtype *tf, + int *nout, booleantype *projerr); +static void InputHelp(); +static int check_retval(void *returnvalue, const char *funcname, int opt); + +/* ----------------------------------------------------------------------------- + * Main Program + * ---------------------------------------------------------------------------*/ + +int main(int argc, char* argv[]) +{ + int i; + int retval; /* reusable return flag */ + int nout = 1; /* number of outputs */ + realtype tol = RCONST(1.0e-5); /* integration tolerance */ + realtype tf = RCONST(30.0); /* final integration time */ + booleantype projerr = SUNTRUE; /* enable error projection */ + + void *cvode_mem = NULL; /* CVODE memory */ + N_Vector yy0 = NULL; /* initial condition vector */ + realtype *yy0data = NULL; /* vector data */ + N_Vector yref = NULL; /* reference solution vector */ + SUNMatrix A = NULL; /* Jacobian matrix */ + SUNLinearSolver LS = NULL; /* linear solver */ + + /* Read command line inputs */ + retval = ReadInputs(&argc, &argv, &tol, &tf, &nout, &projerr); + if (check_retval(&retval, "ReadInputs", 1)) return(1); + + /* Compute reference solution */ + yref = N_VNew_Serial(4); + + retval = RefSol(tf, yref, nout); + if (check_retval(&retval, "RefSol", 1)) return(1); + + /* Create serial vector to store the initial condition */ + yy0 = N_VNew_Serial(4); + if (check_retval((void *)yy0, "N_VNew_Serial", 0)) return(1); + + /* Set the initial condition values */ + yy0data = N_VGetArrayPointer(yy0); + + yy0data[0] = ONE; /* x */ + yy0data[1] = ZERO; /* y */ + yy0data[2] = ZERO; /* xd */ + yy0data[3] = ZERO; /* yd */ + + /* Create CVODE memory */ + cvode_mem = CVodeCreate(CV_BDF); + if (check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1); + + /* Initialize CVODE */ + retval = CVodeInit(cvode_mem, f, ZERO, yy0); + if (check_retval(&retval, "CVodeInit", 1)) return(1); + + /* Set integration tolerances */ + retval = CVodeSStolerances(cvode_mem, tol, tol); + if (check_retval(&retval, "CVodeSStolerances", 1)) return(1); + + /* Create dense SUNMatrix for use in linear solves */ + A = SUNDenseMatrix(4, 4); + if(check_retval((void *)A, "SUNDenseMatrix", 0)) return(1); + + /* Create dense SUNLinearSolver object */ + LS = SUNLinSol_Dense(yy0, A); + if(check_retval((void *)LS, "SUNLinSol_Dense", 0)) return(1); + + /* Attach the matrix and linear solver to CVODE */ + retval = CVodeSetLinearSolver(cvode_mem, LS, A); + if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1); + + /* Set a user-supplied projection function */ + retval = CVodeSetProjFn(cvode_mem, proj); + if(check_retval(&retval, "CVodeSetProjFn", 1)) return(1); + + /* Set maximum number of steps between outputs */ + retval = CVodeSetMaxNumSteps(cvode_mem, 50000); + if (check_retval(&retval, "CVodeSetMaxNumSteps", 1)) return(1); + + /* Compute the solution with various tolerances */ + for (i = 0; i < 5; i++) { + + /* Output tolerance and output header for this run */ + printf("\n\nTol = %8.2" ESYM"\n", tol); + printf("Project x y"); + printf(" x' y' | g | "); + printf("nst rhs eval setups (J eval) | cf ef\n"); + + /* Compute solution with projection */ + retval = GetSol(cvode_mem, yy0, tol, tf, nout, SUNTRUE, projerr, yref); + if (check_retval(&retval, "GetSol", 1)) return(1); + + /* Compute solution without projection */ + retval = GetSol(cvode_mem, yy0, tol, tf, nout, SUNFALSE, SUNFALSE, yref); + if (check_retval(&retval, "GetSol", 1)) return(1); + + /* Reduce tolerance for next run */ + tol /= RCONST(10.0); + } + + /* Free memory */ + N_VDestroy_Serial(yref); + N_VDestroy_Serial(yy0); + SUNMatDestroy(A); + SUNLinSolFree(LS); + CVodeFree(&cvode_mem); + + return(0); +} + + +/* ----------------------------------------------------------------------------- + * Functions to integrate the Cartesian and reference systems + * ---------------------------------------------------------------------------*/ + + +/* Compute the Cartesian system solution */ +int GetSol(void *cvode_mem, N_Vector yy0, realtype tol, realtype tf, int nout, + booleantype proj, booleantype projerr, N_Vector yref) +{ + char outname[100]; /* output file name */ + FILE *FID = NULL; /* output file */ + N_Vector yy = NULL; /* solution vector */ + realtype *yydata = NULL; /* vector data */ + + int retval; /* reusable return flag */ + int out; /* output counter */ + realtype dtout; /* output frequency */ + realtype tout; /* output time */ + realtype t; /* return time */ + realtype x, y; /* position values */ + realtype xd, yd; /* velocity values */ + realtype g; /* constraint value */ + + /* Integrator stats */ + long int nst, nfe, nsetups, nje, nfeLS, ncfn, netf; + + /* Enable or disable projection */ + if (proj) + { + printf(" YES "); + retval = CVodeSetProjFrequency(cvode_mem, 1); + if(check_retval(&retval, "CVodeSetProjFrequency", 1)) return(1); + + /* Enable or disable error projection */ + retval = CVodeSetProjErrEst(cvode_mem, projerr); + if(check_retval(&retval, "CVodeSetProjErrEst", 1)) return(1); + } + else + { + retval = CVodeSetProjFrequency(cvode_mem, 0); + if(check_retval(&retval, "CVodeSetProjFrequency", 1)) return(1); + printf(" NO "); + } + + /* Create vector to store the solution */ + yy = N_VNew_Serial(4); + + /* Copy initial condition into solution vector */ + N_VScale(ONE, yy0, yy); + + /* Get pointer to vector data */ + yydata = N_VGetArrayPointer(yy); + + /* Reinitialize CVODE for this run */ + retval = CVodeReInit(cvode_mem, ZERO, yy0); + if (check_retval(&retval, "CVodeReInit", 1)) + { + N_VDestroy_Serial(yy); + return(retval); + } + + /* Set integration tolerances for this run */ + retval = CVodeSStolerances(cvode_mem, tol, tol); + if (check_retval(&retval, "CVodeSStolerances", 1)) + { + N_VDestroy_Serial(yy); + return(retval); + } + + /* Open output file */ + if (proj) + { + sprintf(outname, "cvPendulum_dns_tol_%03.2" ESYM"_proj.txt", tol); + } + else + { + sprintf(outname, "cvPendulum_dns_tol_%03.2" ESYM".txt", tol); + } + FID = fopen(outname, "w"); + + /* Output initial condition */ + fprintf(FID, + "%0.4" ESYM" %14.6" ESYM" %14.6" ESYM" %14.6" ESYM" %14.6" ESYM"\n", + ZERO, yydata[0], yydata[1], yydata[2], yydata[3]); + + /* Integrate to tf and peridoically output the solution */ + dtout = tf / nout; + tout = dtout; + + for (out = 0; out < nout; out++) + { + /* Set stop time (do not interpolate output) */ + retval = CVodeSetStopTime(cvode_mem, tout); + if (check_retval(&retval, "CVodeSetStopTime", 1)) + { + N_VDestroy_Serial(yy); + fclose(FID); + return(retval); + } + + /* Integrate to tout */ + retval = CVode(cvode_mem, tout, yy, &t, CV_NORMAL); + if (check_retval(&retval, "CVode", 1)) + { + N_VDestroy_Serial(yy); + fclose(FID); + return(retval); + } + + /* Write output */ + fprintf(FID, + "%0.4" ESYM" %14.6" ESYM" %14.6" ESYM" %14.6" ESYM" %14.6" ESYM"\n", + t, yydata[0], yydata[1], yydata[2], yydata[3]); + + /* Update output time */ + if (out < nout - 1) + { + tout += dtout; + } + else + { + tout = tf; + } + } + + /* Close output file */ + fclose(FID); + + /* Compute the constraint violation */ + x = yydata[0]; + y = yydata[1]; + g = ABS(x*x + y*y - ONE); + + /* Compute the absolute error compared to the reference solution */ + N_VLinearSum(ONE, yy, -ONE, yref, yy); + N_VAbs(yy, yy); + + x = yydata[0]; + y = yydata[1]; + xd = yydata[2]; + yd = yydata[3]; + + /* Output errors */ + printf("%8.2" ESYM" %8.2" ESYM" %8.2" ESYM" %8.2" ESYM" | %8.2" ESYM" |", + x, y, xd, yd, g); + + /* Free solution vector */ + N_VDestroy_Serial(yy); + + /* Get integrator stats */ + retval = CVodeGetNumSteps(cvode_mem, &nst); + if (check_retval(&retval, "CVodeGetNumSteps", 1)) return(retval); + + retval = CVodeGetNumRhsEvals(cvode_mem, &nfe); + if (check_retval(&retval, "CVodeGetNumFctEvals", 1)) return(retval); + + retval = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups); + if (check_retval(&retval, "CVodeGetNumLinSolvSetups", 1)) return(retval); + + retval = CVodeGetNumErrTestFails(cvode_mem, &netf); + if (check_retval(&retval, "CVodeGetNumErrTestFails", 1)) return(retval); + + retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn); + if (check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1)) return(retval); + + retval = CVodeGetNumJacEvals(cvode_mem, &nje); + if (check_retval(&retval, "CVodeGetNumJacEvals", 1)) return(retval); + + retval = CVodeGetNumLinRhsEvals(cvode_mem, &nfeLS); + if (check_retval(&retval, "CVodeGetNumLinRhsEvals", 1)) return(retval); + + /* Output stats */ + printf(" %6ld %6ld+%-4ld %4ld (%3ld) | %3ld %3ld\n", + nst, nfe, nfeLS, nsetups, nje, ncfn, netf); + + return(0); +} + + +/* Compute the reference system solution */ +int RefSol(realtype tf, N_Vector yref, int nout) +{ + FILE *FID = NULL; /* output file */ + void *cvode_mem = NULL; /* CVODE memory */ + N_Vector yy = NULL; /* solution vector */ + realtype *yydata = NULL; /* vector data */ + SUNMatrix A = NULL; /* Jacobian matrix */ + SUNLinearSolver LS = NULL; /* linear solver */ + + int retval; /* reusable return flag */ + int out; /* output counter */ + realtype dtout; /* output frequency */ + realtype tout; /* output time */ + realtype t; /* return time */ + realtype th, thd; /* theta and theta dot */ + realtype tol = RCONST(1.0e-14); /* integration tolerance */ + + /* Create the solution vector */ + yy = N_VNew_Serial(2); + if (check_retval((void *)yy, "N_VNew_Serial", 0)) return(-1); + + /* Set the initial condition */ + yydata = N_VGetArrayPointer(yy); + + yydata[0] = ZERO; /* theta */ + yydata[1] = ZERO; /* theta' */ + + /* Create CVODE memory */ + cvode_mem = CVodeCreate(CV_BDF); + if (check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1); + + /* Initialize CVODE */ + retval = CVodeInit(cvode_mem, fref, ZERO, yy); + if (check_retval(&retval, "CVodeInit", 1)) return(1); + + /* Set integration tolerances */ + retval = CVodeSStolerances(cvode_mem, tol, tol); + if (check_retval(&retval, "CVodeSStolerances", 1)) return(1); + + /* Create dense SUNMatrix for use in linear solves */ + A = SUNDenseMatrix(2, 2); + if(check_retval((void *)A, "SUNDenseMatrix", 0)) return(1); + + /* Create dense SUNLinearSolver object */ + LS = SUNLinSol_Dense(yy, A); + if(check_retval((void *)LS, "SUNLinSol_Dense", 0)) return(1); + + /* Attach the matrix and linear solver to CVODE */ + retval = CVodeSetLinearSolver(cvode_mem, LS, A); + if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1); + + /* Set CVODE optional inputs */ + retval = CVodeSetMaxNumSteps(cvode_mem, 100000); + if (check_retval(&retval, "CVodeSetMaxNumSteps", 1)) return(1); + + retval = CVodeSetStopTime(cvode_mem, tf); + if (check_retval(&retval, "CVodeSetStopTime", 1)) return(1); + + /* Open output file */ + FID = fopen("cvPendulum_dns_ref.txt", "w"); + + /* Output initial condition */ + th = yydata[0]; + thd = yydata[1]; + fprintf(FID, + "%0.4" ESYM" %14.6" ESYM" %14.6" ESYM" %14.6" ESYM" %14.6" ESYM"\n", + ZERO, COS(th), SIN(th), -thd * SIN(th), thd * COS(th)); + + /* Integrate to tf and periodically output the solution */ + dtout = tf / nout; + tout = dtout; + + for (out = 0; out < nout; out++) + { + /* Set stop time (do not interpolate output) */ + retval = CVodeSetStopTime(cvode_mem, tout); + if (check_retval(&retval, "CVodeSetStopTime", 1)) + { + N_VDestroy_Serial(yy); + SUNMatDestroy(A); + SUNLinSolFree(LS); + CVodeFree(&cvode_mem); + fclose(FID); + return(retval); + } + + /* Integrate to tout */ + retval = CVode(cvode_mem, tf, yy, &t, CV_NORMAL); + if (check_retval(&retval, "CVode", 1)) + { + N_VDestroy_Serial(yy); + SUNMatDestroy(A); + SUNLinSolFree(LS); + CVodeFree(&cvode_mem); + fclose(FID); + return(retval); + } + + /* Write output */ + th = yydata[0]; + thd = yydata[1]; + fprintf(FID, + "%0.4" ESYM" %14.6" ESYM" %14.6" ESYM" %14.6" ESYM" %14.6" ESYM"\n", + t, COS(th), SIN(th), -thd * SIN(th), thd * COS(th)); + + /* Update output time */ + if (out < nout - 1) + { + tout += dtout; + } + else + { + tout = tf; + } + } + + /* Close output file */ + fclose(FID); + + /* Get solution components */ + th = yydata[0]; + thd = yydata[1]; + + /* Convert to Cartesian reference solution */ + yydata = N_VGetArrayPointer(yref); + + yydata[0] = COS(th); + yydata[1] = SIN(th); + yydata[2] = -thd * SIN(th); + yydata[3] = thd * COS(th); + + /* Free memory */ + N_VDestroy_Serial(yy); + SUNMatDestroy(A); + SUNLinSolFree(LS); + CVodeFree(&cvode_mem); + + return(0); +} + + +/* ----------------------------------------------------------------------------- + * Functions provided to CVODE + * ---------------------------------------------------------------------------*/ + + +/* ODE RHS function for the reference system */ +static int fref(realtype t, N_Vector yy, N_Vector fy, void *f_data) +{ + realtype *yydata = NULL; /* yy vector data */ + realtype *fydata = NULL; /* fy vector data */ + + /* Get vector array pointers */ + yydata = N_VGetArrayPointer(yy); + fydata = N_VGetArrayPointer(fy); + + fydata[0] = yydata[1]; /* theta' */ + fydata[1] = -GRAV * COS(yydata[0]); /* -g * cos(theta) */ + return(0); +} + + +/* ODE RHS function for the Cartesian system */ +static int f(realtype t, N_Vector yy, N_Vector fy, void *f_data) +{ + realtype *yydata = NULL; /* yy vector data */ + realtype *fydata = NULL; /* fy vector data */ + + realtype x, y; /* positions */ + realtype xd, yd; /* velocities */ + realtype tmp; + + /* Get vector array pointers */ + yydata = N_VGetArrayPointer(yy); + fydata = N_VGetArrayPointer(fy); + + /* Get vector components */ + x = yydata[0]; + y = yydata[1]; + xd = yydata[2]; + yd = yydata[3]; + + /* Compute tension */ + tmp = xd * xd + yd * yd - GRAV * y; + + /* Compute RHS */ + fydata[0] = xd; + fydata[1] = yd; + fydata[2] = -x * tmp; + fydata[3] = -y * tmp - GRAV; + + return(0); +} + + +/* Projection function */ +static int proj(realtype t, N_Vector yy, N_Vector corr, + realtype epsProj, N_Vector err, void *pdata) +{ + realtype *yydata = NULL; /* yy vector data */ + realtype *cdata = NULL; /* corr vector data */ + realtype *edata = NULL; /* err vector data */ + + realtype x, y, x_new, y_new; /* positions */ + realtype xd, yd, xd_new, yd_new; /* velocities */ + + realtype e1, e2, e3, e4; + realtype e1_new, e2_new, e3_new, e4_new; + realtype R; + + /* Get vector array pointers */ + + yydata = N_VGetArrayPointer(yy); + cdata = N_VGetArrayPointer(corr); + + /* Extract current solution */ + + x = yydata[0]; + y = yydata[1]; + xd = yydata[2]; + yd = yydata[3]; + + /* Project positions */ + + R = SQRT(x * x + y * y); + + x_new = x / R; + y_new = y / R; + + /* Project velocities + * + * +- -+ +- -+ + * | y*y -x*y | | xd | + * P v = | | | | + * | -x*y x*x | | yd | + * +- -+ +- -+ + */ + + xd_new = xd * y_new * y_new - yd * x_new * y_new; + yd_new = - xd * x_new * y_new + yd * x_new * x_new; + + /* Return position and velocity corrections */ + + cdata[0] = x_new - x; + cdata[1] = y_new - y; + cdata[2] = xd_new - xd; + cdata[3] = yd_new - yd; + + /* Project error P * err */ + if (err != NULL) + { + edata = N_VGetArrayPointer(err); + + e1 = edata[0]; + e2 = edata[1]; + e3 = edata[2]; + e4 = edata[3]; + + e1_new = y_new * y_new * e1 - x_new * y_new * e2; + e2_new = -x_new * y_new * e1 + x_new * x_new * e2; + + e3_new = y_new * y_new * e3 - x_new * y_new * e4; + e4_new = -x_new * y_new * e3 + x_new * x_new * e4; + + edata[0] = e1_new; + edata[1] = e2_new; + edata[2] = e3_new; + edata[3] = e4_new; + } + + return(0); +} + + +/* ----------------------------------------------------------------------------- + * Private helper functions + * ---------------------------------------------------------------------------*/ + + +/* Read command line unputs */ +static int ReadInputs(int *argc, char ***argv, realtype *tol, realtype *tf, + int *nout, booleantype *projerr) +{ + int arg_idx = 1; + + /* check for input args */ + while (arg_idx < (*argc)) + { + if (strcmp((*argv)[arg_idx],"--tol") == 0) + { + arg_idx++; + *tol = atof((*argv)[arg_idx++]); + } + else if (strcmp((*argv)[arg_idx],"--tf") == 0) + { + arg_idx++; + *tf = atof((*argv)[arg_idx++]); + } + else if (strcmp((*argv)[arg_idx],"--nout") == 0) + { + arg_idx++; + *nout = atoi((*argv)[arg_idx++]); + } + else if (strcmp((*argv)[arg_idx],"--noerrproj") == 0) + { + arg_idx++; + *projerr = SUNFALSE; + } + else if (strcmp((*argv)[arg_idx],"--help") == 0 ) + { + InputHelp(); + return(-1); + } + else + { + fprintf(stderr, "ERROR: Invalid input %s",(*argv)[arg_idx]); + InputHelp(); + return(-1); + } + } + + return(0); +} + + +/* Print command line options */ +static void InputHelp() +{ + printf("\nCommand line options:\n"); + printf(" --tol : relative and absolute tolerance\n"); + printf(" --tf