SUNDIALS cleanup

Exec/GNUmakefile

@@ -11,43 +11,14 @@ DIM          = 3
 USE_FFT      = TRUE
 
 USE_SUNDIALS = FALSE
+SUNDIALS_HOME ?= ../../sundials/instdir
 
 include $(AMREX_HOME)/Tools/GNUMake/Make.defs
 
 include ../Source/Make.package
 VPATH_LOCATIONS  += ../Source
 INCLUDE_LOCATIONS += ../Source
 
-ifeq ($(USE_SUNDIALS),TRUE)
-  ifeq ($(USE_CUDA),TRUE)
-    SUNDIALS_ROOT ?= $(TOP)../../sundials/instdir_cuda
-  else
-    SUNDIALS_ROOT ?= $(TOP)../../sundials/instdir
-  endif
-  ifeq ($(NERSC_HOST),perlmutter)
-    SUNDIALS_LIB_DIR ?= $(SUNDIALS_ROOT)/lib64
-  else
-    SUNDIALS_LIB_DIR ?= $(SUNDIALS_ROOT)/lib
-  endif
-
-  USE_CVODE_LIBS ?= TRUE
-  USE_ARKODE_LIBS ?= TRUE
-
-  DEFINES += -DAMREX_USE_SUNDIALS
-  INCLUDE_LOCATIONS += $(SUNDIALS_ROOT)/include
-  LIBRARY_LOCATIONS += $(SUNDIALS_LIB_DIR)
-
-  LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_cvode
-  LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_arkode
-  LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_nvecmanyvector
-  LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_core
-
-  ifeq ($(USE_CUDA),TRUE)
-    LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_nveccuda
-  endif
-
-endif
-
 include $(AMREX_HOME)/Src/Base/Make.package
 include $(AMREX_HOME)/Src/FFT/Make.package
 include $(AMREX_HOME)/Src/Boundary/Make.package 

Exec/README_sundials

@@ -1,94 +1,4 @@
-SUNDIALS installation guide:
-https://computing.llnl.gov/projects/sundials/faq#inst
+Refer to https://amrex-codes.github.io/amrex/docs_html/TimeIntegration_Chapter.html
+for SNUDIALS installation, build, and usage instructions.
 
-Installation
-
-# You need SUNDIALS v7.1.1 or later.
-# Check https://computing.llnl.gov/projects/sundials/sundials-software to see if it's available for download.
-# If so, download sundials-x.y.z.tar.gz and extract it at the same level as amrex using
->> tar -xzvf sundials-x.y.z.tar.gz # where x.y.z is the version of sundials
->> mv sundials-x.y.z sundials-src
-
-# If v7.1.1. is not available on the website, clone the git repo directly and use the latest version
-# At the same level that amrex is cloned, do:
-
->> git clone https://github.com/LLNL/sundials.git
->> mv sundials sundials-src
-
-# Next
-
->> mkdir sundials
->> cd sundials
-
-######################
-HOST BUILD 
-######################
-
->> mkdir instdir
->> mkdir builddir
->> cd builddir
-
->> cmake -DCMAKE_INSTALL_PREFIX=/pathto/sundials/instdir -DEXAMPLES_INSTALL_PATH=/pathto/sundials/instdir/examples -DENABLE_MPI=ON ../../sundials-src
-
-######################
-NVIDIA/CUDA BUILD
-######################
-
-# Navigate back to the 'sundials' directory and do:
-
->> mkdir instdir_cuda
->> mkdir builddir_cuda
->> cd builddir_cuda
-
->> cmake -DCMAKE_INSTALL_PREFIX=/pathto/sundials/instdir_cuda -DEXAMPLES_INSTALL_PATH=/pathto/sundials/instdir_cuda/examples -DENABLE_CUDA=ON -DENABLE_MPI=ON ../../sundials-src
-
-######################
-
->> make -j4
->> make install
-
-# in your .bashrc or preferred configuration file, add the following (and then "source ~/.bashrc")
-
-# If you have a CPU build:
-
->> export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/pathto/sundials/instdir/lib/
-
-# If you have a NVIDIA/CUDA build:
-
->> export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/pathto/sundials/instdir_cuda/lib/
-
-# now you are ready to compile MagneX with:
-
->> make -j4 USE_SUNDIALS=TRUE # optional to have 'USE_CUDA=TRUE' as well
-
-# in your inputs file, you will need to have:
-
-TimeIntegratorOption = 4 #amrex/sundials backend integrators
-
-# INTEGRATION
-## integration.type can take on the following values:
-## 0 or "ForwardEuler" => Native AMReX Forward Euler integrator
-## 1 or "RungeKutta"   => Native AMReX Explicit Runge Kutta controlled by integration.rk.type
-## 2 or "SUNDIALS"     => SUNDIALS backend controlled by integration.sundials.strategy
-integration.type = SUNDIALS
-
-## Native AMReX Explicit Runge-Kutta parameters
-#
-## integration.rk.type can take the following values:
-### 0 = User-specified Butcher Tableau
-### 1 = Forward Euler
-### 2 = Trapezoid Method
-### 3 = SSPRK3 Method
-### 4 = RK4 Method
-integration.rk.type = 1
-
-# Set the SUNDIALS method type:
-# ERK      = Explicit Runge-Kutta method
-# DIRK     = Diagonally Implicit Runge-Kutta method
-#
-# Optionally select a specific SUNDIALS method by name, see the SUNDIALS
-# documentation for the supported method names
-
-# Use forward Euler (fixed step sizes only)
-integration.sundials.type = ERK
-integration.sundials.method = ARKODE_FORWARD_EULER_1_1
+Make sure that SUNDIALS_HOME in the GNUmakefile points to the installation directory.

Exec/standard_problem_inputs/inputs_std4_IMEX

@@ -0,0 +1,120 @@
+amrex.use_gpu_aware_mpi=1
+
+##############
+# 0.78125nm case
+##############
+n_cell = 640 160 4
+max_grid_size_x = 640
+max_grid_size_y = 160
+max_grid_size_z = 4
+dt = 5.e-15 # IMEX, 1.e-14 unstable
+
+stop_time = 4.e-9
+
+plot_int = 25000
+chk_int = -1
+restart = -1
+diag_type = 4
+
+prob_lo = 0. 0. 0.
+prob_hi = 500.e-9 125.e-9 3.125e-9
+
+mu0 = 1.25663e-6
+
+Mx_parser(x,y,z) = "8.e5 * (x>0.)*(x<500.e-9)*(y>0.)*(y<125.e-9)*(z>0.)*(z<3.125e-9)"
+My_parser(x,y,z) = "0."
+Mz_parser(x,y,z) = "0."
+
+# Field 1:  mu_0 Hx=-24.6 mT, mu_0 Hy= 4.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 25 mT, directed 170 degrees counterclockwise from the positive x axis
+timedependent_Hbias = 1
+Hx_bias_parser(x,y,z,t) = "(t<=2.e-11)*1.e5 + (t>2.e-11)*(t<=3.e-11)*(3.e-11-t)*1.e16 - (t>1.e-9)*19576"
+Hy_bias_parser(x,y,z,t) = "(t<=2.e-11)*1.e5 + (t>2.e-11)*(t<=3.e-11)*(3.e-11-t)*1.e16 + (t>1.e-9)*3422"
+Hz_bias_parser(x,y,z,t) = "(t<=2.e-11)*1.e5 + (t>2.e-11)*(t<=3.e-11)*(3.e-11-t)*1.e16"
+
+# Field 2:  mu_0 Hx=-35.5 mT, mu_0 Hy=-6.3 mT, mu_0 Hz= 0.0 mT
+# which is a field approximately 36 mT, directed 190 degrees counterclockwise from the positive x axis
+#timedependent_Hbias = 1
+#Hx_bias_parser(x,y,z,t) = "(t<=2.e-11)*1.e5 + (t>2.e-11)*(t<=3.e-11)*(3.e-11-t)*1.e16 - (t>1.e-9)*28259"
+#Hy_bias_parser(x,y,z,t) = "(t<=2.e-11)*1.e5 + (t>2.e-11)*(t<=3.e-11)*(3.e-11-t)*1.e16 - (t>1.e-9)*5013"
+#Hz_bias_parser(x,y,z,t) = "(t<=2.e-11)*1.e5 + (t>2.e-11)*(t<=3.e-11)*(3.e-11-t)*1.e16"
+
+timedependent_alpha = 1
+alpha_parser(x,y,z,t) =  " (x>0.)*(x<500.e-9)*(y>0.)*(y<125.e-9)*(z>0.)*(z<3.125e-9) * (0.5*(t<=1.e-9) + 0.02*(t>1.e-9))"
+Ms_parser(x,y,z) =       " (x>0.)*(x<500.e-9)*(y>0.)*(y<125.e-9)*(z>0.)*(z<3.125e-9) * 8.e5"
+gamma_parser(x,y,z) =    " (x>0.)*(x<500.e-9)*(y>0.)*(y<125.e-9)*(z>0.)*(z<3.125e-9) * -1.759e11"
+exchange_parser(x,y,z) = " (x>0.)*(x<500.e-9)*(y>0.)*(y<125.e-9)*(z>0.)*(z<3.125e-9)* 1.3e-11"
+anisotropy_parser(x,y,z) = "0."
+DMI_parser(x,y,z) = "0."
+
+precession = 1
+demag_coupling = 1
+FFT_solver = 1
+M_normalization = 1  # 0 = unsaturated case; 1 = saturated case
+exchange_coupling = 1
+anisotropy_coupling = 0
+anisotropy_axis = 0.0 1.0 0.0
+DMI_coupling = 0
+
+# INTEGRATION
+#TimeIntegratorOption = 1 #Forward Euler
+#TimeIntegratorOption = 2 #Predictor-corrector
+#TimeIntegratorOption = 3 #2nd order artemis way
+TimeIntegratorOption = 4 #amrex/sundials backend integrators
+
+# tolerance threshold (L_inf change between iterations) for TimeIntegrationOption 2 and 3
+iterative_tolerance = 1.e-9
+
+# for TimeIntegratorOption = 4, integration.type can take on the following values:
+## 0 or "ForwardEuler" => Native AMReX Forward Euler integrator
+## 1 or "RungeKutta"   => Native AMReX Explicit Runge Kutta controlled by integration.rk.type
+## 2 or "SUNDIALS"     => SUNDIALS backend controlled by integration.sundials.type
+integration.type = SUNDIALS
+
+## for integration.type = RungeKutta, integration.rk.type can take the following values:
+### 0 = User-specified Butcher Tableau
+### 1 = Forward Euler
+### 2 = Trapezoid Method
+### 3 = SSPRK3 Method
+### 4 = RK4 Method
+integration.rk.type = 4
+
+## If using a user-specified Butcher Tableau, then
+## set nodes, weights, and table entries here:
+#
+## The Butcher Tableau is read as a flattened,
+## lower triangular matrix (but including the diagonal)
+## in row major format.
+integration.rk.weights = 1
+integration.rk.nodes = 0
+integration.rk.tableau = 0.0
+
+## for integration.type = SUNDIALS, set the SUNDIALS method type:
+# ERK      = Explicit Runge-Kutta method
+# DIRK     = Diagonally Implicit Runge-Kutta method
+# IMEX-RK  = Implicit-Explicit Additive Runge-Kutta method
+# EX-MRI   = Explicit Multirate Infatesimal method
+# IM-MRI   = Implicit Multirate Infatesimal method
+# IMEX-MRI = Implicit-Explicit Multirate Infatesimal method
+integration.sundials.type = IMEX-RK
+
+# *** Select a specific SUNDIALS ERK or DIRK method ***
+# https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html#explicit-butcher-tables
+# https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html#implicit-butcher-tables
+#integration.sundials.method = ARKODE_FORWARD_EULER_1_1
+#integration.sundials.method = ARKODE_BACKWARD_EULER_1_1
+
+# *** Select a specific SUNDIALS IMEX-RK method ***
+# https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html#additive-butcher-tables
+integration.sundials.method_i = ARKODE_ARK2_DIRK_3_1_2
+integration.sundials.method_e = ARKODE_ARK2_ERK_3_1_2
+
+# *** Select a specific SUNDIALS EX-MRI, IM-MRI, or IMEX-MRI method ***
+# https://sundials.readthedocs.io/en/latest/arkode/Usage/MRIStep/MRIStepCoupling.html#mri-coupling-tables
+#  *** Select a specific SUNDIALS MRI slow method coupling table ***
+#integration.sundials.method = ARKODE_MIS_KW3
+#  *** Select a specific SUNDIALS ERK or DIRK fast method ***
+#integration.sundials.fast_type = ERK  # ERK or DIRK
+#integration.sundials.fast_method = ARKODE_KNOTH_WOLKE_3_3
+
+implicit_exchange = 1

Exec/standard_problem_inputs/inputs_std4_MRI

@@ -7,7 +7,7 @@ n_cell = 640 160 4
 max_grid_size_x = 640
 max_grid_size_y = 160
 max_grid_size_z = 4
-dt = 1.25e-13 # SUNDIALS MRI
+dt = 1.25e-13 # SUNDIALS MRI, 2.5e-13 unstable
 
 stop_time = 4.e-9
 

Exec/standard_problem_inputs/inputs_std4_RK4

@@ -7,11 +7,11 @@ n_cell = 640 160 4
 max_grid_size_x = 640
 max_grid_size_y = 160
 max_grid_size_z = 4
-dt = 2.5e-14 # RK4 (438s) 5.e-14 unstable
+dt = 2.5e-14 # RK4, 5.e-14 unstable
 
 stop_time = 4.e-9
 
-plot_int = 10000
+plot_int = 5000
 chk_int = -1
 restart = -1
 diag_type = 4

Source/main.cpp

@@ -558,39 +558,47 @@ void main_main ()
                 }
 
                 // exchange
-                if ( (using_MRI==1 && fast_exchange) || (using_IMEX==1 && implicit_exchange) ) {
-                    for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
-                        H_exchangefield[idim].setVal(0.);
+                if (exchange_coupling == 1) {
+                    if ( (using_MRI==1 && fast_exchange) || (using_IMEX==1 && implicit_exchange) ) {
+                        for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
+                            H_exchangefield[idim].setVal(0.);
+                        }
+                    } else {
+                        CalculateH_exchange(ar_state, H_exchangefield, Ms, exchange, DMI, geom);
                     }
-                } else {
-                    CalculateH_exchange(ar_state, H_exchangefield, Ms, exchange, DMI, geom);
                 }
 
                 // DMI
-                if ( (using_MRI==1 && fast_DMI) || (using_IMEX==1 && implicit_DMI) ) {
-                    for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
-                        H_DMIfield[idim].setVal(0.);
+                if (DMI_coupling == 1) {
+                    if ( (using_MRI==1 && fast_DMI) || (using_IMEX==1 && implicit_DMI) ) {
+                        for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
+                            H_DMIfield[idim].setVal(0.);
+                        }
+                    } else {
+                        CalculateH_DMI(ar_state, H_DMIfield, Ms, exchange, DMI, geom);
                     }
-                } else {
-                    CalculateH_DMI(ar_state, H_DMIfield, Ms, exchange, DMI, geom);
                 }
 
                 // anisotropy
-                if ( (using_MRI==1 && fast_anisotropy) || (using_IMEX==1 && implicit_anisotropy) ) {
-                    for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
-                        H_anisotropyfield[idim].setVal(.0);
+                if (anisotropy_coupling == 1) {
+                    if ( (using_MRI==1 && fast_anisotropy) || (using_IMEX==1 && implicit_anisotropy) ) {
+                        for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
+                            H_anisotropyfield[idim].setVal(.0);
+                        }
+                    } else {
+                        CalculateH_anisotropy(ar_state, H_anisotropyfield, Ms, anisotropy);
                     }
-                } else {
-                    CalculateH_anisotropy(ar_state, H_anisotropyfield, Ms, anisotropy);
                 }
 
                 // H_demag
-                if ( (using_MRI==1 && fast_demag) || (using_IMEX==1 && implicit_demag) ) {
-                    for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
-                        H_demagfield[idim].setVal(0.);
+                if (demag_coupling == 1) {
+                    if ( (using_MRI==1 && fast_demag) || (using_IMEX==1 && implicit_demag) ) {
+                        for (int idim=0; idim<AMREX_SPACEDIM; ++idim) {
+                            H_demagfield[idim].setVal(0.);
+                        }
+                    } else {
+                        demag_solver.CalculateH_demag(ar_state, H_demagfield);
                     }
-                } else {
-                    demag_solver.CalculateH_demag(ar_state, H_demagfield);
                 }
 
                 // Compute f^n = f(M^n, H^n)