Updating serial version of example 23

Add IMEX functionality to advection-diffusion example (ex23)

.appveyor.yml

@@ -43,7 +43,9 @@ before_build:
 build_script:
 - cmake --build build_parallel
 - cmake --build build_serial
+- cmake --build build_serial --target exec
 
 after_build:
 # - cmake --build build_parallel --target check
-- cmake --build build_serial --target check
+- cmake --build build_serial --target RUN_TESTS
+

.travis.yml

@@ -205,6 +205,7 @@ install:
         else
             echo "Reusing cached hypre-2.10.0b/";
         fi;
+        ln -s hypre-2.10.0b hypre;
      else
         echo "Serial build, not using hypre";
      fi

CHANGELOG

@@ -8,22 +8,21 @@
                                http://mfem.org
 
 
-Version 4.0-RC1, Apr 11, 2019
+Version 4.0-RC2, Apr 24, 2019
 =============================
 
 Requirements and Limitations
 ----------------------------
 - This is a release candidate for mfem-4.0.
-- Use at your own risk -- not everything will work, the API may change.
+- Use at your own risk -- not everything will work and the API may change.
 - We are looking for feedback from friendly users.
 - Unlike previous MFEM releases, this version requires a C++11 compiler.
 
 - GPU-related limitations:
-  * NVCC is not supported in the CMake build system yet.
-  * Element batching is currently ignored.
+  * Hypre preconditioners are not yet available in GPU mode.
+  * Only constant coefficients are currently supported on GPUs.
   * Full-assembly (on device), element assembly, and matrix-free bilinear forms
-    are not supported yet.
-  * FunctionCoefficients do not currently work on GPUs.
+    are not supported yet. Element batching is currently ignored.
   * Partial assembly kernels are not implemented yet for simplices.
 
 GPU support
@@ -145,6 +144,10 @@ New and improved solvers and preconditioners
 
 Miscellaneous
 -------------
+- In SparseMatrix added the option to perform MultTranspose() by matvec with
+  computed and stored transpose matrix. This is required for deterministic
+  results when using devices such as CUDA and OpenMP.
+
 - Added unit tests based on the Catch++ library.
 
 - Renamed the option MFEM_USE_OPENMP to MFEM_USE_LEGACY_OPENMP. This legacy

CMakeLists.txt

@@ -86,6 +86,13 @@ include("${CMAKE_CURRENT_SOURCE_DIR}/config/XSDKDefaults.cmake")
 
 # Enable languages.
 enable_language(CXX)
+if (MFEM_USE_CUDA)
+   # MFEM_USE_CUDA requires CMake 3.8 or newer (for direct CUDA support)
+   cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
+   enable_language(CUDA)
+   message(STATUS "Using CUDA architecture: ${CUDA_ARCH}")
+endif()
+
 if (XSDK_ENABLE_C)
    enable_language(C)
 endif()
@@ -266,6 +273,32 @@ if (MFEM_USE_PUMI)
   endif()
 endif()
 
+# CUDA
+if (MFEM_USE_CUDA)
+   set(CMAKE_CUDA_STANDARD 11)
+   set(CMAKE_CUDA_STANDARD_REQUIRED ON)
+   set(CMAKE_CUDA_EXTENSIONS OFF)
+   set(CMAKE_CUDA_FLAGS "-arch=${CUDA_ARCH} --expt-extended-lambda"
+       CACHE STRING "CUDA flags set for MFEM" FORCE)
+   if (MFEM_USE_MPI)
+      set(CUDA_CCBIN_COMPILER ${MPI_CXX_COMPILER})
+   else()
+      set(CUDA_CCBIN_COMPILER ${CMAKE_CXX_COMPILER})
+   endif()
+   string(APPEND CMAKE_CUDA_FLAGS " -ccbin ${CUDA_CCBIN_COMPILER}")
+   set(MFEM_USE_MM YES CACHE BOOL "Enable MFEM's memory manager" FORCE)
+endif()
+
+# OCCA
+if (MFEM_USE_OCCA)
+   find_package(OCCA REQUIRED)
+endif()
+
+# RAJA
+if (MFEM_USE_RAJA)
+   find_package(RAJA REQUIRED)
+endif()
+
 # MFEM_TIMER_TYPE
 if (NOT DEFINED MFEM_TIMER_TYPE)
   if (APPLE)
@@ -291,7 +324,7 @@ endif()
 #    be before SuiteSparse.
 set(MFEM_TPLS MPI_CXX OPENMP BLAS LAPACK METIS HYPRE SuiteSparse SUNDIALS PETSC
     MESQUITE SuperLUDist STRUMPACK AXOM CONDUIT GECKO GNUTLS NETCDF MPFR PUMI
-    POSIXCLOCKS MFEMBacktrace ZLIB)
+    POSIXCLOCKS MFEMBacktrace ZLIB OCCA RAJA)
 # Add all *_FOUND libraries in the variable TPL_LIBRARIES.
 set(TPL_LIBRARIES "")
 set(TPL_INCLUDE_DIRS "")
@@ -327,6 +360,13 @@ set(MFEM_SOURCE_DIRS general linalg mesh fem)
 foreach(DIR IN LISTS MFEM_SOURCE_DIRS)
   add_subdirectory(${DIR})
 endforeach()
+
+if (MFEM_USE_CUDA)
+   foreach(file IN LISTS SOURCES)
+      set_property(SOURCE ${file} PROPERTY LANGUAGE CUDA)
+   endforeach()
+endif()
+
 add_subdirectory(config)
 set(MASTER_HEADERS
   ${PROJECT_SOURCE_DIR}/mfem.hpp
@@ -337,6 +377,11 @@ set(CMAKE_INSTALL_RPATH_USE_LINK_PATH ON CACHE BOOL "")
 set(CMAKE_INSTALL_RPATH "${_lib_path}" CACHE PATH "")
 set(CMAKE_INSTALL_NAME_DIR "${_lib_path}" CACHE PATH "")
 
+set(MFEM_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR} CACHE PATH
+    "The MFEM source directory" FORCE)
+set(MFEM_INSTALL_DIR ${CMAKE_INSTALL_PREFIX} CACHE PATH
+    "The MFEM install directory" FORCE)
+
 # Declaring the library
 add_library(mfem ${SOURCES} ${HEADERS} ${MASTER_HEADERS})
 # message(STATUS "TPL_LIBRARIES = ${TPL_LIBRARIES}")
@@ -434,12 +479,12 @@ endif()
 # Add 'check' target - quick test
 if (NOT MFEM_USE_MPI)
   add_custom_target(check
-    ${CMAKE_CTEST_COMMAND} -R '^ex1_ser' -C ${CMAKE_CFG_INTDIR}
+    ${CMAKE_CTEST_COMMAND} -R \"^ex1_ser\" -C ${CMAKE_CFG_INTDIR}
     USES_TERMINAL)
   add_dependencies(check ex1)
 else()
   add_custom_target(check
-    ${CMAKE_CTEST_COMMAND} -R '^ex1p' -C ${CMAKE_CFG_INTDIR}
+    ${CMAKE_CTEST_COMMAND} -R \"^ex1p\" -C ${CMAKE_CFG_INTDIR}
     USES_TERMINAL)
   add_dependencies(check ex1p)
 endif()
@@ -484,6 +529,13 @@ install(DIRECTORY ${MFEM_SOURCE_DIRS}
   DESTINATION ${INSTALL_INCLUDE_DIR}/mfem
   FILES_MATCHING PATTERN "*.hpp")
 
+# Install the okl files
+if (MFEM_USE_OCCA)
+  install(DIRECTORY ${MFEM_SOURCE_DIRS}
+    DESTINATION ${INSTALL_INCLUDE_DIR}/mfem
+    FILES_MATCHING PATTERN "*.okl")
+endif()
+
 # Install ${HEADERS}
 # ---
 # foreach (HDR ${HEADERS})

CONTRIBUTING.md

@@ -142,6 +142,16 @@ Origin](#developers-certificate-of-origin-11) at the end of this file.*
   + [`HypreParMatrix`](http://mfem.github.io/doxygen/html/classmfem_1_1HypreParMatrix.html) and [`HypreParVector`](http://mfem.github.io/doxygen/html/classmfem_1_1HypreParVector.html)
   + [`HypreSolver`](http://mfem.github.io/doxygen/html/classmfem_1_1HypreSolver.html) and other [hypre classes](http://mfem.github.io/doxygen/html/hypre_8hpp.html)
 
+- GPU and multi-core CPU support is based on device kernels supporting different
+  backends (CUDA, OCCA, RAJA, OpenMP, etc.) and an internal lightweight
+  device/host memory manager.
+
+- The main device-relevant classes and sources are:
+  + [`Device`](http://mfem.github.io/doxygen/html/device_8hpp.html)
+  + [`MemoryManager`](http://mfem.github.io/doxygen/html/mem_manager_8hpp.html)
+  + the [`MFEM_FORALL`](http://mfem.github.io/doxygen/html/forall_8hpp.html) macro
+  + the [`cuda.hpp`](http://mfem.github.io/doxygen/html/cuda_8hpp.html) and [`occa.hpp`](http://mfem.github.io/doxygen/html/occa_8hpp.html) files
+
 - The `general/` directory contains C++ classes that serve as utilities for
   communication, error handling, arrays, (Boolean) tables, timing, etc.
 

INSTALL

@@ -13,11 +13,17 @@ of MFEM is a (modern) C++ compiler, such as g++. The parallel version of MFEM
 requires an MPI C++ compiler, as well as the following external libraries:
 
 - hypre (a library of high-performance preconditioners)
-  http://www.llnl.gov/CASC/hypre
+  https://github.com/hypre-space/hypre
 
 - METIS (a family of multilevel partitioning algorithms)
   http://glaros.dtc.umn.edu/gkhome/metis/metis/overview
 
+The hypre dependency can be downloaded as a tarball from GitHub or from the
+project webpage https://www.llnl.gov/casc/hypre. For example, the 2.16.0 release
+of hypre is available at
+
+   https://github.com/hypre-space/hypre/archive/v2.16.0.tar.gz
+
 The METIS dependency can be disabled but that is not generally recommended, see
 the option MFEM_USE_METIS.
 
@@ -48,7 +54,7 @@ following package managers:
 
 - Spack, https://github.com/spack/spack
 - OpenHPC, http://openhpc.community
-- Homebrew/Science, https://github.com/Homebrew/homebrew-science
+- Homebrew/Science, https://github.com/Homebrew/homebrew-science (deprecated)
 
 We also recommend downloading and building the MFEM-based GLVis visualization
 tool which can be used to visualize the meshes and solution in MFEM's examples
@@ -60,9 +66,9 @@ Serial build:
    make serial -j 4
 
 Parallel build:
-   (download hypre 2.10.0b and METIS 4 from above URLs)
+   (download hypre and METIS 4 from above URLs)
    (build METIS 4 in ../metis-4.0 relative to mfem/)
-   (build hypre 2.10.0b in ../hypre-2.10.0b relative to mfem/)
+   (build hypre in ../hypre relative to mfem/)
    make parallel -j 4
 
 CUDA build:
@@ -87,13 +93,19 @@ Serial build:
    make -j 4   (assuming "UNIX Makefiles" generator)
 
 Parallel build:
-   (download hypre 2.10.0b and METIS 4 from above URLs)
+   (download hypre and METIS 4 from above URLs)
    (build METIS 4 in ../metis-4.0 relative to mfem/)
-   (build hypre 2.10.0b in ../hypre-2.10.0b relative to mfem/)
+   (build hypre in ../hypre relative to mfem/)
    mkdir <mfem-build-dir> ; cd <mfem-build-dir>
    cmake <mfem-source-dir> -DMFEM_USE_MPI=YES
    make -j 4
 
+CUDA build:
+   (this build requires CMake 3.8 or newer)
+   mkdir <mfem-build-dir> ; cd <mfem-build-dir>
+   cmake <mfem-source-dir> -DMFEM_USE_CUDA=YES
+   make -j 4
+
 Example codes (serial/parallel, depending on the build):
    make examples -j 4
 
@@ -278,6 +290,7 @@ MFEM_THREAD_SAFE = YES/NO
 
 MFEM_USE_LEGACY_OPENMP = YES/NO
    Enable (basic) experimental OpenMP support. Requires MFEM_THREAD_SAFE.
+   This option is deprecated.
 
 MFEM_USE_OPENMP = YES/NO
    Enable the OpenMP backend.
@@ -393,7 +406,8 @@ MFEM_USE_PUMI = YES/NO
 
 MFEM_USE_MM = YES/NO
    Enables support for the MFEM's memory manager (MM), which is required to
-   support devices with different memory spaces.
+   support devices with different memory spaces. This option is required when
+   CUDA support is enabled, i.e. when MFEM_USE_CUDA=YES.
 
 MFEM_USE_CUDA = YES/NO
    Enables support for CUDA devices in MFEM. CUDA is a parallel computing
@@ -406,13 +420,15 @@ MFEM_USE_CUDA = YES/NO
 MFEM_USE_RAJA = YES/NO
    Enable support for the RAJA performance portability layer in MFEM. RAJA
    provides a portable abstraction for loops, supporting different programming
-   model backends. When using the RAJA CUDA backend, MFEM_USE_MM is required.
+   model backends. When using RAJA built with CUDA support, CUDA support must be
+   also enabled in MFEM, i.e. MFEM_USE_CUDA=YES must be set.
 
 MFEM_USE_OCCA = YES/NO
    Enables support for the OCCA library in MFEM. OCCA is an open-source library
    which aims to make it easy to program different types of devices (e.g. CPU,
    GPU, FPGA) by providing an unified API for interacting with JIT-compiled
-   backends. When using the OCCA CUDA backend, MFEM_USE_MM is required.
+   backends. In order to use the OCCA CUDA backend, CUDA support must be enabled
+   in MFEM as well, i.e. MFEM_USE_CUDA=YES must be set.
 
 MFEM_BUILD_TAG = (any value)
    An optional tag to characterize the build. Exported to config/config.mk.
@@ -435,7 +451,7 @@ directory and use the string @MFEM_DIR@, e.g. HYPRE_OPT = -I@MFEM_DIR@/../hypre.
 The specific libraries and their options are:
 
 - HYPRE, required for the parallel build, i.e. when MFEM_USE_MPI = YES.
-  URL: http://www.llnl.gov/CASC/hypre
+  URL: https://github.com/hypre-space/hypre and https://www.llnl.gov/casc/hypre
   Options: HYPRE_OPT, HYPRE_LIB.
 
 - METIS, used when MFEM_USE_METIS = YES. If using METIS 5, set
@@ -645,6 +661,8 @@ Configuration variables (CMake)
 ===============================
 See the configuration file config/defaults.cmake for the default settings.
 
+Note: the option MFEM_USE_CUDA requires CMake version 3.8 or newer!
+
 Non-standard CMake variables for compilers:
    CXX    - If set, overwrite the auto-detected C++ compiler, serial build
    MPICXX - If set, overwrite the auto-detected MPI C++ compiler, parallel build
@@ -675,9 +693,6 @@ MFEM_USE_NETCDF
 MFEM_USE_MPFR
 MFEM_USE_GZSTREAM
 MFEM_USE_PUMI
-
-The following GNU make options are not supported with CMake yet:
-
 MFEM_USE_CUDA
 MFEM_USE_OCCA
 MFEM_USE_RAJA
@@ -728,6 +743,8 @@ The CMake build system adds auto-detection for the following packages/libraries:
  - LIBUNWIND
  - POSIXCLOCKS
  - PUMI
+ - OCCA
+ - RAJA
 
 The following built-in CMake packages are also used:
 

README

@@ -8,9 +8,9 @@
                                http://mfem.org
 
 MFEM is a modular parallel C++ library for finite element methods. Its goal is
-to enable the research and development of scalable finite element discretization
-and solver algorithms through general finite element abstractions, accurate and
-flexible visualization, and tight integration with the hypre library.
+to enable high-performance scalable finite element discretization research and
+application development on a wide variety of platforms, ranging from laptops to
+supercomputers.
 
 * For building instructions, see the file INSTALL, or type "make help".
 
@@ -39,23 +39,24 @@ conforming and non-conforming (AMR) adaptive refinement. Arbitrary element
 transformations, allowing for high-order mesh elements with curved boundaries,
 are also supported.
 
-MFEM is commonly used as a "finite element to linear algebra translator", since
-it can take a problem described in terms of finite element-type objects, and
-produce the corresponding linear algebra vectors and sparse matrices. In order
-to facilitate this, MFEM uses compressed sparse row (CSR) sparse matrix storage
-and includes simple smoothers and Krylov solvers, such as PCG, MINRES and GMRES,
-as well as support for sequential sparse direct solvers from the SuiteSparse
+When used as a "finite element to linear algebra translator", MFEM can take a
+problem described in terms of finite element-type objects, and produce the
+corresponding linear algebra vectors and fully or partially assembled operators,
+e.g. in the form of global sparse matrices or matrix-free operators. The library
+includes simple smoothers and Krylov solvers, such as PCG, MINRES and GMRES, as
+well as support for sequential sparse direct solvers from the SuiteSparse
 library. Nonlinear solvers (the Newton method), eigensolvers (LOBPCG), and
 several explicit and implicit Runge-Kutta time integrators are also available.
 
 MFEM supports MPI-based parallelism throughout the library, and can readily be
-used as a scalable unstructured finite element problem generator. MFEM-based
-applications require minimal changes to transition from a serial to a
-high-performing parallel version of the code, where they can take advantage of
-the integrated scalable linear solvers from the hypre library. Comprehensive
-support for other external packages, e.g. PETSc and SUNDIALS is also included,
-giving access to many additional linear and nonlinear solvers, preconditioners,
-time integrators, etc.
+used as a scalable unstructured finite element problem generator. As of version
+4.0, MFEM offers initial support for GPU acceleration, and programming models,
+such as CUDA, OCCA, RAJA and OpenMP. MFEM-based applications require minimal
+changes to switch from a serial to a high-performing MPI-parallel version of the
+code, where they can take advantage of the integrated linear solvers from the
+hypre library.  Comprehensive support for other external packages, e.g. PETSc
+and SUNDIALS is also included, giving access to many additional linear and
+nonlinear solvers, preconditioners, time integrators, etc.
 
 For examples of using MFEM, see the examples/ and miniapps/ directories, as well
 as the OpenGL visualization tool GLVis which is available at http://glvis.org.

config/XSDKDefaults.cmake

@@ -74,7 +74,7 @@
 
 IF (NOT COMMAND PRINT_VAR)
   FUNCTION(PRINT_VAR  VAR_NAME)
-    MESSAGE("-- " "${VAR_NAME} = '${${VAR_NAME}}'")
+    MESSAGE(STATUS "${VAR_NAME} = '${${VAR_NAME}}'")
   ENDFUNCTION()
 ENDIF()
 
@@ -166,14 +166,14 @@ IF (USE_XSDK_DEFAULTS)
     ENDIF()
     XSDK_HANDLE_LANG_DEFAULTS(Fortran  FC  "FFLAGS;FCFLAGS")
   ENDIF()
-  
+
   # Set XSDK defaults for other CMake variables
-  
+
   IF ("${BUILD_SHARED_LIBS}"  STREQUAL  "")
     MESSAGE("-- " "XSDK: Setting default BUILD_SHARED_LIBS=TRUE")
     SET(BUILD_SHARED_LIBS  TRUE  CACHE  BOOL  "Set by default in XSDK mode")
   ENDIF()
-  
+
   IF ("${CMAKE_BUILD_TYPE}"  STREQUAL  "")
     MESSAGE("-- " "XSDK: Setting default CMAKE_BUILD_TYPE=DEBUG")
     SET(CMAKE_BUILD_TYPE  DEBUG  CACHE  STRING  "Set by default in XSDK mode")
@@ -181,6 +181,13 @@ IF (USE_XSDK_DEFAULTS)
 
 ENDIF()
 
+
+##################################################################################
+#
+#               MFEM-specific additions: set TPL MFEM_USE_* defaults
+#
+##################################################################################
+
 IF (DEFINED TPL_ENABLE_MPI)
   SET(MFEM_USE_MPI ${TPL_ENABLE_MPI} CACHE BOOL "Enable MPI parallel build" FORCE)
 ENDIF()
@@ -252,3 +259,15 @@ ENDIF()
 IF (DEFINED TPL_ENABLE_PUMI)
   SET(MFEM_USE_PUMI ${TPL_ENABLE_PUMI} CACHE BOOL "Enable PUMI" FORCE)
 ENDIF()
+
+IF (DEFINED TPL_ENABLE_CUDA)
+  SET(MFEM_USE_CUDA ${TPL_ENABLE_CUDA} CACHE BOOL "Enable CUDA" FORCE)
+ENDIF()
+
+IF (DEFINED TPL_ENABLE_OCCA)
+  SET(MFEM_USE_OCCA ${TPL_ENABLE_OCCA} CACHE BOOL "Enable OCCA" FORCE)
+ENDIF()
+
+IF (DEFINED TPL_ENABLE_RAJA)
+  SET(MFEM_USE_RAJA ${TPL_ENABLE_RAJA} CACHE BOOL "Enable RAJA" FORCE)
+ENDIF()

config/cmake/MFEMConfig.cmake.in

@@ -41,6 +41,10 @@ set(MFEM_USE_MPFR @MFEM_USE_MPFR@)
 set(MFEM_USE_SIDRE @MFEM_USE_SIDRE@)
 set(MFEM_USE_CONDUIT @MFEM_USE_CONDUIT@)
 set(MFEM_USE_PUMI @MFEM_USE_PUMI@)
+set(MFEM_USE_MM @MFEM_USE_MM@)
+set(MFEM_USE_CUDA @MFEM_USE_CUDA@)
+set(MFEM_USE_OCCA @MFEM_USE_OCCA@)
+set(MFEM_USE_RAJA @MFEM_USE_RAJA@)
 
 set(MFEM_CXX_COMPILER "@CMAKE_CXX_COMPILER@")
 set(MFEM_CXX_FLAGS "@CMAKE_CXX_FLAGS@")

config/cmake/config.hpp.in

@@ -30,6 +30,12 @@
 #define MFEM_VERSION_MINOR (((MFEM_VERSION)/100)%100)
 #define MFEM_VERSION_PATCH ((MFEM_VERSION)%100)
 
+// MFEM source directory.
+#define MFEM_SOURCE_DIR "@MFEM_SOURCE_DIR@"
+
+// MFEM install directory.
+#define MFEM_INSTALL_DIR "@MFEM_INSTALL_DIR@"
+
 // Description of the git commit used to build MFEM.
 #cmakedefine MFEM_GIT_STRING "@MFEM_GIT_STRING@"
 
@@ -104,6 +110,19 @@
 // Enable MFEM functionality based on the PUMI library
 #cmakedefine MFEM_USE_PUMI
 
+// Build the GPU/CUDA-enabled version of the MFEM library.
+// Requires a CUDA compiler (nvcc).
+#cmakedefine MFEM_USE_CUDA
+
+// Enable MFEM functionality based on the RAJA library
+#cmakedefine MFEM_USE_RAJA
+
+// Enable MFEM functionality based on the OCCA library
+#cmakedefine MFEM_USE_OCCA
+
+// Enable MFEM's internal Memory Manager (needed e.g. for MFEM_USE_CUDA)
+#cmakedefine MFEM_USE_MM
+
 // Which library functions to use in class StopWatch for measuring time.
 // For a list of the available options, see INSTALL.
 // If not defined, an option is selected automatically.

config/cmake/modules/FindOCCA.cmake

@@ -0,0 +1,19 @@
+# Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the
+# Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability see http://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the GNU Lesser General Public License (as published by the Free
+# Software Foundation) version 2.1 dated February 1999.
+
+# Defines the following variables:
+#   - OCCA_FOUND
+#   - OCCA_LIBRARIES
+#   - OCCA_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+mfem_find_package(OCCA OCCA OCCA_DIR "include" "occa.hpp" "lib" "occa"
+  "Paths to headers required by OCCA." "Libraries required by OCCA.")

config/cmake/modules/FindRAJA.cmake

@@ -0,0 +1,30 @@
+# Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the
+# Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability see http://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the GNU Lesser General Public License (as published by the Free
+# Software Foundation) version 2.1 dated February 1999.
+
+# Defines the following variables:
+#   - RAJA_FOUND
+#   - RAJA_LIBRARIES
+#   - RAJA_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+mfem_find_package(RAJA RAJA RAJA_DIR "include" "RAJA/RAJA.hpp" "lib" "RAJA"
+  "Paths to headers required by RAJA." "Libraries required by RAJA.")
+
+if (NOT RAJA_CONFIG_CMAKE)
+   set(RAJA_CONFIG_CMAKE "${RAJA_DIR}/share/raja/cmake/raja-config.cmake")
+endif()
+if (EXISTS "${RAJA_CONFIG_CMAKE}")
+   include("${RAJA_CONFIG_CMAKE}")
+   if (ENABLE_CUDA AND NOT MFEM_USE_CUDA)
+      message(FATAL_ERROR
+         "RAJA is built with CUDA: MFEM_USE_CUDA=YES is required")
+   endif()
+endif()

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -232,10 +232,12 @@ function(mfem_find_package Name Prefix DirVar IncSuffixes Header LibSuffixes
   # If we have the TPL_ versions of _INCLUDE_DIRS and _LIBRARIES then set the
   # standard ${Prefix} versions
   if (TPL_${Prefix}_INCLUDE_DIRS)
-    set(${Prefix}_INCLUDE_DIRS ${TPL_${Prefix}_INCLUDE_DIRS} CACHE STRING "TPL_${Prefix}_INCLUDE_DIRS was found." FORCE)
+    set(${Prefix}_INCLUDE_DIRS ${TPL_${Prefix}_INCLUDE_DIRS} CACHE STRING
+      "TPL_${Prefix}_INCLUDE_DIRS was found." FORCE)
   endif()
   if (TPL_${Prefix}_LIBRARIES)
-    set(${Prefix}_LIBRARIES ${TPL_${Prefix}_LIBRARIES} CACHE STRING "TPL_${Prefix}_LIBRARIES was found." FORCE)
+    set(${Prefix}_LIBRARIES ${TPL_${Prefix}_LIBRARIES} CACHE STRING
+      "TPL_${Prefix}_LIBRARIES was found." FORCE)
   endif()
 
   # Quick return
@@ -718,7 +720,8 @@ function(mfem_export_mk_files)
       MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE
       MFEM_USE_SUPERLU MFEM_USE_STRUMPACK MFEM_USE_GECKO MFEM_USE_GNUTLS
       MFEM_USE_NETCDF MFEM_USE_PETSC MFEM_USE_MPFR MFEM_USE_SIDRE
-      MFEM_USE_CONDUIT MFEM_USE_PUMI)
+      MFEM_USE_CONDUIT MFEM_USE_PUMI MFEM_USE_MM MFEM_USE_CUDA MFEM_USE_OCCA
+      MFEM_USE_RAJA)
   foreach(var ${CONFIG_MK_BOOL_VARS})
     if (${var})
       set(${var} YES)
@@ -726,6 +729,7 @@ function(mfem_export_mk_files)
       set(${var} NO)
     endif()
   endforeach()
+  # TODO: Add support for MFEM_USE_CUDA=YES
   set(MFEM_CXX ${CMAKE_CXX_COMPILER})
   set(MFEM_CPPFLAGS "")
   string(STRIP "${CMAKE_CXX_FLAGS_${BUILD_TYPE}} ${CMAKE_CXX_FLAGS}"

config/config.hpp

@@ -56,4 +56,9 @@
 #endif
 #endif // MFEM_USE_MPI not defined
 
+// CUDA requires the memory manager
+#if defined(MFEM_USE_CUDA) && !defined(MFEM_USE_MM)
+#error Building with CUDA (MFEM_USE_CUDA=YES) requires MFEM_USE_MM=YES
+#endif
+
 #endif // MFEM_CONFIG_HPP

config/defaults.cmake

@@ -42,6 +42,10 @@ option(MFEM_USE_MPFR "Enable MPFR usage." OFF)
 option(MFEM_USE_SIDRE "Enable Axom/Sidre usage" OFF)
 option(MFEM_USE_CONDUIT "Enable Conduit usage" OFF)
 option(MFEM_USE_PUMI "Enable PUMI" OFF)
+option(MFEM_USE_MM "Enable MFEM's memory manager" OFF)
+option(MFEM_USE_CUDA "Enable CUDA" OFF)
+option(MFEM_USE_OCCA "Enable OCCA" OFF)
+option(MFEM_USE_RAJA "Enable RAJA" OFF)
 
 set(MFEM_MPI_NP 4 CACHE STRING "Number of processes used for MPI tests")
 
@@ -59,13 +63,16 @@ option(MFEM_ENABLE_MINIAPPS "Build all of the miniapps" OFF)
 # set(CXX g++)
 # set(MPICXX mpicxx)
 
+# Set the target CUDA architecture
+set(CUDA_ARCH "sm_60" CACHE STRING "Target CUDA architecture.")
+
 set(MFEM_DIR ${CMAKE_CURRENT_SOURCE_DIR})
 
 # The *_DIR paths below will be the first place searched for the corresponding
 # headers and library. If these fail, then standard cmake search is performed.
 # Note: if the variables are already in the cache, they are not overwritten.
 
-set(HYPRE_DIR "${MFEM_DIR}/../hypre-2.10.0b/src/hypre" CACHE PATH
+set(HYPRE_DIR "${MFEM_DIR}/../hypre/src/hypre" CACHE PATH
     "Path to the hypre library.")
 # If hypre was compiled to depend on BLAS and LAPACK:
 # set(HYPRE_REQUIRED_PACKAGES "BLAS" "LAPACK" CACHE STRING
@@ -154,6 +161,9 @@ set(Axom_REQUIRED_PACKAGES "Conduit/relay" CACHE STRING
 set(PUMI_DIR "${MFEM_DIR}/../pumi-2.1.0" CACHE STRING
     "Directory where PUMI is installed")
 
+set(OCCA_DIR "${MFEM_DIR}/../occa" CACHE PATH "Path to OCCA")
+set(RAJA_DIR "${MFEM_DIR}/../raja" CACHE PATH "Path to RAJA")
+
 set(BLAS_INCLUDE_DIRS "" CACHE STRING "Path to BLAS headers.")
 set(BLAS_LIBRARIES "" CACHE STRING "The BLAS library.")
 set(LAPACK_INCLUDE_DIRS "" CACHE STRING "Path to LAPACK headers.")

config/defaults.mk

@@ -136,7 +136,7 @@ LIBUNWIND_OPT = -g
 LIBUNWIND_LIB = $(if $(NOTMAC),-lunwind -ldl,)
 
 # HYPRE library configuration (needed to build the parallel version)
-HYPRE_DIR = @MFEM_DIR@/../hypre-2.10.0b/src/hypre
+HYPRE_DIR = @MFEM_DIR@/../hypre/src/hypre
 HYPRE_OPT = -I$(HYPRE_DIR)/include
 HYPRE_LIB = -L$(HYPRE_DIR)/lib -lHYPRE
 
@@ -300,19 +300,17 @@ PUMI_OPT = -I$(PUMI_DIR)/include
 PUMI_LIB = -L$(PUMI_DIR)/lib -lpumi -lcrv -lma -lmds -lapf -lpcu -lgmi -lparma\
    -llion -lmth -lapf_zoltan -lspr
 
-# CUDA library configuration. Since we compile and link with nvcc (when CUDA is
-# enabled) we only need to explicitly link with the CUDA driver, libcuda.*,
-# which is usually in a system path.
+# CUDA library configuration (currently not needed)
 CUDA_OPT =
-CUDA_LIB = $(if $(NOTMAC),,-L/usr/local/cuda/lib) -lcuda
+CUDA_LIB =
 
 # OCCA library configuration
-OCCA_DIR ?= @MFEM_DIR@/../occa
+OCCA_DIR = @MFEM_DIR@/../occa
 OCCA_OPT = -I$(OCCA_DIR)/include
 OCCA_LIB = $(XLINKER)-rpath,$(OCCA_DIR)/lib -L$(OCCA_DIR)/lib -locca
 
 # RAJA library configuration
-RAJA_DIR ?= @MFEM_DIR@/../raja
+RAJA_DIR = @MFEM_DIR@/../raja
 RAJA_OPT = -I$(RAJA_DIR)/include
 ifdef CUB_DIR
    RAJA_OPT += -I$(CUB_DIR)

config/sample-runs.sh

@@ -18,6 +18,8 @@ run_prefix=""
 run_vg="valgrind --leak-check=full --show-reachable=yes --track-origins=yes"
 run_suffix="-no-vis"
 skip_gen_meshes="yes"
+# filter-out device runs ("no") or non-device runs ("yes"):
+device_runs="no"
 cur_dir="${PWD}"
 mfem_dir="$(cd "$(dirname "$0")"/.. && pwd)"
 mfem_build_dir=""
@@ -148,6 +150,11 @@ function extract_sample_runs()
    if [ "$skip_gen_meshes" == "yes" ]; then
       runs=`printf "%s" "$runs" | grep -v ".* -m .*\.gen"`
    fi
+   if [ "$device_runs" == "yes" ]; then
+      runs=`printf "%s" "$runs" | grep ".* -d .*"`
+   else
+      runs=`printf "%s" "$runs" | grep -v ".* -d .*"`
+   fi
    IFS=$'\n'
    runs=(${runs})
    IFS="${old_IFS}"
@@ -169,6 +176,9 @@ function help_message()
       -g <dir> <pattern>
                   Specify explicitly a group (dir + file pattern) to run; This
                   option can be used multiple times to define multiple groups
+      -dev        configure only sample runs using devices.
+                  To test with a parallel build, the parallel (-p|-par) option
+                  should be set first on the command line.
       -v          Enable valgrind
       -o <dir>    [${output_dir:-"<empty>: output goes to stdout"}]
                   If not empty, save output to files inside <dir>
@@ -253,7 +263,7 @@ case "$1" in
    -h|-help)
       opt_help="yes"
       ;;
-   -p|-parallel)
+   -p|-par)
       mfem_config="MFEM_USE_MPI=YES MFEM_DEBUG=NO"
       ;;
    -g)
@@ -264,6 +274,11 @@ case "$1" in
       groups=("${groups[@]}" "${test_group}")
       shift 2
       ;;
+   -dev)
+       device_runs="yes"
+       mfem_config+=" MFEM_USE_CUDA=YES MFEM_USE_MM=YES \
+MFEM_USE_OCCA=YES MFEM_USE_RAJA=YES MFEM_USE_OPENMP=YES"
+      ;;
    -v)
       valgrind="yes"
       ;;
@@ -294,6 +309,10 @@ case "$1" in
    -n)
       run_prefix="echo"
       ;;
+   -*)
+      echo "unknown option: '$1'"
+      exit 1
+      ;;
    *=*)
       eval $1
       ;;

doc/CodeDocumentation.dox

@@ -35,6 +35,10 @@ namespace mfem {
  * - HypreParMatrix and HypreParVector
  * - HypreSolver and other \link hypre.hpp hypre classes\endlink
  *
+ * <H3>Main GPU classes</H3>
+ * - Device
+ * - MemoryManager
+ *
  * <H3>Example codes</H3>
  * - <a class="el" href="examples_2ex1_8cpp_source.html">Example 1</a>: nodal H1 FEM for the Laplace problem
  * - <a class="el" href="examples_2ex1p_8cpp_source.html">Example 1p</a>: parallel nodal H1 FEM for the Laplace problem

examples/CMakeLists.txt

@@ -27,8 +27,8 @@ list(APPEND ALL_EXE_SRCS
   ex18.cpp
   ex19.cpp
   ex20.cpp
-  ex21.cpp
   ex22.cpp
+  ex23.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -53,8 +53,8 @@ if (MFEM_USE_MPI)
     ex18p.cpp
     ex19p.cpp
     ex20p.cpp
-    ex21p.cpp
     ex22p.cpp
+    ex23p.cpp
     )
 endif()
 

examples/ex1.cpp

@@ -26,12 +26,12 @@
 //               ex1 -m ../data/mobius-strip.mesh -o -1 -sc
 //
 // Device sample runs:
-//             > ex1 -pa -d cuda
-//             > ex1 -pa -d raja-cuda
-//             > ex1 -pa -d occa-cuda
-//             > ex1 -pa -d raja-omp
-//             > ex1 -pa -d occa-omp
-//             > ex1 -m ../data/beam-hex.mesh -pa -d cuda
+//               ex1 -pa -d cuda
+//               ex1 -pa -d raja-cuda
+//               ex1 -pa -d occa-cuda
+//               ex1 -pa -d raja-omp
+//               ex1 -pa -d occa-omp
+//               ex1 -m ../data/beam-hex.mesh -pa -d cuda
 //
 // Description:  This example code demonstrates the use of MFEM to define a
 //               simple finite element discretization of the Laplace problem

examples/ex1p.cpp

@@ -26,9 +26,9 @@
 //               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
 //
 // Device sample runs:
-//             > mpirun -np 4 ex1p -pa -d cuda
-//             > mpirun -np 4 ex1p -pa -d occa-cuda
-//             > mpirun -np 4 ex1p -pa -d raja-omp
+//               mpirun -np 4 ex1p -pa -d cuda
+//               mpirun -np 4 ex1p -pa -d occa-cuda
+//               mpirun -np 4 ex1p -pa -d raja-omp
 //
 // Description:  This example code demonstrates the use of MFEM to define a
 //               simple finite element discretization of the Laplace problem

examples/ex21.cpp

@@ -1,477 +0,0 @@
-//                               MFEM Example 21
-//
-// Compile with: make ex21
-//
-// Sample runs:  ex21 -m ../data/inline-segment.mesh -o 3
-//               ex21 -m ../data/inline-tri.mesh -o 3
-//               ex21 -m ../data/inline-quad.mesh -o 3
-//               ex21 -m ../data/inline-quad.mesh -o 3 -p 1
-//               ex21 -m ../data/inline-quad.mesh -o 3 -p 2
-//               ex21 -m ../data/inline-tet.mesh -o 2
-//               ex21 -m ../data/inline-hex.mesh -o 2
-//               ex21 -m ../data/inline-hex.mesh -o 2 -p 1
-//               ex21 -m ../data/inline-hex.mesh -o 2 -p 2
-//               ex21 -m ../data/star.mesh -o 2 -sigma 10.0
-//
-// Description:  This example code demonstrates the use of MFEM to define and
-//               solve simple complex-valued linear systems.  We implement three
-//               variants of a damped harmonic oscillator:
-//
-//               1) A scalar H1 field
-//                  -Div(a Grad u) - omega^2 b u + i omega c u = 0
-//
-//               2) A vector H(Curl) field
-//                  Curl(a Curl u) - omega^2 b u + i omega c u = 0
-//
-//               3) A vector H(Div) field
-//                  -Grad(a Div u) - omega^2 b u + i omega c u = 0
-//
-//               In each case the field is driven by a forced oscillation, with
-//               angular frequency omega, imposed at the boundary or a portion
-//               of the boundary.
-//
-//               In electromagnetics the coefficients are typically named the
-//               permeability, mu = 1/a, permittivity, epsilon = b, and
-//               conductivity, sigma = c.  The user can specify these constants
-//               using either set of names.
-//
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-static double mu_ = 1.0;
-static double epsilon_ = 1.0;
-static double sigma_ = 20.0;
-static double omega_ = 10.0;
-
-double u0_real_exact(const Vector &);
-double u0_imag_exact(const Vector &);
-
-void u1_real_exact(const Vector &, Vector &);
-void u1_imag_exact(const Vector &, Vector &);
-
-void u2_real_exact(const Vector &, Vector &);
-void u2_imag_exact(const Vector &, Vector &);
-
-bool check_for_inline_mesh(const char * mesh_file);
-
-int main(int argc, char *argv[])
-{
-   // 1. Parse command-line options.
-   const char *mesh_file = "../data/inline-quad.mesh";
-   int ref_levels = 0;
-   int order = 1;
-   int prob = 0;
-   double freq = -1.0;
-   double a_coef = 0.0;
-   bool visualization = 1;
-   bool herm_conv = true;
-   bool exact_sol = true;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ref_levels, "-r", "--refine",
-                  "Number of times to refine the mesh uniformly.");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree).");
-   args.AddOption(&prob, "-p", "--problem-type",
-                  "Choose from 0: H_1, 1: H(Curl), or 2: H(Div) "
-                  "damped harmonic oscillator.");
-   args.AddOption(&a_coef, "-a", "--stiffness-coef",
-                  "Stiffness coefficient (spring constant or 1/mu).");
-   args.AddOption(&epsilon_, "-b", "--mass-coef",
-                  "Mass coefficient (or epsilon).");
-   args.AddOption(&sigma_, "-c", "--damping-coef",
-                  "Damping coefficient (or sigma).");
-   args.AddOption(&mu_, "-mu", "--permeability",
-                  "Permeability of free space (or 1/(spring constant)).");
-   args.AddOption(&epsilon_, "-eps", "--permittivity",
-                  "Permittivity of free space (or mass constant).");
-   args.AddOption(&sigma_, "-sigma", "--conductivity",
-                  "Conductivity (or damping constant).");
-   args.AddOption(&freq, "-f", "--frequency",
-                  "Frequency (in Hz).");
-   args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
-                  "--no-hermitian", "Use convention for Hermitian operators.");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(cout);
-      return 1;
-   }
-   args.PrintOptions(cout);
-
-   if ( a_coef != 0.0 )
-   {
-      mu_ = 1.0 / a_coef;
-   }
-   if ( freq > 0.0 )
-   {
-      omega_ = 2.0 * M_PI * freq;
-   }
-
-   exact_sol = check_for_inline_mesh(mesh_file);
-   if (exact_sol)
-   {
-      cout << "Identified an 'inline' mesh" << endl;
-   }
-
-   ComplexOperator::Convention conv =
-      herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
-
-   // 2. Read the mesh from the given mesh file. We can handle triangular,
-   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes
-   //    with the same code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   int dim = mesh->Dimension();
-
-   // 3. Refine the mesh to increase resolution. In this example we do
-   //    'ref_levels' of uniform refinement where the user specifies
-   //    the number of levels with the '-r' option.
-   for (int l = 0; l < ref_levels; l++)
-   {
-      mesh->UniformRefinement();
-   }
-
-   // 4. Define a finite element space on the mesh. Here we use continuous
-   //    Lagrange, Nedelec, or Raviart-Thomas finite elements of the specified
-   //    order.
-   if (dim == 1 && prob != 0 )
-   {
-      cout << "Switching to problem type 0, H1 basis functions, "
-           << "for 1 dimensional mesh." << endl;
-      prob = 0;
-   }
-
-   FiniteElementCollection *fec;
-   switch (prob)
-   {
-      case 0:  fec = new H1_FECollection(order, dim);      break;
-      case 1:  fec = new ND_FECollection(order, dim);      break;
-      case 2:  fec = new RT_FECollection(order - 1, dim);  break;
-   }
-   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
-   cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
-        << endl;
-
-   // 5. Determine the list of true (i.e. conforming) essential boundary dofs.
-   //    In this example, the boundary conditions are defined based on the type
-   //    of mesh and the problem type.
-   Array<int> ess_tdof_list;
-   Array<int> ess_bdr;
-   if (mesh->bdr_attributes.Size())
-   {
-      ess_bdr.SetSize(mesh->bdr_attributes.Max());
-      ess_bdr = 1;
-      if (exact_sol)
-      {
-         switch (prob)
-         {
-            case 0:   ess_bdr = 0; ess_bdr[0] = 1;  break;
-            default:  ess_bdr = 1; ess_bdr[2] = 0;  break;
-         }
-      }
-      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-   }
-
-   // 6. Set up the linear form b(.) which corresponds to the
-   //    right-hand side of the FEM linear system.
-   ComplexLinearForm b(fespace, conv);
-   b.Vector::operator=(0.0);
-
-   // 7. Define the solution vector u as a finite element grid function
-   //    corresponding to fespace. Initialize u with initial guess of 1+0i
-   //    or the exact solution if it is known.
-   ComplexGridFunction u(fespace);
-   ComplexGridFunction * u_exact = NULL;
-   if (exact_sol) { u_exact = new ComplexGridFunction(fespace); }
-
-   FunctionCoefficient u0_r(u0_real_exact);
-   FunctionCoefficient u0_i(u0_imag_exact);
-   VectorFunctionCoefficient u1_r(dim, u1_real_exact);
-   VectorFunctionCoefficient u1_i(dim, u1_imag_exact);
-   VectorFunctionCoefficient u2_r(dim, u2_real_exact);
-   VectorFunctionCoefficient u2_i(dim, u2_imag_exact);
-
-   ConstantCoefficient zeroCoef(0.0);
-   ConstantCoefficient oneCoef(1.0);
-
-   Vector zeroVec(dim); zeroVec = 0.0;
-   Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
-   VectorConstantCoefficient zeroVecCoef(zeroVec);
-   VectorConstantCoefficient oneVecCoef(oneVec);
-
-   switch (prob)
-   {
-      case 0:
-         u.ProjectBdrCoefficient(oneCoef, zeroCoef, ess_bdr);
-         if (exact_sol) { u_exact->ProjectCoefficient(u0_r, u0_i); }
-         break;
-      case 1:
-         u.ProjectBdrCoefficientTangent(oneVecCoef, zeroVecCoef, ess_bdr);
-         if (exact_sol) { u_exact->ProjectCoefficient(u1_r, u1_i); }
-         break;
-      case 2:
-         u.ProjectBdrCoefficientNormal(oneVecCoef, zeroVecCoef, ess_bdr);
-         if (exact_sol) { u_exact->ProjectCoefficient(u2_r, u2_i); }
-         break;
-   }
-
-   if (visualization && exact_sol)
-   {
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      socketstream sol_sock_r(vishost, visport);
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_r.precision(8);
-      sol_sock_i.precision(8);
-      sol_sock_r << "solution\n" << *mesh << u_exact->real()
-                 << "window_title 'Exact Real Part'" << flush;
-      sol_sock_i << "solution\n" << *mesh << u_exact->imag()
-                 << "window_title 'Exact Imaginary Part'" << flush;
-   }
-
-   // 8. Set up the sesquilinear form a(.,.) on the finite element
-   //    space corresponding to the damped harmonic oscillator operator
-   //    of the appropriate type:
-   //
-   //    0) A scalar H1 field
-   //       -Div(a Grad) - omega^2 b + i omega c
-   //
-   //    1) A vector H(Curl) field
-   //       Curl(a Curl) - omega^2 b + i omega c
-   //
-   //    2) A vector H(Div) field
-   //       -Grad(a Div) - omega^2 b + i omega c
-   //
-   ConstantCoefficient stiffnessCoef(1.0/mu_);
-   ConstantCoefficient massCoef(-omega_ * omega_ * epsilon_);
-   ConstantCoefficient lossCoef(omega_ * sigma_);
-   ConstantCoefficient negMassCoef(omega_ * omega_ * epsilon_);
-
-   SesquilinearForm *a = new SesquilinearForm(fespace, conv);
-   switch (prob)
-   {
-      case 0:
-         a->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new MassIntegrator(massCoef),
-                                new MassIntegrator(lossCoef));
-         break;
-      case 1:
-         a->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
-                                new VectorFEMassIntegrator(lossCoef));
-         break;
-      case 2:
-         a->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
-                                new VectorFEMassIntegrator(lossCoef));
-         break;
-   }
-
-   // 9. Assemble the bilinear form and the corresponding linear
-   //    system, applying any necessary transformations such as:
-   //    assembly, eliminating boundary conditions, applying conforming
-   //    constraints for non-conforming AMR, etc.
-   a->Assemble();
-
-   OperatorHandle A;
-   Vector B, U;
-
-   a->FormLinearSystem(ess_tdof_list, u, b, A, U, B);
-   u = 0.0;
-   U = 0.0;
-
-   {
-      ComplexSparseMatrix * Asp =
-         dynamic_cast<ComplexSparseMatrix*>(A.Ptr());
-
-      cout << "Size of linear system: "
-           << 2 * Asp->real().Width() << endl << endl;
-   }
-
-   // 10. Define and apply a GMRES solver for AU=B.
-   {
-      GMRESSolver gmres;
-      gmres.SetOperator(*A.Ptr());
-      gmres.SetRelTol(1e-12);
-      gmres.SetMaxIter(1000);
-      gmres.SetPrintLevel(1);
-      gmres.Mult(B, U);
-   }
-
-   // 11. Recover the solution as a finite element grid function and
-   //     compute the errors if the exact solution is known.
-   a->RecoverFEMSolution(U, b, u);
-
-   if (exact_sol)
-   {
-      double err_r = -1.0;
-      double err_i = -1.0;
-
-      switch (prob)
-      {
-         case 0:
-            err_r = u.real().ComputeL2Error(u0_r);
-            err_i = u.imag().ComputeL2Error(u0_i);
-            break;
-         case 1:
-            err_r = u.real().ComputeL2Error(u1_r);
-            err_i = u.imag().ComputeL2Error(u1_i);
-            break;
-         case 2:
-            err_r = u.real().ComputeL2Error(u2_r);
-            err_i = u.imag().ComputeL2Error(u2_i);
-            break;
-      }
-
-      cout << endl;
-      cout << "|| Re (u_h - u) ||_{L^2} = " << err_r << endl;
-      cout << "|| Im (u_h - u) ||_{L^2} = " << err_i << endl;
-      cout << endl;
-   }
-
-   // 12. Save the refined mesh and the solution. This output can be
-   //     viewed later using GLVis: "glvis -m mesh -g sol".
-   {
-      ofstream mesh_ofs("refined.mesh");
-      mesh_ofs.precision(8);
-      mesh->Print(mesh_ofs);
-
-      ofstream sol_r_ofs("sol_r.gf");
-      ofstream sol_i_ofs("sol_i.gf");
-      sol_r_ofs.precision(8);
-      sol_i_ofs.precision(8);
-      u.real().Save(sol_r_ofs);
-      u.imag().Save(sol_i_ofs);
-   }
-
-   // 13. Send the solution by socket to a GLVis server.
-   if (visualization)
-   {
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      socketstream sol_sock_r(vishost, visport);
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_r.precision(8);
-      sol_sock_i.precision(8);
-      sol_sock_r << "solution\n" << *mesh << u.real()
-                 << "window_title 'Comp Real Part'" << flush;
-      sol_sock_i << "solution\n" << *mesh << u.imag()
-                 << "window_title 'Comp Imaginary Part'" << flush;
-   }
-   if (visualization && exact_sol)
-   {
-      *u_exact -= u;
-
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      socketstream sol_sock_r(vishost, visport);
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_r.precision(8);
-      sol_sock_i.precision(8);
-      sol_sock_r << "solution\n" << *mesh << u_exact->real()
-                 << "window_title 'Exact-Comp Real Part'" << flush;
-      sol_sock_i << "solution\n" << *mesh << u_exact->imag()
-                 << "window_title 'Exact-Comp Imaginary Part'" << flush;
-   }
-   if (visualization)
-   {
-      GridFunction u_t(fespace);
-      u_t = u.real();
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      socketstream sol_sock(vishost, visport);
-      sol_sock.precision(8);
-      sol_sock << "solution\n" << *mesh << u_t
-               << "window_title 'Harmonic Solution (t = 0.0 T)'"
-               << "pause\n" << flush;
-
-      cout << "GLVis visualization paused."
-           << " Press space (in the GLVis window) to resume it.\n";
-      int num_frames = 32;
-      int i = 0;
-      while (sol_sock)
-      {
-         double t = (double)(i % num_frames) / num_frames;
-         ostringstream oss;
-         oss << "Harmonic Solution (t = " << t << " T)";
-
-         add(cos( 2.0 * M_PI * t), u.real(),
-             sin(-2.0 * M_PI * t), u.imag(), u_t);
-         sol_sock << "solution\n" << *mesh << u_t
-                  << "window_title '" << oss.str() << "'" << flush;
-         i++;
-      }
-   }
-
-   // 14. Free the used memory.
-   delete a;
-   delete u_exact;
-   delete fespace;
-   delete fec;
-   delete mesh;
-
-   return 0;
-}
-
-bool check_for_inline_mesh(const char * mesh_file)
-{
-   string file(mesh_file);
-   size_t p0 = file.find_last_of("/");
-   string s0 = file.substr((p0==string::npos)?0:(p0+1),7);
-   return s0 == "inline-";
-}
-
-complex<double> u0_exact(const Vector &x)
-{
-   int dim = x.Size();
-   complex<double> i(0.0, 1.0);
-   complex<double> alpha = (epsilon_ * omega_ - i * sigma_);
-   complex<double> kappa = std::sqrt(mu_ * omega_* alpha);
-   return std::exp(-i * kappa * x[dim - 1]);
-}
-
-double u0_real_exact(const Vector &x)
-{
-   return u0_exact(x).real();
-}
-
-double u0_imag_exact(const Vector &x)
-{
-   return u0_exact(x).imag();
-}
-
-void u1_real_exact(const Vector &x, Vector &v)
-{
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[0] = u0_real_exact(x);
-}
-
-void u1_imag_exact(const Vector &x, Vector &v)
-{
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[0] = u0_imag_exact(x);
-}
-
-void u2_real_exact(const Vector &x, Vector &v)
-{
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[dim-1] = u0_real_exact(x);
-}
-
-void u2_imag_exact(const Vector &x, Vector &v)
-{
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[dim-1] = u0_imag_exact(x);
-}

examples/ex21p.cpp

@@ -1,658 +0,0 @@
-//                       MFEM Example 21 - Parallel Version
-//
-// Compile with: make ex21p
-//
-// Sample runs:  mpirun -np 4 ex21p -m ../data/inline-segment.mesh -o 3
-//               mpirun -np 4 ex21p -m ../data/inline-tri.mesh -o 3
-//               mpirun -np 4 ex21p -m ../data/inline-quad.mesh -o 3
-//               mpirun -np 4 ex21p -m ../data/inline-quad.mesh -o 3 -p 1
-//               mpirun -np 4 ex21p -m ../data/inline-quad.mesh -o 3 -p 2
-//               mpirun -np 4 ex21p -m ../data/inline-tet.mesh -o 2
-//               mpirun -np 4 ex21p -m ../data/inline-hex.mesh -o 2
-//               mpirun -np 4 ex21p -m ../data/inline-hex.mesh -o 2 -p 1
-//               mpirun -np 4 ex21p -m ../data/inline-hex.mesh -o 2 -p 2
-//               mpirun -np 4 ex21p -m ../data/star.mesh -o 2 -sigma 10.0
-//
-// Description:  This example code demonstrates the use of MFEM to define and
-//               solve simple complex-valued linear systems.  We implement three
-//               variants of a damped harmonic oscillator:
-//
-//               1) A scalar H1 field
-//                  -Div(a Grad u) - omega^2 b u + i omega c u = 0
-//
-//               2) A vector H(Curl) field
-//                  Curl(a Curl u) - omega^2 b u + i omega c u = 0
-//
-//               3) A vector H(Div) field
-//                  -Grad(a Div u) - omega^2 b u + i omega c u = 0
-//
-//               In each case the field is driven by a forced oscillation, with
-//               angular frequency omega, imposed at the boundary or a portion
-//               of the boundary.
-//
-//               In electromagnetics the coefficients are typically named the
-//               permeability, mu = 1/a, permittivity, epsilon = b, and
-//               conductivity, sigma = c.  The user can specify these constants
-//               using either set of names.
-//
-//#define MFEM_STRUMPACK_SRC
-
-#include <fstream>
-#include <iostream>
-#include "mfem.hpp"
-
-using namespace std;
-using namespace mfem;
-
-static double mu_ = 1.0;
-static double epsilon_ = 1.0;
-static double sigma_ = 20.0;
-static double omega_ = 10.0;
-
-double u0_real_exact(const Vector &);
-double u0_imag_exact(const Vector &);
-
-void u1_real_exact(const Vector &, Vector &);
-void u1_imag_exact(const Vector &, Vector &);
-
-void u2_real_exact(const Vector &, Vector &);
-void u2_imag_exact(const Vector &, Vector &);
-
-bool check_for_inline_mesh(const char * mesh_file);
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI.
-   int num_procs, myid;
-   MPI_Init(&argc, &argv);
-   MPI_Comm comm = MPI_COMM_WORLD;
-   MPI_Comm_size(comm, &num_procs);
-   MPI_Comm_rank(comm, &myid);
-
-   // 2. Parse command-line options.
-   const char *mesh_file = "../data/inline-quad.mesh";
-   int ser_ref_levels = 1;
-   int par_ref_levels = 1;
-   int order = 1;
-   int prob = 0;
-   double freq = -1.0;
-   double a_coef = 0.0;
-   bool visualization = 1;
-   bool herm_conv = true;
-   bool exact_sol = true;
-#ifdef MFEM_USE_STRUMPACK
-   bool strumpack = false;
-#endif
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
-                  "Number of times to refine the mesh uniformly in serial.");
-   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
-                  "Number of times to refine the mesh uniformly in parallel.");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree).");
-   args.AddOption(&prob, "-p", "--problem-type",
-                  "Choose from 0: H_1, 1: H(Curl), or 2: H(Div) "
-                  "damped harmonic oscillator.");
-   args.AddOption(&a_coef, "-a", "--stiffness-coef",
-                  "Stiffness coefficient (spring constant or 1/mu).");
-   args.AddOption(&epsilon_, "-b", "--mass-coef",
-                  "Mass coefficient (or epsilon).");
-   args.AddOption(&sigma_, "-c", "--damping-coef",
-                  "Damping coefficient (or sigma).");
-   args.AddOption(&mu_, "-mu", "--permeability",
-                  "Permeability of free space (or 1/(spring constant)).");
-   args.AddOption(&epsilon_, "-eps", "--permittivity",
-                  "Permittivity of free space (or mass constant).");
-   args.AddOption(&sigma_, "-sigma", "--conductivity",
-                  "Conductivity (or damping constant).");
-   args.AddOption(&freq, "-f", "--frequency",
-                  "Frequency (in Hz).");
-#ifdef MFEM_USE_STRUMPACK
-   args.AddOption(&strumpack, "-strumpack", "--strumpack-solver",
-                  "-no-strumpack", "--no-strumpack-solver",
-                  "Use STRUMPACK's double complex linear solver.");
-#endif
-   args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
-                  "--no-hermitian", "Use convention for Hermitian operators.");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      MPI_Finalize();
-      return 1;
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
-
-   if ( a_coef != 0.0 )
-   {
-      mu_ = 1.0 / a_coef;
-   }
-   if ( freq > 0.0 )
-   {
-      omega_ = 2.0 * M_PI * freq;
-   }
-
-   exact_sol = check_for_inline_mesh(mesh_file);
-   if (myid == 0 && exact_sol)
-   {
-      cout << "Identified an 'inline' mesh" << endl;
-   }
-
-   ComplexOperator::Convention conv =
-      herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
-
-   // 3. Read the (serial) mesh from the given mesh file on all processors.  We
-   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
-   //    and volume meshes with the same code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   int dim = mesh->Dimension();
-
-   // 4. Refine the serial mesh on all processors to increase the resolution.
-   for (int l = 0; l < ser_ref_levels; l++)
-   {
-      mesh->UniformRefinement();
-   }
-
-   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh further in parallel to increase the resolution. Once the
-   //    parallel mesh is defined, the serial mesh can be deleted.
-   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
-   delete mesh;
-   for (int l = 0; l < par_ref_levels; l++)
-   {
-      pmesh->UniformRefinement();
-   }
-
-   // 6. Define a parallel finite element space on the parallel
-   //    mesh. Here we use continuous Lagrange, Nedelec, or
-   //    Raviart-Thomas finite elements of the specified order.
-   if (dim == 1 && prob != 0 )
-   {
-      if (myid == 0)
-      {
-         cout << "Switching to problem type 0, H1 basis functions, "
-              << "for 1 dimensional mesh." << endl;
-      }
-      prob = 0;
-   }
-
-   FiniteElementCollection *fec;
-   switch (prob)
-   {
-      case 0:  fec = new H1_FECollection(order, dim);      break;
-      case 1:  fec = new ND_FECollection(order, dim);      break;
-      case 2:  fec = new RT_FECollection(order - 1, dim);  break;
-   }
-   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
-   HYPRE_Int size = fespace->GlobalTrueVSize();
-   if (myid == 0)
-   {
-      cout << "Number of finite element unknowns: " << size << endl;
-   }
-
-   // 7. Determine the list of true (i.e. parallel conforming) essential
-   //    boundary dofs. In this example, the boundary conditions are defined
-   //    based on the type of mesh and the problem type.
-   Array<int> ess_tdof_list;
-   Array<int> ess_bdr;
-   if (pmesh->bdr_attributes.Size())
-   {
-      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
-      ess_bdr = 1;
-      if (exact_sol)
-      {
-         switch (prob)
-         {
-            case 0:   ess_bdr = 0; ess_bdr[0] = 1;  break;
-            default:  ess_bdr = 1; ess_bdr[2] = 0;  break;
-         }
-      }
-      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-   }
-
-   // 8. Set up the parallel linear form b(.) which corresponds to the
-   //    right-hand side of the FEM linear system.
-   ParComplexLinearForm b(fespace, conv);
-   b.Vector::operator=(0.0);
-
-   // 9. Define the solution vector u as a parallel finite element
-   //    grid function corresponding to fespace. Initialize u with
-   //    initial guess of 1+0i or the exact solution if it is known.
-   ParComplexGridFunction u(fespace);
-   ParComplexGridFunction * u_exact = NULL;
-   if (exact_sol) { u_exact = new ParComplexGridFunction(fespace); }
-
-   FunctionCoefficient u0_r(u0_real_exact);
-   FunctionCoefficient u0_i(u0_imag_exact);
-   VectorFunctionCoefficient u1_r(dim, u1_real_exact);
-   VectorFunctionCoefficient u1_i(dim, u1_imag_exact);
-   VectorFunctionCoefficient u2_r(dim, u2_real_exact);
-   VectorFunctionCoefficient u2_i(dim, u2_imag_exact);
-
-   ConstantCoefficient zeroCoef(0.0);
-   ConstantCoefficient oneCoef(1.0);
-
-   Vector zeroVec(dim); zeroVec = 0.0;
-   Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
-   VectorConstantCoefficient zeroVecCoef(zeroVec);
-   VectorConstantCoefficient oneVecCoef(oneVec);
-
-   switch (prob)
-   {
-      case 0:
-         u.ProjectBdrCoefficient(oneCoef, zeroCoef, ess_bdr);
-         if (exact_sol) { u_exact->ProjectCoefficient(u0_r, u0_i); }
-         break;
-      case 1:
-         u.ProjectBdrCoefficientTangent(oneVecCoef, zeroVecCoef, ess_bdr);
-         if (exact_sol) { u_exact->ProjectCoefficient(u1_r, u1_i); }
-         break;
-      case 2:
-         u.ProjectBdrCoefficientNormal(oneVecCoef, zeroVecCoef, ess_bdr);
-         if (exact_sol) { u_exact->ProjectCoefficient(u2_r, u2_i); }
-         break;
-   }
-
-   if (visualization && exact_sol)
-   {
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      socketstream sol_sock_r(vishost, visport);
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_r.precision(8);
-      sol_sock_i.precision(8);
-      sol_sock_r << "solution\n" << *pmesh << u_exact->real()
-                 << "window_title 'Exact Real Part'" << flush;
-      sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
-                 << "window_title 'Exact Imaginary Part'" << flush;
-   }
-
-   // 10. Set up the parallel sesquilinear form a(.,.) on the finite element
-   //     space corresponding to the damped harmonic oscillator operator
-   //     of the appropriate type:
-   //
-   //     0) A scalar H1 field
-   //        -Div(a Grad) - omega^2 b + i omega c
-   //
-   //     1) A vector H(Curl) field
-   //        Curl(a Curl) - omega^2 b + i omega c
-   //
-   //     2) A vector H(Div) field
-   //        -Grad(a Div) - omega^2 b + i omega c
-   //
-   ConstantCoefficient stiffnessCoef(1.0/mu_);
-   ConstantCoefficient massCoef(-omega_ * omega_ * epsilon_);
-   ConstantCoefficient lossCoef(omega_ * sigma_);
-   ConstantCoefficient negMassCoef(omega_ * omega_ * epsilon_);
-
-   ParSesquilinearForm *a = new ParSesquilinearForm(fespace, conv);
-   switch (prob)
-   {
-      case 0:
-         a->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new MassIntegrator(massCoef),
-                                new MassIntegrator(lossCoef));
-         break;
-      case 1:
-         a->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
-                                new VectorFEMassIntegrator(lossCoef));
-         break;
-      case 2:
-         a->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef),
-                                NULL);
-         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
-                                new VectorFEMassIntegrator(lossCoef));
-         break;
-   }
-
-   // 10a. Set up the parallel bilinear form for the preconditioner
-   //      corresponding to the appropriate operator if the STRUMPACK solver
-   //      has not been selected.
-   //
-   //      0) A scalar H1 field
-   //         -Div(a Grad) - omega^2 b + omega c
-   //
-   //      1) A vector H(Curl) field
-   //         Curl(a Curl) + omega^2 b + omega c
-   //
-   //      2) A vector H(Div) field
-   //         -Grad(a Div) - omega^2 b + omega c
-   //
-   ParBilinearForm *pcOp = NULL;
-#ifdef MFEM_USE_STRUMPACK
-   if (!strumpack)
-#endif
-   {
-      pcOp = new ParBilinearForm(fespace);
-      switch (prob)
-      {
-         case 0:
-            pcOp->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef));
-            pcOp->AddDomainIntegrator(new MassIntegrator(massCoef));
-            pcOp->AddDomainIntegrator(new MassIntegrator(lossCoef));
-            break;
-         case 1:
-            pcOp->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef));
-            pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(negMassCoef));
-            pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(lossCoef));
-            break;
-         case 2:
-            pcOp->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef));
-            pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef));
-            pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(lossCoef));
-            break;
-      }
-   }
-
-   // 11. Assemble the parallel bilinear form and the corresponding linear
-   //     system, applying any necessary transformations such as: parallel
-   //     assembly, eliminating boundary conditions, applying conforming
-   //     constraints for non-conforming AMR, etc.
-   a->Assemble();
-   if (pcOp) { pcOp->Assemble(); }
-
-   OperatorHandle A;
-   Vector B, U;
-
-   a->FormLinearSystem(ess_tdof_list, u, b, A, U, B);
-   u = 0.0;
-   U = 0.0;
-
-   OperatorHandle PCOp;
-   if (pcOp) { pcOp->FormSystemMatrix(ess_tdof_list, PCOp); }
-
-   if (myid == 0)
-   {
-      ComplexHypreParMatrix * Ahyp =
-         dynamic_cast<ComplexHypreParMatrix*>(A.Ptr());
-
-      cout << "Size of linear system: "
-           << 2 * Ahyp->real().GetGlobalNumRows() << endl << endl;
-   }
-
-   // 12. Define and apply a parallel FGMRES solver for AU=B with a
-   //     block diagonal preconditioner based on the appropriate multigrid
-   //     preconditioner from hypre or simply use STRUMPACK.
-#ifdef MFEM_USE_STRUMPACK
-   if (!strumpack)
-#endif
-   {
-      Array<HYPRE_Int> blockTrueOffsets;
-      blockTrueOffsets.SetSize(3);
-      blockTrueOffsets[0] = 0;
-      blockTrueOffsets[1] = PCOp.Ptr()->Height();
-      blockTrueOffsets[2] = PCOp.Ptr()->Height();
-      blockTrueOffsets.PartialSum();
-
-      BlockDiagonalPreconditioner BDP(blockTrueOffsets);
-
-      Operator * pc_r = NULL;
-      Operator * pc_i = NULL;
-
-      switch (prob)
-      {
-         case 0:
-            pc_r =
-               new HypreBoomerAMG(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()));
-            pc_i = new ScaledOperator(pc_r,
-                                      (conv == ComplexOperator::HERMITIAN) ?
-                                      1.0:-1.0);
-            break;
-         case 1:
-            pc_r = new HypreAMS(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()),
-                                fespace);
-            pc_i = new ScaledOperator(pc_r,
-                                      (conv == ComplexOperator::HERMITIAN) ?
-                                      1.0:-1.0);
-            break;
-         case 2:
-            if (dim == 2 )
-            {
-               pc_r = new HypreAMS(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()),
-                                   fespace);
-            }
-            else
-            {
-               pc_r = new HypreADS(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()),
-                                   fespace);
-            }
-            pc_i = new ScaledOperator(pc_r,
-                                      (conv == ComplexOperator::HERMITIAN) ?
-                                      1.0:-1.0);
-            break;
-      }
-      BDP.SetDiagonalBlock(0, pc_r);
-      BDP.SetDiagonalBlock(1, pc_i);
-      BDP.owns_blocks = 0;
-
-      FGMRESSolver fgmres(MPI_COMM_WORLD);
-      fgmres.SetPreconditioner(BDP);
-      fgmres.SetOperator(*A.Ptr());
-      fgmres.SetRelTol(1e-12);
-      fgmres.SetMaxIter(1000);
-      fgmres.SetPrintLevel(1);
-      fgmres.Mult(B, U);
-   }
-#ifdef MFEM_USE_STRUMPACK
-   else
-   {
-      ComplexHypreParMatrix * Ahyp =
-         dynamic_cast<ComplexHypreParMatrix*>(A.Ptr());
-
-      STRUMPACKRowLocCmplxMatrix A_strmp(Ahyp->real(), Ahyp->imag());
-
-      STRUMPACKCmplxSolver strmp(argc, argv, comm);
-
-      strmp.SetPrintFactorStatistics(true);
-      strmp.SetPrintSolveStatistics(true);
-      // strmp.SetKrylovSolver(strumpack::KrylovSolver::AUTO); // core dump
-      strmp.SetKrylovSolver(strumpack::KrylovSolver::DIRECT); // core dump
-      // strmp.SetKrylovSolver(strumpack::KrylovSolver::REFINE); // core dump
-      // strmp.SetKrylovSolver(strumpack::KrylovSolver::PREC_GMRES); // index out of range asserts from strumpack::DenseMatrix
-      // strmp.SetKrylovSolver(strumpack::KrylovSolver::GMRES); // WORKS
-      strmp.SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
-      strmp.SetOperator(A_strmp);
-      strmp.SetFromCommandLine();
-      strmp.Mult(B, U);
-   }
-#endif
-
-   // 13. Recover the parallel grid function corresponding to U. This is the
-   //     local finite element solution on each processor.
-   a->RecoverFEMSolution(U, b, u);
-
-   if (exact_sol)
-   {
-      double err_r = -1.0;
-      double err_i = -1.0;
-
-      switch (prob)
-      {
-         case 0:
-            err_r = u.real().ComputeL2Error(u0_r);
-            err_i = u.imag().ComputeL2Error(u0_i);
-            break;
-         case 1:
-            err_r = u.real().ComputeL2Error(u1_r);
-            err_i = u.imag().ComputeL2Error(u1_i);
-            break;
-         case 2:
-            err_r = u.real().ComputeL2Error(u2_r);
-            err_i = u.imag().ComputeL2Error(u2_i);
-            break;
-      }
-
-      if ( myid == 0 )
-      {
-         cout << endl;
-         cout << "|| Re (u_h - u) ||_{L^2} = " << err_r << endl;
-         cout << "|| Im (u_h - u) ||_{L^2} = " << err_i << endl;
-         cout << endl;
-      }
-   }
-
-   // 14. Save the refined mesh and the solution in parallel. This output can be
-   //     viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
-   {
-      ostringstream mesh_name, sol_r_name, sol_i_name;
-      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
-      sol_r_name << "sol_r." << setfill('0') << setw(6) << myid;
-      sol_i_name << "sol_i." << setfill('0') << setw(6) << myid;
-
-      ofstream mesh_ofs(mesh_name.str().c_str());
-      mesh_ofs.precision(8);
-      pmesh->Print(mesh_ofs);
-
-      ofstream sol_r_ofs(sol_r_name.str().c_str());
-      ofstream sol_i_ofs(sol_i_name.str().c_str());
-      sol_r_ofs.precision(8);
-      sol_i_ofs.precision(8);
-      u.real().Save(sol_r_ofs);
-      u.imag().Save(sol_i_ofs);
-   }
-
-   // 15. Send the solution by socket to a GLVis server.
-   if (visualization)
-   {
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      socketstream sol_sock_r(vishost, visport);
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_r.precision(8);
-      sol_sock_i.precision(8);
-      sol_sock_r << "solution\n" << *pmesh << u.real()
-                 << "window_title 'Comp Real Part'" << flush;
-      sol_sock_i << "solution\n" << *pmesh << u.imag()
-                 << "window_title 'Comp Imaginary Part'" << flush;
-   }
-   if (visualization && exact_sol)
-   {
-      *u_exact -= u;
-
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      socketstream sol_sock_r(vishost, visport);
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_r.precision(8);
-      sol_sock_i.precision(8);
-      sol_sock_r << "solution\n" << *pmesh << u_exact->real()
-                 << "window_title 'Exact-Comp Real Part'" << flush;
-      sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
-                 << "window_title 'Exact-Comp Imaginary Part'" << flush;
-   }
-   if (visualization)
-   {
-      ParGridFunction u_t(fespace);
-      u_t = u.real();
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      socketstream sol_sock(vishost, visport);
-      sol_sock << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock.precision(8);
-      sol_sock << "solution\n" << *pmesh << u_t
-               << "window_title 'Harmonic Solution (t = 0.0 T)'"
-               << "pause\n" << flush;
-      if (myid == 0)
-         cout << "GLVis visualization paused."
-              << " Press space (in the GLVis window) to resume it.\n";
-      int num_frames = 32;
-      int i = 0;
-      while (sol_sock)
-      {
-         double t = (double)(i % num_frames) / num_frames;
-         ostringstream oss;
-         oss << "Harmonic Solution (t = " << t << " T)";
-
-         add(cos( 2.0 * M_PI * t), u.real(),
-             sin(-2.0 * M_PI * t), u.imag(), u_t);
-         sol_sock << "parallel " << num_procs << " " << myid << "\n";
-         sol_sock << "solution\n" << *pmesh << u_t
-                  << "window_title '" << oss.str() << "'" << flush;
-         i++;
-      }
-   }
-
-   // 16. Free the used memory.
-   delete a;
-   delete u_exact;
-   delete pcOp;
-   delete fespace;
-   delete fec;
-   delete pmesh;
-
-   MPI_Finalize();
-
-   return 0;
-}
-
-bool check_for_inline_mesh(const char * mesh_file)
-{
-   string file(mesh_file);
-   size_t p0 = file.find_last_of("/");
-   string s0 = file.substr((p0==string::npos)?0:(p0+1),7);
-   return s0 == "inline-";
-}
-
-complex<double> u0_exact(const Vector &x)
-{
-   int dim = x.Size();
-   complex<double> i(0.0, 1.0);
-   complex<double> alpha = (epsilon_ * omega_ - i * sigma_);
-   complex<double> kappa = std::sqrt(mu_ * omega_* alpha);
-   return std::exp(-i * kappa * x[dim - 1]);
-}
-
-double u0_real_exact(const Vector &x)
-{
-   return u0_exact(x).real();
-}
-
-double u0_imag_exact(const Vector &x)
-{
-   return u0_exact(x).imag();
-}
-
-void u1_real_exact(const Vector &x, Vector &v)
-{
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[0] = u0_real_exact(x);
-}
-
-void u1_imag_exact(const Vector &x, Vector &v)
-{
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[0] = u0_imag_exact(x);
-}
-
-void u2_real_exact(const Vector &x, Vector &v)
-{
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[dim-1] = u0_real_exact(x);
-}
-
-void u2_imag_exact(const Vector &x, Vector &v)
-{
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[dim-1] = u0_imag_exact(x);
-}

examples/ex23.cpp

@@ -0,0 +1,734 @@
+//                               MFEM Example 23
+//
+// Compile with: make ex23
+//
+// Sample runs:
+//    ex23 -m ../data/periodic-segment.mesh -p 0 -s 2 -dt 0.001 -vs 50
+//    ex23 -m ../data/periodic-segment.mesh -p 0 -s 12 -dt 0.01
+//    ex23 -m ../data/periodic-segment.mesh -p 0 -s 22 -dt 0.01
+//    ex23 -m ../data/periodic-segment.mesh -p 0 -s 32 -dt 0.005 -vs 10
+//    ex23 -m ../data/periodic-square.mesh -p 0 -dt 0.01
+//    ex23 -m ../data/periodic-square.mesh -p 0 -s 32 -dt 0.01
+//    ex23 -m ../data/periodic-hexagon.mesh -p 0 -d 0.001 -s 12 -dt 0.02
+//    ex23 -m ../data/periodic-hexagon.mesh -p 0 -d 0.001 -s 32 -dt 0.009 -vs 10
+//    ex23 -m ../data/periodic-square.mesh -p 1 -dt 0.01 -tf 9
+//    ex23 -m ../data/periodic-hexagon.mesh -p 1 -dt 0.01 -tf 9
+//    ex23 -m ../data/amr-quad.mesh -p 1 -dt 0.01 -tf 9 -vs 2
+//    ex23 -m ../data/disc-nurbs.mesh -p 1 -r 3 -dt 0.01 -tf 9
+//    ex23 -m ../data/disc-nurbs.mesh -p 2 -r 3 -dt 0.01 -tf 9
+//    ex23 -m ../data/disc-nurbs.mesh -p 3 -r 3 -dt 0.01 -tf 9 -d 0.02
+//    ex23 -m ../data/periodic-square.mesh -p 3 -r 3 -dt 0.025 -tf 9
+//    ex23 -m ../data/periodic-cube.mesh -p 0 -o 2 -dt 0.025 -tf 8
+//
+// Description:  This example code solves the time-dependent advection-diffusion
+//               equation
+//               du/dt - div(D grad(u)) + v.grad(u) = 0, where
+//               D is a diffusion coefficient,
+//               v is a given fluid velocity, and
+//               u0(x)=u(0,x) is a given initial condition.
+//
+//               The example demonstrates the use of Discontinuous Galerkin (DG)
+//               bilinear forms in MFEM (face integrators), the use of explicit,
+//               implicit, and implicit-explicit ODE time integrators, the
+//               definition of periodic boundary conditions through periodic
+//               meshes, as well as the use of GLVis for persistent
+//               visualization of a time-evolving solution. The saving of
+//               time-dependent data files for external visualization with
+//               VisIt (visit.llnl.gov) is also illustrated.
+//
+//               This example is a merger of examples 9 and 14.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Choice for the problem setup. The fluid velocity, initial condition and
+// boundary condition are chosen based on this parameter.
+int problem;
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v);
+
+// Initial condition
+double u0_function(const Vector &x);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+
+/** A time-dependent operator for the right-hand side of the ODE for use with
+    explicit ODE solvers. The DG weak form of du/dt = div(D grad(u))-v.grad(u) is
+    M du/dt = - S u + K u + b, where M, S, and K are the mass,
+    stiffness, and advection matrices, and b describes sources and the flow on
+    the boundary.
+    This can be written as a general ODE,
+    du/dt = M^{-1} (-S u + K u + b), and this class is used to compute the RHS
+    and perform the solve for du/dt. */
+class EX_Evolution : public TimeDependentOperator
+{
+private:
+   SparseMatrix &M, &S, &K;
+   const Vector &b;
+
+   DSmoother M_prec;
+   CGSolver M_solver;
+
+   mutable Vector z;
+
+   void initA(double dt);
+
+public:
+   EX_Evolution(SparseMatrix &_M, SparseMatrix &_S, SparseMatrix &_K,
+                const Vector &_b);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual ~EX_Evolution() {}
+};
+
+/** A time-dependent operator for the right-hand side of the ODE for use with
+    implicit ODE solvers. The DG weak form of du/dt = div(D grad(u))-v.grad(u) is
+    [M + dt (S - K)] du/dt = - S u + K u + b, where M, S, and K are the mass,
+    stiffness, and advection matrices, and b describes sources and the flow on
+    the boundary.
+    This can be written as a general ODE,
+    du/dt = A^{-1} (-S u + K u + b) with A = [M + dt (S - K)], and this class is
+    used to perform the fully implicit solve for du/dt. */
+class IM_Evolution : public TimeDependentOperator
+{
+private:
+   SparseMatrix &M, &S, &K;
+   SparseMatrix *A;
+   const Vector &b;
+
+   DSmoother M_prec;
+   CGSolver M_solver;
+
+   DSmoother *A_prec;
+   GMRESSolver *A_solver;
+   double dt;
+
+   mutable Vector z;
+
+   void initA(double dt);
+
+public:
+   IM_Evolution(SparseMatrix &_M, SparseMatrix &_S, SparseMatrix &_K,
+                const Vector &_b);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &y);
+
+   virtual ~IM_Evolution() { delete A_solver; delete A_prec; delete A; }
+};
+
+/** A time-dependent operator for the right-hand side of the ODE for use with
+    IMEX (Implicit-Explicit) ODE solvers. The DG weak form of
+    du/dt = div(D grad(u))-v.grad(u) is
+    [M + dt S] du/dt = - S u + K u + b, where M, S, and K are the mass,
+    stiffness, and advection matrices, and b describes sources and the flow on
+    the boundary.
+    This can be written as a general ODE,
+    du/dt = A^{-1} (-S u + K u + b) with A = [M + dt (S - K)], and this class is
+    used to perform the implicit or explicit solve for du/dt. */
+class IMEX_Evolution : public TimeDependentOperator
+{
+private:
+   SparseMatrix &M, &S, &K;
+   SparseMatrix *A;
+   const Vector &b;
+
+   DSmoother M_prec;
+   CGSolver M_solver;
+
+   DSmoother *A_prec;
+   CGSolver *A_solver;
+   double dt;
+
+   mutable Vector z;
+
+   void initA(double dt);
+
+public:
+   IMEX_Evolution(SparseMatrix &_M, SparseMatrix &_S, SparseMatrix &_K,
+                  const Vector &_b);
+
+   virtual void ExplicitMult(const Vector &x, Vector &y) const;
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &y);
+
+   virtual ~IMEX_Evolution() { delete A_solver; delete A_prec; delete A; }
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   problem = 0;
+   const char *mesh_file = "../data/periodic-hexagon.mesh";
+   int ref_levels = 2;
+   int order = 3;
+   int ode_solver_type = 12;
+   double t_final = 10.0;
+   double d_coef = 0.01;
+   double dt = 0.01;
+   double sigma = -1.0;
+   double kappa = -1.0;
+   bool visualization = true;
+   bool visit = false;
+   bool binary = false;
+   int vis_steps = 5;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver: 1 - Forward Euler, 2 - RK2, 3 - RK3 SSP,"
+                  " 4 - RK4, 5 - Generalized Alpha,\n\t"
+                  "11 - Backward Euler, 12 - SDIRK2, 13 - SDIRK3,\n\t"
+                  "22 - Implicit Midpoint, 23 SDIRK23, 24 - SDIRK34,\n\t"
+                  "31 - IMEX BE/FE, 32 - IMEX RK2.");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&d_coef, "-d", "--diff-coef",
+                  "Diffusion coefficient.");
+   args.AddOption(&sigma, "-s", "--sigma",
+                  "One of the two DG penalty parameters, typically +1/-1."
+                  " See the documentation of class DGDiffusionIntegrator.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "One of the two DG penalty parameters, should be positive."
+                  " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
+                  "--no-visit-datafiles",
+                  "Save data files for VisIt (visit.llnl.gov) visualization.");
+   args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
+                  "--ascii-datafiles",
+                  "Use binary (Sidre) or ascii format for VisIt data files.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   if (kappa < 0)
+   {
+      kappa = (order+1)*(order+1);
+   }
+   args.PrintOptions(cout);
+
+   // 2. Define the ODE solver used for time integration. Several explicit
+   //    Runge-Kutta methods are available.
+   ODESolver *ode_solver = NULL;
+   switch (ode_solver_type)
+   {
+      // Explicit methods
+      case 1:  ode_solver = new ForwardEulerSolver; break;
+      case 2:  ode_solver = new RK2Solver(0.5); break; // midpoint method
+      case 3:  ode_solver = new RK3SSPSolver; break;
+      case 4:  ode_solver = new RK4Solver; break;
+      case 5:  ode_solver = new GeneralizedAlphaSolver(0.5); break;
+      // Implicit L-stable methods
+      case 11: ode_solver = new BackwardEulerSolver; break;
+      case 12: ode_solver = new SDIRK23Solver(2); break;
+      case 13: ode_solver = new SDIRK33Solver; break;
+      // Implicit A-stable methods (not L-stable)
+      case 22: ode_solver = new ImplicitMidpointSolver; break;
+      case 23: ode_solver = new SDIRK23Solver; break;
+      case 24: ode_solver = new SDIRK34Solver; break;
+      // Implicit-Explicit methods
+      case 31: ode_solver = new IMEX_BE_FE; break;
+      case 32: ode_solver = new IMEXRK2; break;
+      default:
+         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         return 3;
+   }
+
+   // 3. Read the serial mesh from the given mesh file on all processors. We can
+   //    handle geometrically periodic meshes in this code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 4. Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter. If the mesh is of NURBS type, we convert it
+   //    to a (piecewise-polynomial) high-order mesh.
+   for (int lev = 0; lev < ref_levels; lev++)
+   {
+      mesh.UniformRefinement();
+   }
+   if (mesh.NURBSext)
+   {
+      mesh.SetCurvature(max(order, 1));
+   }
+   mesh.GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 5. Define the parallel discontinuous DG finite element space on the
+   //    parallel refined mesh of the given polynomial order.
+   DG_FECollection fec(order, dim);
+   FiniteElementSpace fes(&mesh, &fec);
+
+   cout << "Number of unknowns: " << fes.GetVSize() << endl;
+
+   // 6. Set up and assemble the parallel bilinear and linear forms (and the
+   //    parallel hypre matrices) corresponding to the DG discretization. The
+   //    DGTraceIntegrator involves integrals over mesh interior faces.
+   ConstantCoefficient diff_coef(d_coef);
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient u0(u0_function);
+
+   BilinearForm m(&fes);
+   m.AddDomainIntegrator(new MassIntegrator);
+
+   BilinearForm s(&fes);
+   s.AddDomainIntegrator(new DiffusionIntegrator(diff_coef));
+   s.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma,
+                                                         kappa));
+   s.AddBdrFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma, kappa));
+
+   BilinearForm k(&fes);
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   k.AddInteriorFaceIntegrator(
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+   k.AddBdrFaceIntegrator(
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+
+   LinearForm b(&fes);
+   b.AddBdrFaceIntegrator(
+      new DGDirichletLFIntegrator(u0, diff_coef, sigma, kappa));
+
+   int skip_zeros = 0;
+   m.Assemble(skip_zeros);
+   m.Finalize(skip_zeros);
+   s.Assemble(skip_zeros);
+   s.Finalize(skip_zeros);
+   k.Assemble(skip_zeros);
+   k.Finalize(skip_zeros);
+   b.Assemble();
+
+   // 7. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   GridFunction u(&fes);
+   u.ProjectCoefficient(u0);
+
+   {
+      ofstream omesh("ex23.mesh");
+      omesh.precision(precision);
+      mesh.Print(omesh);
+      ofstream osol("ex23-init.gf");
+      osol.precision(precision);
+      u.Save(osol);
+   }
+
+   // Create data collection for solution output: either VisItDataCollection for
+   // ascii data files, or SidreDataCollection for binary data files.
+   DataCollection *dc = NULL;
+   if (visit)
+   {
+      if (binary)
+      {
+#ifdef MFEM_USE_SIDRE
+         dc = new SidreDataCollection("Example23", &mesh);
+#else
+         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+#endif
+      }
+      else
+      {
+         dc = new VisItDataCollection("Example23", &mesh);
+         dc->SetPrecision(precision);
+      }
+      dc->RegisterField("solution", &u);
+      dc->SetCycle(0);
+      dc->SetTime(0.0);
+      dc->Save();
+   }
+
+   socketstream sout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      if (!sout)
+      {
+         cout << "Unable to connect to GLVis server at "
+              << vishost << ':' << visport << endl;
+         visualization = false;
+         cout << "GLVis visualization disabled.\n";
+      }
+      else
+      {
+         sout.precision(precision);
+         sout << "solution\n" << mesh << u;
+         sout << "pause\n";
+         sout << flush;
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+   // 8. Define the time-dependent evolution operator describing the ODE
+   //    right-hand side, and perform time-integration (looping over the time
+   //    iterations, ti, with a time-step dt).
+
+   TimeDependentOperator *adv = NULL;
+   if (ode_solver_type < 10)
+   {
+      adv = new EX_Evolution(m.SpMat(), s.SpMat(), k.SpMat(), b);
+   }
+   else if (ode_solver_type < 30)
+   {
+      adv = new IM_Evolution(m.SpMat(), s.SpMat(), k.SpMat(), b);
+   }
+   else
+   {
+      adv = new IMEX_Evolution(m.SpMat(), s.SpMat(), k.SpMat(), b);
+   }
+
+   double t = 0.0;
+   adv->SetTime(t);
+   ode_solver->Init(*adv);
+
+   int n_steps = (int)ceil(t_final / dt);
+   double dt_real = t_final / n_steps;
+
+   for (int ti = 0; ti < n_steps; )
+   {
+      ode_solver->Step(u, t, dt_real);
+      ti++;
+
+      if (ti % vis_steps == 0 || ti == n_steps)
+      {
+         cout << "time step: " << ti << ", time: " << t << endl;
+
+         if (visualization)
+         {
+            sout << "solution\n" << mesh << u << flush;
+         }
+
+         if (visit)
+         {
+            dc->SetCycle(ti);
+            dc->SetTime(t);
+            dc->Save();
+         }
+      }
+   }
+
+   // 9. Save the final solution in parallel. This output can be viewed later
+   //    using GLVis: "glvis -np <np> -m ex23-mesh -g ex23-final".
+   {
+      ofstream osol("ex23-final.gf");
+      osol.precision(precision);
+      u.Save(osol);
+   }
+
+   // 10. Free the used memory.
+   delete ode_solver;
+   delete adv;
+   delete dc;
+
+   return 0;
+}
+
+
+// Implementation of class EX_Evolution
+EX_Evolution::EX_Evolution(SparseMatrix &_M, SparseMatrix &_S,
+                           SparseMatrix &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     M(_M), S(_S), K(_K), b(_b), z(_M.Height())
+{
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void EX_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (-S x + K x + b)
+   K.Mult(x, z);
+   S.AddMult(x, z, -1.0);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+// Implementation of class IM_Evolution
+IM_Evolution::IM_Evolution(SparseMatrix &_M, SparseMatrix &_S,
+                           SparseMatrix &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     M(_M), S(_S), K(_K), A(NULL), b(_b),
+     A_prec(NULL), A_solver(NULL), dt(-1.0), z(_M.Height())
+{
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void IM_Evolution::initA(double _dt)
+{
+   if (fabs(dt - _dt) > 1e-4 * _dt)
+   {
+      delete A_solver;
+      delete A_prec;
+      delete A;
+
+      SparseMatrix * SK = Add(1.0, S, -1.0, K);
+      A = Add(1.0, M, _dt, *SK);
+      delete SK;
+      dt = _dt;
+
+      A_prec = new DSmoother(*A);
+      A_solver = new GMRESSolver;
+      A_solver->SetOperator(*A);
+      A_solver->SetPreconditioner(*A_prec);
+
+      A_solver->iterative_mode = false;
+      A_solver->SetRelTol(1e-9);
+      A_solver->SetAbsTol(0.0);
+      A_solver->SetMaxIter(100);
+      A_solver->SetPrintLevel(0);
+   }
+}
+
+void IM_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (-S x + K x + b)
+   K.Mult(x, z);
+   S.AddMult(x, z, -1.0);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IM_Evolution::ImplicitSolve(const double _dt, const Vector &x, Vector &y)
+{
+   this->initA(_dt);
+
+   // y = (M + dt S - dt K)^{-1} (-S x + K x + b)
+   K.Mult(x, z);
+   S.AddMult(x, z, -1.0);
+   z += b;
+   A_solver->Mult(z, y);
+}
+
+// Implementation of class IMEX_Evolution
+IMEX_Evolution::IMEX_Evolution(SparseMatrix &_M, SparseMatrix &_S,
+                               SparseMatrix &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     M(_M), S(_S), K(_K), A(NULL), b(_b),
+     A_prec(NULL), A_solver(NULL), dt(-1.0), z(_M.Height())
+{
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void IMEX_Evolution::initA(double _dt)
+{
+   if (fabs(dt - _dt) > 1e-4 * _dt)
+   {
+      delete A_solver;
+      delete A_prec;
+      delete A;
+
+      A = Add(_dt, S, 1.0, M); // A = M + dt * S
+      dt = _dt;
+
+      A_prec = new DSmoother(*A);
+      A_solver = new CGSolver;
+      A_solver->SetOperator(*A);
+      A_solver->SetPreconditioner(*A_prec);
+
+      A_solver->iterative_mode = false;
+      A_solver->SetRelTol(1e-9);
+      A_solver->SetAbsTol(0.0);
+      A_solver->SetMaxIter(100);
+      A_solver->SetPrintLevel(0);
+   }
+}
+
+void IMEX_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (-S x + K x + b)
+   K.Mult(x, z);
+   S.AddMult(x, z, -1.0);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IMEX_Evolution::ExplicitMult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (K x + b)
+   K.Mult(x, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IMEX_Evolution::ImplicitSolve(const double _dt, const Vector &x, Vector &y)
+{
+   this->initA(_dt);
+   // y = (M + dt S)^{-1} (-S x + b)
+   S.Mult(x, z);
+   z *= -1.0;
+   z += b;
+   A_solver->Mult(z, y);
+}
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      {
+         // Translations in 1D, 2D, and 3D
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
+            case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
+               break;
+         }
+         break;
+      }
+      case 1:
+      case 2:
+      {
+         // Clockwise rotation in 2D around the origin
+         const double w = M_PI/2;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
+            case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+      case 3:
+      {
+         // Clockwise twisting rotation in 2D around the origin
+         const double w = M_PI/2;
+         double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
+         d = d*d;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
+            case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+   }
+}
+
+// Initial condition
+double u0_function(const Vector &x)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      {
+         switch (dim)
+         {
+            case 1:
+               return exp(-40.*pow(X(0)-0.5,2));
+            case 2:
+            case 3:
+            {
+               double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
+               if (dim == 3)
+               {
+                  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
+                  rx *= s;
+                  ry *= s;
+               }
+               return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
+                        erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
+            }
+         }
+      }
+      case 2:
+      {
+         double x_ = X(0), y_ = X(1), rho, phi;
+         rho = hypot(x_, y_);
+         phi = atan2(y_, x_);
+         return pow(sin(M_PI*rho),2)*sin(3*phi);
+      }
+      case 3:
+      {
+         const double f = M_PI;
+         return sin(f*X(0))*sin(f*X(1));
+      }
+   }
+   return 0.0;
+}
+
+// Inflow boundary condition (zero for the problems considered in this example)
+double inflow_function(const Vector &x)
+{
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      case 2:
+      case 3: return 0.0;
+   }
+   return 0.0;
+}

examples/ex23p.cpp

@@ -0,0 +1,797 @@
+//                       MFEM Example 23 - Parallel Version
+//
+// Compile with: make ex23p
+//
+// Sample runs:
+//    mpirun -np 4 ex23p -m ../data/periodic-segment.mesh -p 0 -s 2 -dt 0.001 -vs 50
+//    mpirun -np 4 ex23p -m ../data/periodic-segment.mesh -p 0 -s 12 -dt 0.01
+//    mpirun -np 4 ex23p -m ../data/periodic-segment.mesh -p 0 -s 22 -dt 0.01
+//    mpirun -np 4 ex23p -m ../data/periodic-segment.mesh -p 0 -s 32 -dt 0.005 -vs 10
+//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 0 -dt 0.01
+//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 0 -s 32 -dt 0.01
+//    mpirun -np 4 ex23p -m ../data/periodic-hexagon.mesh -p 0 -d 0.001 -s 12 -dt 0.02
+//    mpirun -np 4 ex23p -m ../data/periodic-hexagon.mesh -p 0 -d 0.001 -s 32 -dt 0.009 -vs 10
+//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 1 -dt 0.01 -tf 9
+//    mpirun -np 4 ex23p -m ../data/periodic-hexagon.mesh -p 1 -dt 0.01 -tf 9
+//    mpirun -np 4 ex23p -m ../data/amr-quad.mesh -p 1 -dt 0.01 -tf 9 -vs 2
+//    mpirun -np 4 ex23p -m ../data/disc-nurbs.mesh -p 1 -rp 1 -dt 0.01 -tf 9
+//    mpirun -np 4 ex23p -m ../data/disc-nurbs.mesh -p 2 -rp 1 -dt 0.01 -tf 9
+//    mpirun -np 4 ex23p -m ../data/disc-nurbs.mesh -p 3 -rp 1 -dt 0.01 -tf 9 -d 0.02
+//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 3 -rp 1 -dt 0.025 -tf 9
+//    mpirun -np 4 ex23p -m ../data/periodic-cube.mesh -p 0 -o 2 -dt 0.025 -tf 8
+//
+// Description:  This example code solves the time-dependent advection-diffusion
+//               equation
+//               du/dt - div(D grad(u)) + v.grad(u) = 0, where
+//               D is a diffusion coefficient,
+//               v is a given fluid velocity, and
+//               u0(x)=u(0,x) is a given initial condition.
+//
+//               The example demonstrates the use of Discontinuous Galerkin (DG)
+//               bilinear forms in MFEM (face integrators), the use of explicit,
+//               implicit, and implicit-explicit ODE time integrators, the
+//               definition of periodic boundary conditions through periodic
+//               meshes, as well as the use of GLVis for persistent
+//               visualization of a time-evolving solution. The saving of
+//               time-dependent data files for external visualization with
+//               VisIt (visit.llnl.gov) is also illustrated.
+//
+//               This example is a merger of examples 9 and 14.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Choice for the problem setup. The fluid velocity, initial condition and
+// boundary condition are chosen based on this parameter.
+int problem;
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v);
+
+// Initial condition
+double u0_function(const Vector &x);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+/** A time-dependent operator for the right-hand side of the ODE for use with
+    explicit ODE solvers. The DG weak form of du/dt = div(D grad(u))-v.grad(u) is
+    M du/dt = - S u + K u + b, where M, S, and K are the mass,
+    stiffness, and advection matrices, and b describes sources and the flow on
+    the boundary.
+    This can be written as a general ODE,
+    du/dt = M^{-1} (-S u + K u + b), and this class is used to compute the RHS
+    and perform the solve for du/dt. */
+class EX_Evolution : public TimeDependentOperator
+{
+private:
+   HypreParMatrix &M, &S, &K;
+   const Vector &b;
+
+   HypreSmoother M_prec;
+   CGSolver M_solver;
+
+   mutable Vector z;
+
+   void initA(double dt);
+
+public:
+   EX_Evolution(HypreParMatrix &_M, HypreParMatrix &_S, HypreParMatrix &_K,
+                const Vector &_b);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual ~EX_Evolution() {}
+};
+
+/** A time-dependent operator for the right-hand side of the ODE for use with
+    implicit ODE solvers. The DG weak form of du/dt = div(D grad(u))-v.grad(u) is
+    [M + dt (S - K)] du/dt = - S u + K u + b, where M, S, and K are the mass,
+    stiffness, and advection matrices, and b describes sources and the flow on
+    the boundary.
+    This can be written as a general ODE,
+    du/dt = A^{-1} (-S u + K u + b) with A = [M + dt (S - K)], and this class is
+    used to perform the fully implicit solve for du/dt. */
+class IM_Evolution : public TimeDependentOperator
+{
+private:
+   HypreParMatrix &M, &S, &K;
+   HypreParMatrix *A;
+   const Vector &b;
+
+   HypreSmoother M_prec;
+   CGSolver M_solver;
+
+   HypreBoomerAMG *A_prec;
+   GMRESSolver *A_solver;
+   double dt;
+
+   mutable Vector z;
+
+   void initA(double dt);
+
+public:
+   IM_Evolution(HypreParMatrix &_M, HypreParMatrix &_S, HypreParMatrix &_K,
+                const Vector &_b);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &y);
+
+   virtual ~IM_Evolution() { delete A_solver; delete A_prec; delete A; }
+};
+
+/** A time-dependent operator for the right-hand side of the ODE for use with
+    IMEX (Implicit-Explicit) ODE solvers. The DG weak form of
+    du/dt = div(D grad(u))-v.grad(u) is
+    [M + dt S] du/dt = - S u + K u + b, where M, S, and K are the mass,
+    stiffness, and advection matrices, and b describes sources and the flow on
+    the boundary.
+    This can be written as a general ODE,
+    du/dt = A^{-1} (-S u + K u + b) with A = [M + dt (S - K)], and this class is
+    used to perform the implicit or explicit solve for du/dt. */
+class IMEX_Evolution : public TimeDependentOperator
+{
+private:
+   HypreParMatrix &M, &S, &K;
+   HypreParMatrix *A;
+   const Vector &b;
+
+   HypreSmoother M_prec;
+   CGSolver M_solver;
+
+   HypreBoomerAMG *A_prec;
+   CGSolver *A_solver;
+   double dt;
+
+   mutable Vector z;
+
+   void initA(double dt);
+
+public:
+   IMEX_Evolution(HypreParMatrix &_M, HypreParMatrix &_S, HypreParMatrix &_K,
+                  const Vector &_b);
+
+   virtual void ExplicitMult(const Vector &x, Vector &y) const;
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &y);
+
+   virtual ~IMEX_Evolution() { delete A_solver; delete A_prec; delete A; }
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   problem = 0;
+   const char *mesh_file = "../data/periodic-hexagon.mesh";
+   int ser_ref_levels = 2;
+   int par_ref_levels = 0;
+   int order = 3;
+   int ode_solver_type = 12;
+   double t_final = 10.0;
+   double d_coef = 0.01;
+   double dt = 0.01;
+   double sigma = -1.0;
+   double kappa = -1.0;
+   bool visualization = true;
+   bool visit = false;
+   bool binary = false;
+   int vis_steps = 5;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver: 1 - Forward Euler, 2 - RK2, 3 - RK3 SSP,"
+                  " 4 - RK4, 5 - Generalized Alpha,\n\t"
+                  "11 - Backward Euler, 12 - SDIRK2, 13 - SDIRK3,\n\t"
+                  "22 - Implicit Midpoint, 23 SDIRK23, 24 - SDIRK34,\n\t"
+                  "31 - IMEX BE/FE, 32 - IMEX RK2.");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&d_coef, "-d", "--diff-coef",
+                  "Diffusion coefficient.");
+   args.AddOption(&sigma, "-s", "--sigma",
+                  "One of the two DG penalty parameters, typically +1/-1."
+                  " See the documentation of class DGDiffusionIntegrator.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "One of the two DG penalty parameters, should be positive."
+                  " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
+                  "--no-visit-datafiles",
+                  "Save data files for VisIt (visit.llnl.gov) visualization.");
+   args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
+                  "--ascii-datafiles",
+                  "Use binary (Sidre) or ascii format for VisIt data files.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (kappa < 0)
+   {
+      kappa = (order+1)*(order+1);
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Define the ODE solver used for time integration. Several explicit,
+   //    implicitit, and implicit-explicit Runge-Kutta methods are available.
+   ODESolver *ode_solver = NULL;
+
+   switch (ode_solver_type)
+   {
+      // Explicit methods
+      case 1:  ode_solver = new ForwardEulerSolver; break;
+      case 2:  ode_solver = new RK2Solver(0.5); break; // midpoint method
+      case 3:  ode_solver = new RK3SSPSolver; break;
+      case 4:  ode_solver = new RK4Solver; break;
+      case 5:  ode_solver = new GeneralizedAlphaSolver(0.5); break;
+      // Implicit L-stable methods
+      case 11: ode_solver = new BackwardEulerSolver; break;
+      case 12: ode_solver = new SDIRK23Solver(2); break;
+      case 13: ode_solver = new SDIRK33Solver; break;
+      // Implicit A-stable methods (not L-stable)
+      case 22: ode_solver = new ImplicitMidpointSolver; break;
+      case 23: ode_solver = new SDIRK23Solver; break;
+      case 24: ode_solver = new SDIRK34Solver; break;
+      // Implicit-Explicit methods
+      case 31: ode_solver = new IMEX_BE_FE; break;
+      case 32: ode_solver = new IMEXRK2; break;
+      default:
+         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         return 3;
+   }
+
+   // 4. Read the serial mesh from the given mesh file on all processors. We can
+   //    handle geometrically periodic meshes in this code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 5. Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter. If the mesh is of NURBS type, we convert it
+   //    to a (piecewise-polynomial) high-order mesh.
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+   if (mesh->NURBSext)
+   {
+      mesh->SetCurvature(max(order, 1));
+   }
+   mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // 7. Define the parallel discontinuous DG finite element space on the
+   //    parallel refined mesh of the given polynomial order.
+   DG_FECollection fec(order, dim);
+   ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
+
+   HYPRE_Int global_vSize = fes->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of unknowns: " << global_vSize << endl;
+   }
+
+   // 8. Set up and assemble the parallel bilinear and linear forms (and the
+   //    parallel hypre matrices) corresponding to the DG discretization. The
+   //    DGTraceIntegrator involves integrals over mesh interior faces.
+   ConstantCoefficient diff_coef(d_coef);
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient u0(u0_function);
+
+   ParBilinearForm *m = new ParBilinearForm(fes);
+   m->AddDomainIntegrator(new MassIntegrator);
+
+   ParBilinearForm *s = new ParBilinearForm(fes);
+   s->AddDomainIntegrator(new DiffusionIntegrator(diff_coef));
+   s->AddInteriorFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma,
+                                                          kappa));
+   s->AddBdrFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma, kappa));
+
+   ParBilinearForm *k = new ParBilinearForm(fes);
+   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   k->AddInteriorFaceIntegrator(
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+   k->AddBdrFaceIntegrator(
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+
+   ParLinearForm *b = new ParLinearForm(fes);
+   b->AddBdrFaceIntegrator(
+      new DGDirichletLFIntegrator(u0, diff_coef, sigma, kappa));
+
+   int skip_zeros = 0;
+   m->Assemble(skip_zeros);
+   m->Finalize(skip_zeros);
+   s->Assemble(skip_zeros);
+   s->Finalize(skip_zeros);
+   k->Assemble(skip_zeros);
+   k->Finalize(skip_zeros);
+   b->Assemble();
+
+   HypreParMatrix *M = m->ParallelAssemble();
+   HypreParMatrix *S = s->ParallelAssemble();
+   HypreParMatrix *K = k->ParallelAssemble();
+   HypreParVector *B = b->ParallelAssemble();
+
+   // 9. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   ParGridFunction *u = new ParGridFunction(fes);
+   u->ProjectCoefficient(u0);
+   HypreParVector *U = u->GetTrueDofs();
+
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "ex23-mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "ex23-init." << setfill('0') << setw(6) << myid;
+      ofstream omesh(mesh_name.str().c_str());
+      omesh.precision(precision);
+      pmesh->Print(omesh);
+      ofstream osol(sol_name.str().c_str());
+      osol.precision(precision);
+      u->Save(osol);
+   }
+
+   // Create data collection for solution output: either VisItDataCollection for
+   // ascii data files, or SidreDataCollection for binary data files.
+   DataCollection *dc = NULL;
+   if (visit)
+   {
+      if (binary)
+      {
+#ifdef MFEM_USE_SIDRE
+         dc = new SidreDataCollection("Example23-Parallel", pmesh);
+#else
+         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+#endif
+      }
+      else
+      {
+         dc = new VisItDataCollection("Example23-Parallel", pmesh);
+         dc->SetPrecision(precision);
+         // To save the mesh using MFEM's parallel mesh format:
+         // dc->SetFormat(DataCollection::PARALLEL_FORMAT);
+      }
+      dc->RegisterField("solution", u);
+      dc->SetCycle(0);
+      dc->SetTime(0.0);
+      dc->Save();
+   }
+
+   socketstream sout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      if (!sout)
+      {
+         if (myid == 0)
+            cout << "Unable to connect to GLVis server at "
+                 << vishost << ':' << visport << endl;
+         visualization = false;
+         if (myid == 0)
+         {
+            cout << "GLVis visualization disabled.\n";
+         }
+      }
+      else
+      {
+         sout << "parallel " << num_procs << " " << myid << "\n";
+         sout.precision(precision);
+         sout << "solution\n" << *pmesh << *u;
+         sout << "pause\n";
+         sout << flush;
+         if (myid == 0)
+            cout << "GLVis visualization paused."
+                 << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+   // 10. Define the time-dependent evolution operator describing the ODE
+   //     right-hand side, and perform time-integration (looping over the time
+   //     iterations, ti, with a time-step dt).
+
+   TimeDependentOperator *adv = NULL;
+   if (ode_solver_type < 10)
+   {
+      adv = new EX_Evolution(*M, *S, *K, *B);
+   }
+   else if (ode_solver_type < 30)
+   {
+      adv = new IM_Evolution(*M, *S, *K, *B);
+   }
+   else
+   {
+      adv = new IMEX_Evolution(*M, *S, *K, *B);
+   }
+
+   double t = 0.0;
+   adv->SetTime(t);
+   ode_solver->Init(*adv);
+
+   int n_steps = (int)ceil(t_final / dt);
+   double dt_real = t_final / n_steps;
+
+   for (int ti = 0; ti < n_steps; )
+   {
+      ode_solver->Step(*U, t, dt_real);
+      ti++;
+
+      if (ti % vis_steps == 0 || ti == n_steps)
+      {
+         if (myid == 0)
+         {
+            cout << "time step: " << ti << ", time: " << t << endl;
+         }
+
+         // 11. Extract the parallel grid function corresponding to the finite
+         //     element approximation U (the local solution on each processor).
+         *u = *U;
+
+         if (visualization)
+         {
+            sout << "parallel " << num_procs << " " << myid << "\n";
+            sout << "solution\n" << *pmesh << *u << flush;
+         }
+
+         if (visit)
+         {
+            dc->SetCycle(ti);
+            dc->SetTime(t);
+            dc->Save();
+         }
+      }
+   }
+
+   // 12. Save the final solution in parallel. This output can be viewed later
+   //     using GLVis: "glvis -np <np> -m ex23-mesh -g ex23-final".
+   {
+      *u = *U;
+      ostringstream sol_name;
+      sol_name << "ex23-final." << setfill('0') << setw(6) << myid;
+      ofstream osol(sol_name.str().c_str());
+      osol.precision(precision);
+      u->Save(osol);
+   }
+
+   // 13. Free the used memory.
+   delete U;
+   delete u;
+   delete B;
+   delete b;
+   delete K;
+   delete k;
+   delete S;
+   delete s;
+   delete M;
+   delete m;
+   delete fes;
+   delete pmesh;
+   delete ode_solver;
+   delete adv;
+   delete dc;
+
+   MPI_Finalize();
+   return 0;
+}
+
+
+// Implementation of class EX_Evolution
+EX_Evolution::EX_Evolution(HypreParMatrix &_M, HypreParMatrix &_S,
+                           HypreParMatrix &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     M(_M), S(_S), K(_K), b(_b),
+     M_prec(M), M_solver(M.GetComm()), z(M.Height())
+{
+   M_prec.SetType(HypreSmoother::Jacobi);
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void EX_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (-S x + K x + b)
+   S.Mult(-1.0, x, 0.0, z);
+   K.Mult(1.0, x, 1.0, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+// Implementation of class IM_Evolution
+IM_Evolution::IM_Evolution(HypreParMatrix &_M, HypreParMatrix &_S,
+                           HypreParMatrix &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     M(_M), S(_S), K(_K), A(NULL), b(_b),
+     M_prec(M), M_solver(M.GetComm()),
+     A_prec(NULL), A_solver(NULL), dt(-1.0), z(M.Height())
+{
+   M_prec.SetType(HypreSmoother::Jacobi);
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void IM_Evolution::initA(double _dt)
+{
+   if (fabs(dt - _dt) > 1e-4 * _dt)
+   {
+      delete A_solver;
+      delete A_prec;
+      delete A;
+
+      HypreParMatrix * SK = Add(1.0, S, -1.0, K); // SK = S - K
+      A = Add(_dt, *SK, 1.0, M);                  // A = M + dt * (S - K)
+      delete SK;
+      dt = _dt;
+
+      A_prec = new HypreBoomerAMG(*A);
+      A_solver = new GMRESSolver(A->GetComm());
+      A_solver->SetOperator(*A);
+      A_solver->SetPreconditioner(*A_prec);
+
+      A_solver->iterative_mode = false;
+      A_solver->SetRelTol(1e-9);
+      A_solver->SetAbsTol(0.0);
+      A_solver->SetMaxIter(100);
+      A_solver->SetPrintLevel(0);
+   }
+}
+
+void IM_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (-S x + K x + b)
+   S.Mult(-1.0, x, 0.0, z);
+   K.Mult(1.0, x, 1.0, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IM_Evolution::ImplicitSolve(const double _dt, const Vector &x, Vector &y)
+{
+   this->initA(_dt);
+
+   // y = (M + dt S - dt K)^{-1} (-S x + K x + b)
+   S.Mult(-1.0, x, 0.0, z);
+   K.Mult(1.0, x, 1.0, z);
+   z += b;
+   A_solver->Mult(z, y);
+}
+
+// Implementation of class IMEX_Evolution
+IMEX_Evolution::IMEX_Evolution(HypreParMatrix &_M, HypreParMatrix &_S,
+                               HypreParMatrix &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     M(_M), S(_S), K(_K), A(NULL), b(_b),
+     M_prec(M), M_solver(M.GetComm()),
+     A_prec(NULL), A_solver(NULL), dt(-1.0), z(M.Height())
+{
+   M_prec.SetType(HypreSmoother::Jacobi);
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void IMEX_Evolution::initA(double _dt)
+{
+   if (fabs(dt - _dt) > 1e-4 * _dt)
+   {
+      delete A_solver;
+      delete A_prec;
+      delete A;
+
+      A = Add(_dt, S, 1.0, M); // A = M + dt * S
+      dt = _dt;
+
+      A_prec = new HypreBoomerAMG(*A);
+      A_solver = new CGSolver(A->GetComm());
+      A_solver->SetOperator(*A);
+      A_solver->SetPreconditioner(*A_prec);
+
+      A_solver->iterative_mode = false;
+      A_solver->SetRelTol(1e-9);
+      A_solver->SetAbsTol(0.0);
+      A_solver->SetMaxIter(100);
+      A_solver->SetPrintLevel(0);
+   }
+}
+
+void IMEX_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (-S x + K x + b)
+   S.Mult(-1.0, x, 0.0, z);
+   K.Mult(1.0, x, 1.0, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IMEX_Evolution::ExplicitMult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (K x + b)
+   K.Mult(1.0, x, 0.0, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IMEX_Evolution::ImplicitSolve(const double _dt, const Vector &x, Vector &y)
+{
+   this->initA(_dt);
+   // y = (M + dt S)^{-1} (-S x + b)
+   S.Mult(-1.0, x, 0.0, z);
+   z += b;
+   A_solver->Mult(z, y);
+}
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      {
+         // Translations in 1D, 2D, and 3D
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
+            case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
+               break;
+         }
+         break;
+      }
+      case 1:
+      case 2:
+      {
+         // Clockwise rotation in 2D around the origin
+         const double w = M_PI/2;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
+            case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+      case 3:
+      {
+         // Clockwise twisting rotation in 2D around the origin
+         const double w = M_PI/2;
+         double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
+         d = d*d;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
+            case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+   }
+}
+
+// Initial condition
+double u0_function(const Vector &x)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      {
+         switch (dim)
+         {
+            case 1:
+               return exp(-40.*pow(X(0)-0.5,2));
+            case 2:
+            case 3:
+            {
+               double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
+               if (dim == 3)
+               {
+                  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
+                  rx *= s;
+                  ry *= s;
+               }
+               return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
+                        erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
+            }
+         }
+      }
+      case 2:
+      {
+         double x_ = X(0), y_ = X(1), rho, phi;
+         rho = hypot(x_, y_);
+         phi = atan2(y_, x_);
+         return pow(sin(M_PI*rho),2)*sin(3*phi);
+      }
+      case 3:
+      {
+         const double f = M_PI;
+         return sin(f*X(0))*sin(f*X(1));
+      }
+   }
+   return 0.0;
+}

examples/ex3p_sp.cpp

@@ -1,334 +0,0 @@
-//                       MFEM Example 3 - Parallel Version
-//
-// Compile with: make ex3p
-//
-// Sample runs:  mpirun -np 4 ex3p -m ../data/star.mesh
-//               mpirun -np 4 ex3p -m ../data/square-disc.mesh -o 2
-//               mpirun -np 4 ex3p -m ../data/beam-tet.mesh
-//               mpirun -np 4 ex3p -m ../data/beam-hex.mesh
-//               mpirun -np 4 ex3p -m ../data/escher.mesh
-//               mpirun -np 4 ex3p -m ../data/escher.mesh -o 2
-//               mpirun -np 4 ex3p -m ../data/fichera.mesh
-//               mpirun -np 4 ex3p -m ../data/fichera-q2.vtk
-//               mpirun -np 4 ex3p -m ../data/fichera-q3.mesh
-//               mpirun -np 4 ex3p -m ../data/square-disc-nurbs.mesh
-//               mpirun -np 4 ex3p -m ../data/beam-hex-nurbs.mesh
-//               mpirun -np 4 ex3p -m ../data/amr-quad.mesh -o 2
-//               mpirun -np 4 ex3p -m ../data/amr-hex.mesh
-//               mpirun -np 4 ex3p -m ../data/star-surf.mesh -o 2
-//               mpirun -np 4 ex3p -m ../data/mobius-strip.mesh -o 2 -f 0.1
-//               mpirun -np 4 ex3p -m ../data/klein-bottle.mesh -o 2 -f 0.1
-//
-// Description:  This example code solves a simple electromagnetic diffusion
-//               problem corresponding to the second order definite Maxwell
-//               equation curl curl E + E = f with boundary condition
-//               E x n = <given tangential field>. Here, we use a given exact
-//               solution E and compute the corresponding r.h.s. f.
-//               We discretize with Nedelec finite elements in 2D or 3D.
-//
-//               The example demonstrates the use of H(curl) finite element
-//               spaces with the curl-curl and the (vector finite element) mass
-//               bilinear form, as well as the computation of discretization
-//               error when the exact solution is known. Static condensation is
-//               also illustrated.
-//
-//               We recommend viewing examples 1-2 before viewing this example.
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-// Exact solution, E, and r.h.s., f. See below for implementation.
-void E_exact(const Vector &, Vector &);
-void f_exact(const Vector &, Vector &);
-double freq = 1.0, kappa;
-int dim;
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI.
-   int num_procs, myid;
-   MPI_Init(&argc, &argv);
-   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
-   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
-
-   // 2. Parse command-line options.
-   const char *mesh_file = "../data/beam-tet.mesh";
-   int order = 1;
-   bool static_cond = false;
-   bool visualization = 1;
-#ifdef MFEM_USE_STRUMPACK
-   bool use_strumpack = false;
-#endif
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree).");
-   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
-                  " solution.");
-   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                  "--no-static-condensation", "Enable static condensation.");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-#ifdef MFEM_USE_STRUMPACK
-   args.AddOption(&use_strumpack, "-strumpack", "--strumpack-solver",
-                  "-no-strumpack", "--no-strumpack-solver",
-                  "Use STRUMPACK's double complex linear solver.");
-#endif
-
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      MPI_Finalize();
-      return 1;
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
-   kappa = freq * M_PI;
-
-   // 3. Read the (serial) mesh from the given mesh file on all processors.  We
-   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
-   //    and volume meshes with the same code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   dim = mesh->Dimension();
-   int sdim = mesh->SpaceDimension();
-
-   // 4. Refine the serial mesh on all processors to increase the resolution. In
-   //    this example we do 'ref_levels' of uniform refinement. We choose
-   //    'ref_levels' to be the largest number that gives a final mesh with no
-   //    more than 1,000 elements.
-   {
-      int ref_levels =
-         (int)floor(log(100000./mesh->GetNE())/log(2.)/dim);
-      for (int l = 0; l < ref_levels; l++)
-      {
-         mesh->UniformRefinement();
-      }
-   }
-
-   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh further in parallel to increase the resolution. Once the
-   //    parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
-   //    meshes need to be reoriented before we can define high-order Nedelec
-   //    spaces on them.
-   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
-   delete mesh;
-   {
-      int par_ref_levels = 2;
-      for (int l = 0; l < par_ref_levels; l++)
-      {
-         pmesh->UniformRefinement();
-      }
-   }
-   pmesh->ReorientTetMesh();
-
-   // 6. Define a parallel finite element space on the parallel mesh. Here we
-   //    use the Nedelec finite elements of the specified order.
-   FiniteElementCollection *fec = new ND_FECollection(order, dim);
-   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
-   HYPRE_Int size = fespace->GlobalTrueVSize();
-   if (myid == 0)
-   {
-      cout << "Number of finite element unknowns: " << size << endl;
-   }
-
-   // 7. Determine the list of true (i.e. parallel conforming) essential
-   //    boundary dofs. In this example, the boundary conditions are defined
-   //    by marking all the boundary attributes from the mesh as essential
-   //    (Dirichlet) and converting them to a list of true dofs.
-   Array<int> ess_tdof_list;
-   if (pmesh->bdr_attributes.Size())
-   {
-      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
-      ess_bdr = 1;
-      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-   }
-
-   // 8. Set up the parallel linear form b(.) which corresponds to the
-   //    right-hand side of the FEM linear system, which in this case is
-   //    (f,phi_i) where f is given by the function f_exact and phi_i are the
-   //    basis functions in the finite element fespace.
-   VectorFunctionCoefficient f(sdim, f_exact);
-   ParLinearForm *b = new ParLinearForm(fespace);
-   b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
-   b->Assemble();
-
-   // 9. Define the solution vector x as a parallel finite element grid function
-   //    corresponding to fespace. Initialize x by projecting the exact
-   //    solution. Note that only values from the boundary edges will be used
-   //    when eliminating the non-homogeneous boundary condition to modify the
-   //    r.h.s. vector b.
-   ParGridFunction x(fespace);
-   VectorFunctionCoefficient E(sdim, E_exact);
-   x.ProjectCoefficient(E);
-
-   // 10. Set up the parallel bilinear form corresponding to the EM diffusion
-   //     operator curl muinv curl + sigma I, by adding the curl-curl and the
-   //     mass domain integrators.
-   Coefficient *muinv = new ConstantCoefficient(1.0);
-   Coefficient *sigma = new ConstantCoefficient(-1.0);
-   ParBilinearForm *a = new ParBilinearForm(fespace);
-   a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
-   a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma));
-
-   // 11. Assemble the parallel bilinear form and the corresponding linear
-   //     system, applying any necessary transformations such as: parallel
-   //     assembly, eliminating boundary conditions, applying conforming
-   //     constraints for non-conforming AMR, static condensation, etc.
-   if (static_cond) { a->EnableStaticCondensation(); }
-   a->Assemble();
-
-   HypreParMatrix A;
-   Vector B, X;
-   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
-
-   if (myid == 0)
-   {
-      cout << "Size of linear system: " << A.GetGlobalNumRows() << endl;
-   }
-
-   StopWatch chrono;
-   chrono.Clear();
-   chrono.Start();
-
-#ifdef MFEM_USE_STRUMPACK
-   if (use_strumpack)
-   {
-      Operator * Arow = new STRUMPACKRowLocMatrix(A);
-
-      STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv, MPI_COMM_WORLD);
-      strumpack->SetPrintFactorStatistics(true);
-      strumpack->SetPrintSolveStatistics(false);
-      strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
-      strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
-      // strumpack->SetMC64Job(strumpack::MC64Job::NONE);
-      // strumpack->SetSymmetricPattern(true);
-      strumpack->SetOperator(*Arow);
-      strumpack->SetFromCommandLine();
-      //Solver * precond = strumpack;
-
-      strumpack->Mult(B, X);
-
-      delete strumpack;
-      delete Arow;
-   }
-   else
-#endif
-   {
-      // 12. Define and apply a parallel PCG solver for AX=B with the AMS
-      //     preconditioner from hypre.
-      ParFiniteElementSpace *prec_fespace =
-         (a->StaticCondensationIsEnabled() ? a->SCParFESpace() : fespace);
-      HypreSolver *ams = new HypreAMS(A, prec_fespace);
-      HyprePCG *pcg = new HyprePCG(A);
-      pcg->SetTol(1e-12);
-      pcg->SetMaxIter(500);
-      pcg->SetPrintLevel(2);
-      pcg->SetPreconditioner(*ams);
-      pcg->Mult(B, X);
-
-      delete pcg;
-      delete ams;
-   }
-
-   chrono.Stop();
-   cout << "Solver time " << chrono.RealTime() << endl;
-
-   // 13. Recover the parallel grid function corresponding to X. This is the
-   //     local finite element solution on each processor.
-   a->RecoverFEMSolution(X, *b, x);
-
-   // 14. Compute and print the L^2 norm of the error.
-   {
-      double err = x.ComputeL2Error(E);
-      if (myid == 0)
-      {
-         cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
-      }
-   }
-
-   // 15. Save the refined mesh and the solution in parallel. This output can
-   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
-   {
-      ostringstream mesh_name, sol_name;
-      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
-      sol_name << "sol." << setfill('0') << setw(6) << myid;
-
-      ofstream mesh_ofs(mesh_name.str().c_str());
-      mesh_ofs.precision(8);
-      pmesh->Print(mesh_ofs);
-
-      ofstream sol_ofs(sol_name.str().c_str());
-      sol_ofs.precision(8);
-      x.Save(sol_ofs);
-   }
-
-   // 16. Send the solution by socket to a GLVis server.
-   if (visualization)
-   {
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      socketstream sol_sock(vishost, visport);
-      sol_sock << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock.precision(8);
-      sol_sock << "solution\n" << *pmesh << x << flush;
-   }
-
-   // 17. Free the used memory.
-   delete a;
-   delete sigma;
-   delete muinv;
-   delete b;
-   delete fespace;
-   delete fec;
-   delete pmesh;
-
-   MPI_Finalize();
-
-   return 0;
-}
-
-
-void E_exact(const Vector &x, Vector &E)
-{
-   if (dim == 3)
-   {
-      E(0) = sin(kappa * x(1));
-      E(1) = sin(kappa * x(2));
-      E(2) = sin(kappa * x(0));
-   }
-   else
-   {
-      E(0) = sin(kappa * x(1));
-      E(1) = sin(kappa * x(0));
-      if (x.Size() == 3) { E(2) = 0.0; }
-   }
-}
-
-void f_exact(const Vector &x, Vector &f)
-{
-   if (dim == 3)
-   {
-      f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
-      f(1) = (1. + kappa * kappa) * sin(kappa * x(2));
-      f(2) = (1. + kappa * kappa) * sin(kappa * x(0));
-   }
-   else
-   {
-      f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
-      f(1) = (1. + kappa * kappa) * sin(kappa * x(0));
-      if (x.Size() == 3) { f(2) = 0.0; }
-   }
-}

examples/ex6.cpp

@@ -16,9 +16,9 @@
 //               ex6 -m ../data/amr-quad.mesh
 //
 // Device sample runs:
-//             > ex6 -pa -d cuda
-//             > ex6 -pa -d occa-cuda
-//             > ex6 -pa -d raja-omp
+//               ex6 -pa -d cuda
+//               ex6 -pa -d occa-cuda
+//               ex6 -pa -d raja-omp
 //
 // Description:  This is a version of Example 1 with a simple adaptive mesh
 //               refinement loop. The problem being solved is again the Laplace

examples/ex6p.cpp

@@ -16,9 +16,9 @@
 //               mpirun -np 4 ex6p -m ../data/amr-quad.mesh
 //
 // Device sample runs:
-//             > mpirun -np 4 ex6p -pa -d cuda
-//             > mpirun -np 4 ex6p -pa -d occa-cuda
-//             > mpirun -np 4 ex6p -pa -d raja-omp
+//               mpirun -np 4 ex6p -pa -d cuda
+//               mpirun -np 4 ex6p -pa -d occa-cuda
+//               mpirun -np 4 ex6p -pa -d raja-omp
 //
 // Description:  This is a version of Example 1 with a simple adaptive mesh
 //               refinement loop. The problem being solved is again the Laplace

examples/makefile

@@ -22,9 +22,9 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_EXAMPLES = ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 ex17\
-  ex18 ex19 ex20 ex21 ex22
+  ex18 ex19 ex20 ex22 ex23
 PAR_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p ex12p\
- ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p
+ ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex22p ex23p
 
 ifeq ($(MFEM_USE_MPI),NO)
    EXAMPLES = $(SEQ_EXAMPLES)
@@ -117,8 +117,8 @@ clean-build:
 
 clean-exec:
 	@rm -f refined.mesh displaced.mesh mesh.* ex5.mesh
-	@rm -rf Example5* Example9* Example15* Example16*
-	@rm -f sphere_refined.* sol.* sol_u.* sol_p.* sol_r.* sol_i.*
+	@rm -rf Example5* Example9* Example15* Example16* Example23*
+	@rm -f sphere_refined.* sol.* sol_u.* sol_p.*
 	@rm -f ex9.mesh ex9-mesh.* ex9-init.* ex9-final.*
 	@rm -f deformed.* velocity.* elastic_energy.* mode_*
 	@rm -f ex16.mesh ex16-mesh.* ex16-init.* ex16-final.*
@@ -126,3 +126,4 @@ clean-exec:
 	@rm -f deformation.* pressure.*
 	@rm -f ex20.dat ex20p_?????.dat gnuplot_ex20.inp gnuplot_ex20p.inp
 	@rm -f ex22*.mesh ex22*.sol ex22p_*.*
+	@rm -f ex23.mesh ex23-mesh.* ex23-init.* ex23-final.*

fem/bilinearform.cpp

@@ -588,14 +588,18 @@ void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
                                     Vector &b, OperatorHandle &A, Vector &X,
                                     Vector &B, int copy_interior)
 {
+   const SparseMatrix *P = fes->GetConformingProlongation();
+
    if (ext)
    {
+      if (P != NULL && assembly != AssemblyLevel::FULL && Device::IsEnabled())
+      {
+         P->BuildTranspose();
+      }
       ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
       return;
    }
 
-   const SparseMatrix *P = fes->GetConformingProlongation();
-
    FormSystemMatrix(ess_tdof_list, A);
 
    // Transform the system and perform the elimination in B, based on the

fem/complex_fem.cpp

@@ -1,784 +0,0 @@
-// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
-// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
-// reserved. See file COPYRIGHT for details.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability see http://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the GNU Lesser General Public License (as published by the Free
-// Software Foundation) version 2.1 dated February 1999.
-
-#include "complex_fem.hpp"
-
-using namespace std;
-
-namespace mfem
-{
-
-ComplexGridFunction::ComplexGridFunction(FiniteElementSpace *fes)
-   : Vector(2*(fes->GetVSize()))
-{
-   gfr_ = new GridFunction(fes, &data[0]);
-   gfi_ = new GridFunction(fes, &data[fes->GetVSize()]);
-}
-
-void
-ComplexGridFunction::Update()
-{
-   FiniteElementSpace * fes = gfr_->FESpace();
-
-   int vsize = fes->GetVSize();
-
-   const Operator *T = fes->GetUpdateOperator();
-   if (T)
-   {
-      // Update the individual GridFunction objects.  This will allocate
-      // new data arrays for each GridFunction.
-      gfr_->Update();
-      gfi_->Update();
-
-      // Our data array now contains old data as well as being the wrong size
-      // so reallocate it.
-      this->SetSize(2 * vsize);
-
-      // Create temporary vectors which point to the new data array
-      Vector gf_r(&data[0], vsize);
-      Vector gf_i(&data[vsize], vsize);
-
-      // Copy the updated GridFunctions into the new data array
-      gf_r = *gfr_;
-      gf_i = *gfi_;
-
-      // Replace the individual data arrays with pointers into the new data array
-      gfr_->NewDataAndSize(&data[0], vsize);
-      gfi_->NewDataAndSize(&data[vsize], vsize);
-   }
-   else
-   {
-      // The existing data will not be transferred to the new GridFunctions
-      // so delete it a allocate a new array
-      this->SetSize(2 * vsize);
-
-      // Point the individual GridFunctions to the new data array
-      gfr_->NewDataAndSize(&data[0], vsize);
-      gfi_->NewDataAndSize(&data[vsize], vsize);
-
-      // These updates will only set the proper 'sequence' value within
-      // the individual GridFunction objects because their sizes are
-      // already correct
-      gfr_->Update();
-      gfi_->Update();
-   }
-}
-
-void
-ComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
-                                        Coefficient &imag_coeff)
-{
-   gfr_->ProjectCoefficient(real_coeff);
-   gfi_->ProjectCoefficient(imag_coeff);
-}
-
-void
-ComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
-                                        VectorCoefficient &imag_vcoeff)
-{
-   gfr_->ProjectCoefficient(real_vcoeff);
-   gfi_->ProjectCoefficient(imag_vcoeff);
-}
-
-void
-ComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
-                                           Coefficient &imag_coeff,
-                                           Array<int> &attr)
-{
-   gfr_->ProjectBdrCoefficient(real_coeff, attr);
-   gfi_->ProjectBdrCoefficient(imag_coeff, attr);
-}
-
-void
-ComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
-                                                 VectorCoefficient &imag_vcoeff,
-                                                 Array<int> &attr)
-{
-   gfr_->ProjectBdrCoefficientNormal(real_vcoeff, attr);
-   gfi_->ProjectBdrCoefficientNormal(imag_vcoeff, attr);
-}
-
-void
-ComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
-                                                  &real_vcoeff,
-                                                  VectorCoefficient
-                                                  &imag_vcoeff,
-                                                  Array<int> &attr)
-{
-   gfr_->ProjectBdrCoefficientTangent(real_vcoeff, attr);
-   gfi_->ProjectBdrCoefficientTangent(imag_vcoeff, attr);
-}
-
-
-ComplexLinearForm::ComplexLinearForm(FiniteElementSpace *f,
-                                     ComplexOperator::Convention convention)
-   : Vector(2*(f->GetVSize())),
-     conv_(convention)
-{
-   lfr_ = new LinearForm(f, &data[0]);
-   lfi_ = new LinearForm(f, &data[f->GetVSize()]);
-}
-
-ComplexLinearForm::~ComplexLinearForm()
-{
-   delete lfr_;
-   delete lfi_;
-}
-
-void
-ComplexLinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi_real,
-                                       LinearFormIntegrator *lfi_imag)
-{
-   if ( lfi_real ) { lfr_->AddDomainIntegrator(lfi_real); }
-   if ( lfi_imag ) { lfi_->AddDomainIntegrator(lfi_imag); }
-}
-
-void
-ComplexLinearForm::Update()
-{
-   FiniteElementSpace *fes = lfr_->FESpace();
-
-   this->Update(fes);
-}
-
-void
-ComplexLinearForm::Update(FiniteElementSpace *fes)
-{
-   int vsize = fes->GetVSize();
-   SetSize(2 * vsize);
-
-   Vector lfr(&data[0], vsize);
-   Vector lfi(&data[vsize], vsize);
-
-   lfr_->Update(fes, lfr, 0);
-   lfi_->Update(fes, lfi, 0);
-}
-
-void
-ComplexLinearForm::Assemble()
-{
-   lfr_->Assemble();
-   lfi_->Assemble();
-   if (conv_ == ComplexOperator::BLOCK_SYMMETRIC)
-   {
-      *lfi_ *= -1.0;
-   }
-}
-
-complex<double>
-ComplexLinearForm::operator()(const ComplexGridFunction &gf) const
-{
-   double s = (conv_ == ComplexOperator::HERMITIAN)?1.0:-1.0;
-   return complex<double>((*lfr_)(gf.real()) - s * (*lfi_)(gf.imag()),
-                          (*lfr_)(gf.imag()) + s * (*lfi_)(gf.real()));
-}
-
-
-SesquilinearForm::SesquilinearForm(FiniteElementSpace *f,
-                                   ComplexOperator::Convention convention)
-   : conv_(convention),
-     blfr_(new BilinearForm(f)),
-     blfi_(new BilinearForm(f))
-{}
-
-SesquilinearForm::~SesquilinearForm()
-{
-   delete blfr_;
-   delete blfi_;
-}
-
-void SesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
-                                           BilinearFormIntegrator *bfi_imag)
-{
-   if (bfi_real) { blfr_->AddDomainIntegrator(bfi_real); }
-   if (bfi_imag) { blfi_->AddDomainIntegrator(bfi_imag); }
-}
-
-void
-SesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
-                                        BilinearFormIntegrator *bfi_imag)
-{
-   if (bfi_real) { blfr_->AddBoundaryIntegrator(bfi_real); }
-   if (bfi_imag) { blfi_->AddBoundaryIntegrator(bfi_imag); }
-}
-
-void
-SesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
-                                        BilinearFormIntegrator *bfi_imag,
-                                        Array<int> & bdr_marker)
-{
-   if (bfi_real) { blfr_->AddBoundaryIntegrator(bfi_real, bdr_marker); }
-   if (bfi_imag) { blfi_->AddBoundaryIntegrator(bfi_imag, bdr_marker); }
-}
-
-void
-SesquilinearForm::Assemble(int skip_zeros)
-{
-   blfr_->Assemble(skip_zeros);
-   blfi_->Assemble(skip_zeros);
-}
-
-void
-SesquilinearForm::Finalize(int skip_zeros)
-{
-   blfr_->Finalize(skip_zeros);
-   blfi_->Finalize(skip_zeros);
-}
-
-ComplexSparseMatrix *
-SesquilinearForm::AssembleCompSpMat()
-{
-   return new ComplexSparseMatrix(&blfr_->SpMat(),
-                                  &blfi_->SpMat(),
-                                  false, false, conv_);
-
-}
-
-void
-SesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
-                                   Vector &x, Vector &b,
-                                   OperatorHandle &A,
-                                   Vector &X, Vector &B,
-                                   int ci)
-{
-   FiniteElementSpace * fes = blfr_->FESpace();
-
-   int vsize  = fes->GetVSize();
-   // int tvsize = pfes->GetTrueVSize();
-
-   double s = (conv_ == ComplexOperator::HERMITIAN)?1.0:-1.0;
-
-   // Allocate temporary vectors
-   Vector b_0(vsize);  b_0 = 0.0;
-   // Vector B_0(tvsize); B_0 = 0.0;
-
-   // Extract the real and imaginary parts of the input vectors
-   MFEM_ASSERT(x.Size() == 2 * vsize, "Input GridFunction of incorrect size!");
-   Vector x_r(x.GetData(), vsize);
-   Vector x_i(&(x.GetData())[vsize], vsize);
-
-   MFEM_ASSERT(b.Size() == 2 * vsize, "Input LinearForm of incorrect size!");
-   Vector b_r(b.GetData(), vsize);
-   Vector b_i(&(b.GetData())[vsize], vsize);
-   b_i *= s;
-   /*
-   X.SetSize(2 * tvsize);
-   Vector X_r(X.GetData(), tvsize);
-   Vector X_i(&(X.GetData())[tvsize], tvsize);
-
-   B.SetSize(2 * tvsize);
-   Vector B_r(B.GetData(), tvsize);
-   Vector B_i(&(B.GetData())[tvsize], tvsize);
-   */
-   SparseMatrix * A_r = new SparseMatrix;
-   SparseMatrix * A_i = new SparseMatrix;
-   Vector X_0, B_0;
-
-   b_0 = b_r;
-   blfr_->FormLinearSystem(ess_tdof_list, x_r, b_r, *A_r, X_0, B_0, ci);
-
-   int tvsize = B_0.Size();
-   X.SetSize(2 * tvsize);
-   B.SetSize(2 * tvsize);
-   Vector X_r(X.GetData(), tvsize);
-   Vector X_i(&(X.GetData())[tvsize], tvsize);
-   Vector B_r(B.GetData(), tvsize);
-   Vector B_i(&(B.GetData())[tvsize], tvsize);
-   X_r = X_0; B_r = B_0;
-
-   b_0 = 0.0;
-   blfi_->FormLinearSystem(ess_tdof_list, x_i, b_0, *A_i, X_0, B_0, false);
-   B_r -= B_0;
-
-   b_0 = b_i;
-   blfr_->FormLinearSystem(ess_tdof_list, x_i, b_0, *A_r, X_0, B_0, ci);
-   X_i = X_0; B_i = B_0;
-
-   b_0 = 0.0;
-   blfi_->FormLinearSystem(ess_tdof_list, x_r, b_0, *A_i, X_0, B_0, false);
-   B_i += B_0;
-
-   B_i *= s;
-   b_i *= s;
-
-   // A = A_r + i A_i
-   A.Clear();
-   ComplexSparseMatrix * A_sp =
-      new ComplexSparseMatrix(A_r, A_i, true, true, conv_);
-   A.Reset<ComplexSparseMatrix>(A_sp, true);
-}
-
-void
-SesquilinearForm::RecoverFEMSolution(const Vector &X, const Vector &b,
-                                     Vector &x)
-{
-   FiniteElementSpace * fes = blfr_->FESpace();
-
-   const SparseMatrix *P = fes->GetConformingProlongation();
-
-   int vsize  = fes->GetVSize();
-   int tvsize = X.Size() / 2;
-
-   Vector X_r(X.GetData(), tvsize);
-   Vector X_i(&(X.GetData())[tvsize], tvsize);
-
-   Vector x_r(x.GetData(), vsize);
-   Vector x_i(&(x.GetData())[vsize], vsize);
-
-   if (!P)
-   {
-      x = X;
-   }
-   else
-   {
-      // Apply conforming prolongation
-      P->Mult(X_r, x_r);
-      P->Mult(X_i, x_i);
-   }
-}
-
-void
-SesquilinearForm::Update(FiniteElementSpace *nfes)
-{
-   if ( blfr_ ) { blfr_->Update(nfes); }
-   if ( blfi_ ) { blfi_->Update(nfes); }
-}
-
-
-#ifdef MFEM_USE_MPI
-
-ParComplexGridFunction::ParComplexGridFunction(ParFiniteElementSpace *pfes)
-   : Vector(2*(pfes->GetVSize()))
-{
-   pgfr_ = new ParGridFunction(pfes, &data[0]);
-   pgfi_ = new ParGridFunction(pfes, &data[pfes->GetVSize()]);
-}
-
-void
-ParComplexGridFunction::Update()
-{
-   ParFiniteElementSpace * pfes = pgfr_->ParFESpace();
-
-   int vsize = pfes->GetVSize();
-
-   const Operator *T = pfes->GetUpdateOperator();
-   if (T)
-   {
-      // Update the individual GridFunction objects.  This will allocate
-      // new data arrays for each GridFunction.
-      pgfr_->Update();
-      pgfi_->Update();
-
-      // Our data array now contains old data as well as being the wrong size
-      // so reallocate it.
-      this->SetSize(2 * vsize);
-
-      // Create temporary vectors which point to the new data array
-      Vector gf_r(&data[0], vsize);
-      Vector gf_i(&data[vsize], vsize);
-
-      // Copy the updated GridFunctions into the new data array
-      gf_r = *pgfr_;
-      gf_i = *pgfi_;
-
-      // Replace the individual data arrays with pointers into the new data array
-      pgfr_->NewDataAndSize(&data[0], vsize);
-      pgfi_->NewDataAndSize(&data[vsize], vsize);
-   }
-   else
-   {
-      // The existing data will not be transferred to the new GridFunctions
-      // so delete it a allocate a new array
-      this->SetSize(2 * vsize);
-
-      // Point the individual GridFunctions to the new data array
-      pgfr_->NewDataAndSize(&data[0], vsize);
-      pgfi_->NewDataAndSize(&data[vsize], vsize);
-
-      // These updates will only set the proper 'sequence' value within
-      // the individual GridFunction objects because their sizes are
-      // already correct
-      pgfr_->Update();
-      pgfi_->Update();
-   }
-}
-
-void
-ParComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
-                                           Coefficient &imag_coeff)
-{
-   pgfr_->ProjectCoefficient(real_coeff);
-   pgfi_->ProjectCoefficient(imag_coeff);
-}
-
-void
-ParComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
-                                           VectorCoefficient &imag_vcoeff)
-{
-   pgfr_->ProjectCoefficient(real_vcoeff);
-   pgfi_->ProjectCoefficient(imag_vcoeff);
-}
-
-void
-ParComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
-                                              Coefficient &imag_coeff,
-                                              Array<int> &attr)
-{
-   pgfr_->ProjectBdrCoefficient(real_coeff, attr);
-   pgfi_->ProjectBdrCoefficient(imag_coeff, attr);
-}
-
-void
-ParComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
-                                                    &real_vcoeff,
-                                                    VectorCoefficient
-                                                    &imag_vcoeff,
-                                                    Array<int> &attr)
-{
-   pgfr_->ProjectBdrCoefficientNormal(real_vcoeff, attr);
-   pgfi_->ProjectBdrCoefficientNormal(imag_vcoeff, attr);
-}
-
-void
-ParComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
-                                                     &real_vcoeff,
-                                                     VectorCoefficient
-                                                     &imag_vcoeff,
-                                                     Array<int> &attr)
-{
-   pgfr_->ProjectBdrCoefficientTangent(real_vcoeff, attr);
-   pgfi_->ProjectBdrCoefficientTangent(imag_vcoeff, attr);
-}
-
-void
-ParComplexGridFunction::Distribute(const Vector *tv)
-{
-   ParFiniteElementSpace * pfes = pgfr_->ParFESpace();
-   HYPRE_Int size = pfes->GetTrueVSize();
-
-   double * tvd = tv->GetData();
-   Vector tvr(tvd, size);
-   Vector tvi(&tvd[size], size);
-
-   pgfr_->Distribute(tvr);
-   pgfi_->Distribute(tvi);
-}
-
-void
-ParComplexGridFunction::ParallelProject(Vector &tv) const
-{
-   ParFiniteElementSpace * pfes = pgfr_->ParFESpace();
-   HYPRE_Int size = pfes->GetTrueVSize();
-
-   double * tvd = tv.GetData();
-   Vector tvr(tvd, size);
-   Vector tvi(&tvd[size], size);
-
-   pgfr_->ParallelProject(tvr);
-   pgfi_->ParallelProject(tvi);
-}
-
-
-ParComplexLinearForm::ParComplexLinearForm(ParFiniteElementSpace *pfes,
-                                           ComplexOperator::Convention
-                                           convention)
-   : Vector(2*(pfes->GetVSize())),
-     conv_(convention)
-{
-   plfr_ = new ParLinearForm(pfes, &data[0]);
-   plfi_ = new ParLinearForm(pfes, &data[pfes->GetVSize()]);
-
-   HYPRE_Int * tdof_offsets = pfes->GetTrueDofOffsets();
-
-   int n = (HYPRE_AssumedPartitionCheck()) ? 2 : pfes->GetNRanks();
-   tdof_offsets_ = new HYPRE_Int[n+1];
-
-   for (int i=0; i<=n; i++)
-   {
-      tdof_offsets_[i] = 2 * tdof_offsets[i];
-   }
-}
-
-ParComplexLinearForm::~ParComplexLinearForm()
-{
-   delete plfr_;
-   delete plfi_;
-   delete [] tdof_offsets_;
-}
-
-void
-ParComplexLinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi_real,
-                                          LinearFormIntegrator *lfi_imag)
-{
-   if ( lfi_real ) { plfr_->AddDomainIntegrator(lfi_real); }
-   if ( lfi_imag ) { plfi_->AddDomainIntegrator(lfi_imag); }
-}
-
-void
-ParComplexLinearForm::Update(ParFiniteElementSpace *pf)
-{
-   ParFiniteElementSpace *pfes = (pf!=NULL)?pf:plfr_->ParFESpace();
-   int vsize = pfes->GetVSize();
-   SetSize(2 * vsize);
-
-   Vector plfr(&data[0], vsize);
-   Vector plfi(&data[vsize], vsize);
-
-   plfr_->Update(pfes, plfr, 0);
-   plfi_->Update(pfes, plfi, 0);
-}
-
-void
-ParComplexLinearForm::Assemble()
-{
-   plfr_->Assemble();
-   plfi_->Assemble();
-   if (conv_ == ComplexOperator::BLOCK_SYMMETRIC)
-   {
-      *plfi_ *= -1.0;
-   }
-}
-
-void
-ParComplexLinearForm::ParallelAssemble(Vector &tv)
-{
-   HYPRE_Int size = plfr_->ParFESpace()->GetTrueVSize();
-
-   double * tvd = tv.GetData();
-   Vector tvr(tvd, size);
-   Vector tvi(&tvd[size], size);
-
-   plfr_->ParallelAssemble(tvr);
-   plfi_->ParallelAssemble(tvi);
-}
-
-HypreParVector *
-ParComplexLinearForm::ParallelAssemble()
-{
-   const ParFiniteElementSpace * pfes = plfr_->ParFESpace();
-
-   HypreParVector * tv = new HypreParVector(pfes->GetComm(),
-                                            2*(pfes->GlobalTrueVSize()),
-                                            tdof_offsets_);
-
-   HYPRE_Int size = pfes->GetTrueVSize();
-
-   double * tvd = tv->GetData();
-   Vector tvr(tvd, size);
-   Vector tvi(&tvd[size], size);
-
-   plfr_->ParallelAssemble(tvr);
-   plfi_->ParallelAssemble(tvi);
-
-   return tv;
-}
-
-complex<double>
-ParComplexLinearForm::operator()(const ParComplexGridFunction &gf) const
-{
-   double s = (conv_ == ComplexOperator::HERMITIAN)?1.0:-1.0;
-   return complex<double>((*plfr_)(gf.real()) - s * (*plfi_)(gf.imag()),
-                          (*plfr_)(gf.imag()) + s * (*plfi_)(gf.real()));
-}
-
-
-ParSesquilinearForm::ParSesquilinearForm(ParFiniteElementSpace *pf,
-                                         ComplexOperator::Convention
-                                         convention)
-   : conv_(convention),
-     pblfr_(new ParBilinearForm(pf)),
-     pblfi_(new ParBilinearForm(pf))
-{}
-
-ParSesquilinearForm::~ParSesquilinearForm()
-{
-   delete pblfr_;
-   delete pblfi_;
-}
-
-void ParSesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
-                                              BilinearFormIntegrator *bfi_imag)
-{
-   if (bfi_real) { pblfr_->AddDomainIntegrator(bfi_real); }
-   if (bfi_imag) { pblfi_->AddDomainIntegrator(bfi_imag); }
-}
-
-void
-ParSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
-                                           BilinearFormIntegrator *bfi_imag)
-{
-   if (bfi_real) { pblfr_->AddBoundaryIntegrator(bfi_real); }
-   if (bfi_imag) { pblfi_->AddBoundaryIntegrator(bfi_imag); }
-}
-
-void
-ParSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
-                                           BilinearFormIntegrator *bfi_imag,
-                                           Array<int> & bdr_marker)
-{
-   if (bfi_real) { pblfr_->AddBoundaryIntegrator(bfi_real, bdr_marker); }
-   if (bfi_imag) { pblfi_->AddBoundaryIntegrator(bfi_imag, bdr_marker); }
-}
-
-void
-ParSesquilinearForm::Assemble(int skip_zeros)
-{
-   pblfr_->Assemble(skip_zeros);
-   pblfi_->Assemble(skip_zeros);
-}
-
-void
-ParSesquilinearForm::Finalize(int skip_zeros)
-{
-   pblfr_->Finalize(skip_zeros);
-   pblfi_->Finalize(skip_zeros);
-}
-
-ComplexHypreParMatrix *
-ParSesquilinearForm::ParallelAssemble()
-{
-   return new ComplexHypreParMatrix(pblfr_->ParallelAssemble(),
-                                    pblfi_->ParallelAssemble(),
-                                    true, true, conv_);
-
-}
-
-void
-ParSesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
-                                      Vector &x, Vector &b,
-                                      OperatorHandle &A,
-                                      Vector &X, Vector &B,
-                                      int ci)
-{
-   ParFiniteElementSpace * pfes = pblfr_->ParFESpace();
-
-   int tvs = pfes->TrueVSize();
-   cout << "TrueVSize returns " << tvs << endl;
-   cout << "GetVSize returns " << pfes->GetVSize() << endl;
-
-   int vsize = x.Size() / 2;
-   // int vsize  = pfes->GetVSize();
-   // int tvsize = pfes->GetTrueVSize();
-
-   cout << "x.Size/2 returns " << vsize << endl;
-
-   double s = (conv_ == ComplexOperator::HERMITIAN)?1.0:-1.0;
-
-   // Allocate temporary vectors
-   Vector b_0(vsize);  b_0 = 0.0;
-   // Vector B_0(tvsize); B_0 = 0.0;
-
-   // Extract the real and imaginary parts of the input vectors
-   // MFEM_ASSERT(x.Size() == 2 * vsize, "Input GridFunction of incorrect size!");
-   Vector x_r(x.GetData(), vsize);
-   Vector x_i(&(x.GetData())[vsize], vsize);
-
-   MFEM_ASSERT(b.Size() == 2 * vsize, "Input LinearForm of incorrect size!");
-   Vector b_r(b.GetData(), vsize);
-   Vector b_i(&(b.GetData())[vsize], vsize);
-   b_i *= s;
-   /*
-   X.SetSize(2 * tvsize);
-   Vector X_r(X.GetData(), tvsize);
-   Vector X_i(&(X.GetData())[tvsize], tvsize);
-
-   B.SetSize(2 * tvsize);
-   Vector B_r(B.GetData(), tvsize);
-   Vector B_i(&(B.GetData())[tvsize], tvsize);
-   */
-   OperatorHandle A_r, A_i;
-   Vector X_0, B_0;
-   cout << "pblfr fls 1" << endl << flush;
-   b_0 = b_r;
-   pblfr_->FormLinearSystem(ess_tdof_list, x_r, b_0, A_r, X_0, B_0, ci);
-
-   int tvsize = B_0.Size();
-   X.SetSize(2 * tvsize);
-   B.SetSize(2 * tvsize);
-   Vector X_r(X.GetData(), tvsize);
-   Vector X_i(&(X.GetData())[tvsize], tvsize);
-   Vector B_r(B.GetData(), tvsize);
-   Vector B_i(&(B.GetData())[tvsize], tvsize);
-   X_r = X_0; B_r = B_0;
-   cout << "pblfi fls 1" << endl << flush;
-   b_0 = 0.0;
-   pblfi_->FormLinearSystem(ess_tdof_list, x_i, b_0, A_i, X_0, B_0, false);
-   B_r -= B_0;
-   cout << "pblfr fls 2" << endl << flush;
-   b_0 = b_i;
-   pblfr_->FormLinearSystem(ess_tdof_list, x_i, b_0, A_r, X_0, B_0, ci);
-   X_i = X_0; B_i = B_0;
-   cout << "pblfi fls 2" << endl << flush;
-   b_0 = 0.0;
-   pblfi_->FormLinearSystem(ess_tdof_list, x_r, b_0, A_i, X_0, B_0, false);
-   B_i += B_0;
-
-   B_i *= s;
-   b_i *= s;
-
-   // A = A_r + i A_i
-   A.Clear();
-   if ( A_r.Type() == Operator::Hypre_ParCSR &&
-        A_i.Type() == Operator::Hypre_ParCSR )
-   {
-      ComplexHypreParMatrix * A_hyp =
-         new ComplexHypreParMatrix(A_r.As<HypreParMatrix>(),
-                                   A_i.As<HypreParMatrix>(),
-                                   A_r.OwnsOperator(),
-                                   A_i.OwnsOperator(),
-                                   conv_);
-      A.Reset<ComplexHypreParMatrix>(A_hyp, true);
-   }
-   else
-   {
-      ComplexOperator * A_op =
-         new ComplexOperator(A_r.As<Operator>(),
-                             A_i.As<Operator>(),
-                             A_r.OwnsOperator(),
-                             A_i.OwnsOperator(),
-                             conv_);
-      A.Reset<ComplexOperator>(A_op, true);
-   }
-}
-
-void
-ParSesquilinearForm::RecoverFEMSolution(const Vector &X, const Vector &b,
-                                        Vector &x)
-{
-   ParFiniteElementSpace * pfes = pblfr_->ParFESpace();
-
-   const Operator &P = *pfes->GetProlongationMatrix();
-
-   int vsize  = pfes->GetVSize();
-   int tvsize = X.Size() / 2;
-
-   Vector X_r(X.GetData(), tvsize);
-   Vector X_i(&(X.GetData())[tvsize], tvsize);
-
-   Vector x_r(x.GetData(), vsize);
-   Vector x_i(&(x.GetData())[vsize], vsize);
-
-   // Apply conforming prolongation
-   P.Mult(X_r, x_r);
-   P.Mult(X_i, x_i);
-}
-
-void
-ParSesquilinearForm::Update(FiniteElementSpace *nfes)
-{
-   if ( pblfr_ ) { pblfr_->Update(nfes); }
-   if ( pblfi_ ) { pblfi_->Update(nfes); }
-}
-
-
-#endif // MFEM_USE_MPI
-
-}

fem/complex_fem.hpp

@@ -1,356 +0,0 @@
-// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
-// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
-// reserved. See file COPYRIGHT for details.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability see http://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the GNU Lesser General Public License (as published by the Free
-// Software Foundation) version 2.1 dated February 1999.
-
-#ifndef MFEM_COMPLEX_FEM
-#define MFEM_COMPLEX_FEM
-
-#include "../linalg/complex_operator.hpp"
-#include "gridfunc.hpp"
-#include "linearform.hpp"
-#include "bilinearform.hpp"
-#ifdef MFEM_USE_MPI
-#include "pgridfunc.hpp"
-#include "plinearform.hpp"
-#include "pbilinearform.hpp"
-#endif
-#include <complex>
-
-namespace mfem
-{
-
-/// Class for complex-valued grid function - Vector with associated FE space.
-class ComplexGridFunction : public Vector
-{
-private:
-
-   GridFunction * gfr_;
-   GridFunction * gfi_;
-
-protected:
-   void Destroy() { delete gfr_; delete gfi_; }
-
-public:
-
-   /* @brief Construct a ComplexGridFunction associated with the
-      FiniteElementSpace @a *f. */
-   ComplexGridFunction(FiniteElementSpace *f);
-
-   void Update();
-
-   /// Assign constant values to the ComplexGridFunction data.
-   ComplexGridFunction &operator=(const std::complex<double> & value)
-   { *gfr_ = value.real(); *gfi_ = value.imag(); return *this; }
-
-   virtual void ProjectCoefficient(Coefficient &real_coeff,
-                                   Coefficient &imag_coeff);
-   virtual void ProjectCoefficient(VectorCoefficient &real_vcoeff,
-                                   VectorCoefficient &imag_vcoeff);
-
-   virtual void ProjectBdrCoefficient(Coefficient &real_coeff,
-                                      Coefficient &imag_coeff,
-                                      Array<int> &attr);
-   virtual void ProjectBdrCoefficientNormal(VectorCoefficient &real_coeff,
-                                            VectorCoefficient &imag_coeff,
-                                            Array<int> &attr);
-   virtual void ProjectBdrCoefficientTangent(VectorCoefficient &real_coeff,
-                                             VectorCoefficient &imag_coeff,
-                                             Array<int> &attr);
-
-   FiniteElementSpace *FESpace() { return gfr_->FESpace(); }
-   const FiniteElementSpace *FESpace() const { return gfr_->FESpace(); }
-
-   GridFunction & real() { return *gfr_; }
-   GridFunction & imag() { return *gfi_; }
-   const GridFunction & real() const { return *gfr_; }
-   const GridFunction & imag() const { return *gfi_; }
-
-   /// Destroys grid function.
-   virtual ~ComplexGridFunction() { Destroy(); }
-
-};
-
-class ComplexLinearForm : public Vector
-{
-private:
-   ComplexOperator::Convention conv_;
-
-protected:
-   LinearForm * lfr_;
-   LinearForm * lfi_;
-
-   // HYPRE_Int * tdof_offsets_;
-
-public:
-
-   ComplexLinearForm(FiniteElementSpace *fes,
-                     ComplexOperator::Convention
-                     convention = ComplexOperator::HERMITIAN);
-
-   virtual ~ComplexLinearForm();
-
-   /// Adds new Domain Integrator.
-   void AddDomainIntegrator(LinearFormIntegrator *lfi_real,
-                            LinearFormIntegrator *lfi_imag);
-
-   FiniteElementSpace *FESpace() const { return lfr_->FESpace(); }
-
-   LinearForm & real() { return *lfr_; }
-   LinearForm & imag() { return *lfi_; }
-   const LinearForm & real() const { return *lfr_; }
-   const LinearForm & imag() const { return *lfi_; }
-
-   void Update();
-   void Update(FiniteElementSpace *f);
-
-   /// Assembles the linear form i.e. sums over all domain/bdr integrators.
-   void Assemble();
-
-   std::complex<double> operator()(const ComplexGridFunction &gf) const;
-
-};
-
-// Class for sesquilinear form
-class SesquilinearForm
-{
-private:
-   ComplexOperator::Convention conv_;
-
-   //protected:
-   BilinearForm *blfr_;
-   BilinearForm *blfi_;
-
-public:
-   SesquilinearForm(FiniteElementSpace *fes,
-                    ComplexOperator::Convention
-                    convention = ComplexOperator::HERMITIAN);
-
-   ComplexOperator::Convention GetConvention() const { return conv_; }
-   void SetConvention(const ComplexOperator::Convention &
-                      convention) { conv_  = convention; }
-
-   BilinearForm & real() { return *blfr_; }
-   BilinearForm & imag() { return *blfi_; }
-   const BilinearForm & real() const { return *blfr_; }
-   const BilinearForm & imag() const { return *blfi_; }
-
-   /// Adds new Domain Integrator.
-   void AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
-                            BilinearFormIntegrator *bfi_imag);
-
-   /// Adds new Boundary Integrator.
-   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
-                              BilinearFormIntegrator *bfi_imag);
-
-   /// Adds new Boundary Integrator, restricted to specific boundary attributes.
-   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
-                              BilinearFormIntegrator *bfi_imag,
-                              Array<int> &bdr_marker);
-
-   /// Assemble the local matrix
-   void Assemble(int skip_zeros = 1);
-
-   /// Finalizes the matrix initialization.
-   void Finalize(int skip_zeros = 1);
-
-   /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
-   /** The returned matrix has to be deleted by the caller. */
-   ComplexSparseMatrix *AssembleCompSpMat();
-
-   /// Return the parallel FE space associated with the ParBilinearForm.
-   FiniteElementSpace *FESpace() const { return blfr_->FESpace(); }
-
-   void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
-                         OperatorHandle &A, Vector &X, Vector &B,
-                         int copy_interior = 0);
-
-   /** Call this method after solving a linear system constructed using the
-       FormLinearSystem method to recover the solution as a ParGridFunction-size
-       vector in x. Use the same arguments as in the FormLinearSystem call. */
-   virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);
-
-   virtual void Update(FiniteElementSpace *nfes = NULL);
-
-   virtual ~SesquilinearForm();
-};
-
-#ifdef MFEM_USE_MPI
-
-/// Class for complex-valued grid function - Vector with associated FE space.
-class ParComplexGridFunction : public Vector
-{
-private:
-
-   ParGridFunction * pgfr_;
-   ParGridFunction * pgfi_;
-
-protected:
-   void Destroy() { delete pgfr_; delete pgfi_; }
-
-public:
-
-   /* @brief Construct a ParComplexGridFunction associated with the
-      ParFiniteElementSpace @a *f. */
-   ParComplexGridFunction(ParFiniteElementSpace *pf);
-
-   void Update();
-
-   /// Assign constant values to the ParComplexGridFunction data.
-   ParComplexGridFunction &operator=(const std::complex<double> & value)
-   { *pgfr_ = value.real(); *pgfi_ = value.imag(); return *this; }
-
-   virtual void ProjectCoefficient(Coefficient &real_coeff,
-                                   Coefficient &imag_coeff);
-   virtual void ProjectCoefficient(VectorCoefficient &real_vcoeff,
-                                   VectorCoefficient &imag_vcoeff);
-
-   virtual void ProjectBdrCoefficient(Coefficient &real_coeff,
-                                      Coefficient &imag_coeff,
-                                      Array<int> &attr);
-   virtual void ProjectBdrCoefficientNormal(VectorCoefficient &real_coeff,
-                                            VectorCoefficient &imag_coeff,
-                                            Array<int> &attr);
-   virtual void ProjectBdrCoefficientTangent(VectorCoefficient &real_coeff,
-                                             VectorCoefficient &imag_coeff,
-                                             Array<int> &attr);
-
-   void Distribute(const Vector *tv);
-   void Distribute(const Vector &tv) { Distribute(&tv); }
-
-   /// Returns the vector restricted to the true dofs.
-   void ParallelProject(Vector &tv) const;
-
-   FiniteElementSpace *FESpace() { return pgfr_->FESpace(); }
-   const FiniteElementSpace *FESpace() const { return pgfr_->FESpace(); }
-
-   ParGridFunction & real() { return *pgfr_; }
-   ParGridFunction & imag() { return *pgfi_; }
-   const ParGridFunction & real() const { return *pgfr_; }
-   const ParGridFunction & imag() const { return *pgfi_; }
-
-   /// Destroys grid function.
-   virtual ~ParComplexGridFunction() { Destroy(); }
-
-};
-
-class ParComplexLinearForm : public Vector
-{
-private:
-   ComplexOperator::Convention conv_;
-
-protected:
-   ParLinearForm * plfr_;
-   ParLinearForm * plfi_;
-
-   HYPRE_Int * tdof_offsets_;
-
-public:
-
-   ParComplexLinearForm(ParFiniteElementSpace *pf,
-                        ComplexOperator::Convention
-                        convention = ComplexOperator::HERMITIAN);
-
-   virtual ~ParComplexLinearForm();
-
-   /// Adds new Domain Integrator.
-   void AddDomainIntegrator(LinearFormIntegrator *lfi_real,
-                            LinearFormIntegrator *lfi_imag);
-
-   ParFiniteElementSpace *ParFESpace() const { return plfr_->ParFESpace(); }
-
-   ParLinearForm & real() { return *plfr_; }
-   ParLinearForm & imag() { return *plfi_; }
-   const ParLinearForm & real() const { return *plfr_; }
-   const ParLinearForm & imag() const { return *plfi_; }
-
-   void Update(ParFiniteElementSpace *pf = NULL);
-
-   /// Assembles the linear form i.e. sums over all domain/bdr integrators.
-   void Assemble();
-
-   /// Assemble the vector on the true dofs, i.e. P^t v.
-   void ParallelAssemble(Vector &tv);
-
-   /// Returns the vector assembled on the true dofs, i.e. P^t v.
-   HypreParVector *ParallelAssemble();
-
-   std::complex<double> operator()(const ParComplexGridFunction &gf) const;
-
-};
-
-// Class for parallel sesquilinear form
-class ParSesquilinearForm
-{
-private:
-   ComplexOperator::Convention conv_;
-
-   //protected:
-   ParBilinearForm *pblfr_;
-   ParBilinearForm *pblfi_;
-
-public:
-   ParSesquilinearForm(ParFiniteElementSpace *pf,
-                       ComplexOperator::Convention
-                       convention = ComplexOperator::HERMITIAN);
-
-   ComplexOperator::Convention GetConvention() const { return conv_; }
-   void SetConvention(const ComplexOperator::Convention &
-                      convention) { conv_  = convention; }
-
-   ParBilinearForm & real() { return *pblfr_; }
-   ParBilinearForm & imag() { return *pblfi_; }
-   const ParBilinearForm & real() const { return *pblfr_; }
-   const ParBilinearForm & imag() const { return *pblfi_; }
-
-   /// Adds new Domain Integrator.
-   void AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
-                            BilinearFormIntegrator *bfi_imag);
-
-   /// Adds new Boundary Integrator.
-   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
-                              BilinearFormIntegrator *bfi_imag);
-
-   /// Adds new Boundary Integrator, restricted to specific boundary attributes.
-   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
-                              BilinearFormIntegrator *bfi_imag,
-                              Array<int> &bdr_marker);
-
-   /// Assemble the local matrix
-   void Assemble(int skip_zeros = 1);
-
-   /// Finalizes the matrix initialization.
-   void Finalize(int skip_zeros = 1);
-
-   /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
-   /** The returned matrix has to be deleted by the caller. */
-   ComplexHypreParMatrix *ParallelAssemble();
-
-   /// Return the parallel FE space associated with the ParBilinearForm.
-   ParFiniteElementSpace *ParFESpace() const { return pblfr_->ParFESpace(); }
-
-   void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
-                         OperatorHandle &A, Vector &X, Vector &B,
-                         int copy_interior = 0);
-
-   /** Call this method after solving a linear system constructed using the
-       FormLinearSystem method to recover the solution as a ParGridFunction-size
-       vector in x. Use the same arguments as in the FormLinearSystem call. */
-   virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);
-
-   virtual void Update(FiniteElementSpace *nfes = NULL);
-
-   virtual ~ParSesquilinearForm();
-};
-
-#endif // MFEM_USE_MPI
-
-}
-
-#endif // MFEM_COMPLEX_FEM

fem/fem.hpp

@@ -18,7 +18,6 @@
 #include "fe_coll.hpp"
 #include "eltrans.hpp"
 #include "coefficient.hpp"
-#include "complex_fem.hpp"
 #include "lininteg.hpp"
 #include "nonlininteg.hpp"
 #include "bilininteg.hpp"

fem/linearform.hpp

@@ -81,15 +81,6 @@ public:
        Update(FiniteElementSpace *, Vector &, int). */
    LinearForm() { fes = NULL; extern_lfs = 0; }
 
-   /// Construct a LinearForm using previously allocated array @a data.
-   /** The LinearForm does not assume ownership of @a data which is assumed to
-       be of size at least `f->GetVSize()`. Similar to the Vector constructor
-       for externally allocated array, the pointer @a data can be NULL. The data
-       array can be replaced later using the method SetData().
-    */
-   LinearForm(FiniteElementSpace *f, double *data) : Vector(data, f->GetVSize())
-   { fes = f; }
-
    /// Copy assignment. Only the data of the base class Vector is copied.
    /** It is assumed that this object and @a rhs use FiniteElementSpace%s that
        have the same size.

fem/plinearform.hpp

@@ -45,16 +45,6 @@ public:
    /** The pointer @a pf is not owned by the newly constructed object. */
    ParLinearForm(ParFiniteElementSpace *pf) : LinearForm(pf) { pfes = pf; }
 
-   /// Construct a ParLinearForm using previously allocated array @a data.
-   /** The ParLinearForm does not assume ownership of @a data which is assumed
-       to be of size at least `pf->GetVSize()`. Similar to the LinearForm and
-       Vector constructors for externally allocated array, the pointer @a data
-       can be NULL. The data array can be replaced later using the method
-       SetData().
-    */
-   ParLinearForm(ParFiniteElementSpace *pf, double *data) :
-      LinearForm(pf, data), pfes(pf) { }
-
    /** @brief Create a ParLinearForm on the ParFiniteElementSpace @a *pf, using
        the same integrators as the ParLinearForm @a *plf.
 

general/cuda.cpp

@@ -10,17 +10,27 @@
 // Software Foundation) version 2.1 dated February 1999.
 
 #include "cuda.hpp"
+#include "globals.hpp"
 
 namespace mfem
 {
 
+#ifdef MFEM_USE_CUDA
+void mfem_cuda_error(cudaError_t err, const char *expr, const char *func,
+                     const char *file, int line)
+{
+   mfem::err << "CUDA error: (" << expr << ") failed with error:\n --> "
+             << cudaGetErrorString(err)
+             << "\n ... in function: " << func
+             << "\n ... in file: " << file << ':' << line << '\n';
+   mfem_error();
+}
+#endif
+
 void* CuMemAlloc(void** dptr, size_t bytes)
 {
 #ifdef MFEM_USE_CUDA
-   if (CUDA_SUCCESS != ::cuMemAlloc((CUdeviceptr*)dptr, bytes))
-   {
-      mfem_error("Error in CuMemAlloc");
-   }
+   MFEM_CUDA_CHECK(cudaMalloc(dptr, bytes));
 #endif
    return *dptr;
 }
@@ -28,10 +38,7 @@ void* CuMemAlloc(void** dptr, size_t bytes)
 void* CuMemFree(void *dptr)
 {
 #ifdef MFEM_USE_CUDA
-   if (CUDA_SUCCESS != ::cuMemFree((CUdeviceptr)dptr))
-   {
-      mfem_error("Error in CuMemFree");
-   }
+   MFEM_CUDA_CHECK(cudaFree(dptr));
 #endif
    return dptr;
 }
@@ -39,22 +46,15 @@ void* CuMemFree(void *dptr)
 void* CuMemcpyHtoD(void* dst, const void* src, size_t bytes)
 {
 #ifdef MFEM_USE_CUDA
-   if (CUDA_SUCCESS != ::cuMemcpyHtoD((CUdeviceptr)dst, src, bytes))
-   {
-      mfem_error("Error in CuMemcpyHtoD");
-   }
+   MFEM_CUDA_CHECK(cudaMemcpy(dst, src, bytes, cudaMemcpyHostToDevice));
 #endif
    return dst;
 }
 
-void* CuMemcpyHtoDAsync(void* dst, const void* src, size_t bytes, void *s)
+void* CuMemcpyHtoDAsync(void* dst, const void* src, size_t bytes)
 {
 #ifdef MFEM_USE_CUDA
-   if (CUDA_SUCCESS !=
-       ::cuMemcpyHtoDAsync((CUdeviceptr)dst, src, bytes, (CUstream)s))
-   {
-      mfem_error("Error in CuMemcpyHtoDAsync");
-   }
+   MFEM_CUDA_CHECK(cudaMemcpyAsync(dst, src, bytes, cudaMemcpyHostToDevice));
 #endif
    return dst;
 }
@@ -62,24 +62,15 @@ void* CuMemcpyHtoDAsync(void* dst, const void* src, size_t bytes, void *s)
 void* CuMemcpyDtoD(void* dst, void* src, size_t bytes)
 {
 #ifdef MFEM_USE_CUDA
-   if (CUDA_SUCCESS !=
-       ::cuMemcpyDtoD((CUdeviceptr)dst, (CUdeviceptr)src, bytes))
-   {
-      mfem_error("Error in CuMemcpyDtoD");
-   }
+   MFEM_CUDA_CHECK(cudaMemcpy(dst, src, bytes, cudaMemcpyDeviceToDevice));
 #endif
    return dst;
 }
 
-void* CuMemcpyDtoDAsync(void* dst, void* src, size_t bytes, void *s)
+void* CuMemcpyDtoDAsync(void* dst, void* src, size_t bytes)
 {
 #ifdef MFEM_USE_CUDA
-   if (CUDA_SUCCESS !=
-       ::cuMemcpyDtoDAsync((CUdeviceptr)dst, (CUdeviceptr)src,
-                           bytes, (CUstream)s))
-   {
-      mfem_error("Error in CuMemcpyDtoDAsync");
-   }
+   MFEM_CUDA_CHECK(cudaMemcpyAsync(dst, src, bytes, cudaMemcpyDeviceToDevice));
 #endif
    return dst;
 }
@@ -87,10 +78,7 @@ void* CuMemcpyDtoDAsync(void* dst, void* src, size_t bytes, void *s)
 void* CuMemcpyDtoH(void *dst, void *src, size_t bytes)
 {
 #ifdef MFEM_USE_CUDA
-   if (CUDA_SUCCESS != ::cuMemcpyDtoH(dst, (CUdeviceptr)src, bytes))
-   {
-      mfem_error("Error in CuMemcpyDtoH");
-   }
+   MFEM_CUDA_CHECK(cudaMemcpy(dst, src, bytes, cudaMemcpyDeviceToHost));
 #endif
    return dst;
 }
@@ -98,11 +86,7 @@ void* CuMemcpyDtoH(void *dst, void *src, size_t bytes)
 void* CuMemcpyDtoHAsync(void* dst, void* src, size_t bytes, void *s)
 {
 #ifdef MFEM_USE_CUDA
-   if (CUDA_SUCCESS !=
-       ::cuMemcpyDtoHAsync(dst, (CUdeviceptr)src, bytes, (CUstream)s))
-   {
-      mfem_error("Error in CuMemcpyDtoHAsync");
-   }
+   MFEM_CUDA_CHECK(cudaMemcpyAsync(dst, src, bytes, cudaMemcpyDeviceToHost));
 #endif
    return dst;
 }

general/cuda.hpp

@@ -26,89 +26,33 @@
 #ifdef MFEM_USE_CUDA
 #define MFEM_ATTR_DEVICE __device__
 #define MFEM_ATTR_HOST_DEVICE __host__ __device__
-// Define the CUDA debug macros:
-// - MFEM_CUDA_CHECK_DRV(x) where 'x' returns/is type 'CUresult'
-// - MFEM_CUDA_CHECK_RT(x)  where 'x' returns/is type 'cudaError_t'
-#ifdef MFEM_DEBUG
-#define MFEM_CUDA_CHECK_DRV(x) \
-   do \
-   { \
-      CUresult err = (x); \
-      if (err != CUDA_SUCCESS) \
-      { \
-         const char *error_string; \
-         cuGetErrorString(err, &error_string); \
-         _MFEM_MESSAGE("CUDA error: (" << #x \
-                       << ") failed with error:\n --> " \
-                       << error_string, 0); \
-      } \
-   } \
-   while (0)
-#define MFEM_CUDA_CHECK_RT(x) \
+// Define a CUDA error check macro, MFEM_CUDA_CHECK(x), where x returns/is of
+// type 'cudaError_t'. This macro evaluates 'x' and raises an error if the
+// result is not cudaSuccess.
+#define MFEM_CUDA_CHECK(x) \
    do \
    { \
       cudaError_t err = (x); \
       if (err != cudaSuccess) \
       { \
-         _MFEM_MESSAGE("CUDA error: (" << #x \
-                       << ") failed with error:\n --> " \
-                       << cudaGetErrorString(err), 0); \
+         mfem_cuda_error(err, #x, _MFEM_FUNC_NAME, __FILE__, __LINE__); \
       } \
    } \
    while (0)
-#else
-#define MFEM_CUDA_CHECK_DRV(x) x
-#define MFEM_CUDA_CHECK_RT(x) x
-#endif
 #else // MFEM_USE_CUDA
 #define MFEM_ATTR_DEVICE
 #define MFEM_ATTR_HOST_DEVICE
-typedef int CUdevice;
-typedef int CUcontext;
-typedef void* CUstream;
 #endif // MFEM_USE_CUDA
 
 
 namespace mfem
 {
 
-// Define 'atomicAdd' function.
-#ifdef __CUDA_ARCH__
-#if __CUDA_ARCH__ < 600
-static __device__ inline double atomicAdd(double* address, double val)
-{
-   unsigned long long int* address_as_ull = (unsigned long long int*)address;
-   unsigned long long int old = *address_as_ull, assumed;
-   do
-   {
-      assumed = old;
-      old =
-         atomicCAS(address_as_ull, assumed,
-                   __double_as_longlong(val +
-                                        __longlong_as_double(assumed)));
-      // Note: uses integer comparison to avoid hang in case of NaN
-      // (since NaN != NaN)
-   }
-   while (assumed != old);
-   return __longlong_as_double(old);
-}
-#endif // __CUDA_ARCH__ < 600
-template<typename T> MFEM_ATTR_DEVICE
-inline T AtomicAdd(T volatile *address, T val)
-{
-   return atomicAdd((T *)address, val);
-}
-#else // __CUDA_ARCH__
-template<typename T> inline T AtomicAdd(T volatile *address, T val)
-{
-#ifdef MFEM_USE_OPENMP
-   #pragma omp atomic
+#ifdef MFEM_USE_CUDA
+// Function used by the macro MFEM_CUDA_CHECK.
+void mfem_cuda_error(cudaError_t err, const char *expr, const char *func,
+                     const char *file, int line);
 #endif
-   *address += val;
-   return *address;
-}
-#endif // __CUDA_ARCH__
-
 
 /// Allocates device memory
 void* CuMemAlloc(void **d_ptr, size_t bytes);
@@ -120,20 +64,19 @@ void* CuMemFree(void *d_ptr);
 void* CuMemcpyHtoD(void *d_dst, const void *h_src, size_t bytes);
 
 /// Copies memory from Host to Device
-void* CuMemcpyHtoDAsync(void *d_dst, const void *h_src,
-                        size_t bytes, void *stream);
+void* CuMemcpyHtoDAsync(void *d_dst, const void *h_src, size_t bytes);
 
 /// Copies memory from Device to Device
 void* CuMemcpyDtoD(void *d_dst, void *d_src, size_t bytes);
 
 /// Copies memory from Device to Device
-void* CuMemcpyDtoDAsync(void *d_dst, void *d_src, size_t bytes, void *stream);
+void* CuMemcpyDtoDAsync(void *d_dst, void *d_src, size_t bytes);
 
 /// Copies memory from Device to Host
 void* CuMemcpyDtoH(void *h_dst, void *d_src, size_t bytes);
 
 /// Copies memory from Device to Host
-void* CuMemcpyDtoHAsync(void *h_dst, void *d_src, size_t bytes, void *stream);
+void* CuMemcpyDtoHAsync(void *h_dst, void *d_src, size_t bytes);
 
 } // namespace mfem
 

general/device.cpp

@@ -24,9 +24,6 @@ namespace mfem
 namespace internal
 {
 
-CUstream *cuStream = NULL;
-static CUdevice cuDevice;
-static CUcontext cuContext;
 OccaDevice occaDevice;
 
 // Backends listed by priority, high to low:
@@ -100,14 +97,9 @@ void Device::Print(std::ostream &out)
 #ifdef MFEM_USE_CUDA
 static void DeviceSetup(const int dev, int &ngpu)
 {
-   cudaGetDeviceCount(&ngpu);
-   MFEM_VERIFY(ngpu>0, "No CUDA device found!");
-   cuInit(0);
-   cuDeviceGet(&internal::cuDevice, dev);
-   cuCtxCreate(&internal::cuContext, CU_CTX_SCHED_AUTO, internal::cuDevice);
-   internal::cuStream = new CUstream;
-   MFEM_VERIFY(internal::cuStream, "CUDA stream could not be created!");
-   cuStreamCreate(internal::cuStream, CU_STREAM_DEFAULT);
+   MFEM_CUDA_CHECK(cudaGetDeviceCount(&ngpu));
+   MFEM_VERIFY(ngpu > 0, "No CUDA device found!");
+   MFEM_CUDA_CHECK(cudaSetDevice(dev));
 }
 #endif
 
@@ -125,7 +117,7 @@ static void RajaDeviceSetup(const int dev, int &ngpu)
 #endif
 }
 
-static void OccaDeviceSetup(CUdevice cu_dev, CUcontext cu_ctx)
+static void OccaDeviceSetup(const int dev)
 {
 #ifdef MFEM_USE_OCCA
    const int cpu  = Device::Allows(Backend::OCCA_CPU);
@@ -138,7 +130,8 @@ static void OccaDeviceSetup(CUdevice cu_dev, CUcontext cu_ctx)
    if (cuda)
    {
 #if OCCA_CUDA_ENABLED
-      internal::occaDevice = occa::cuda::wrapDevice(cu_dev, cu_ctx);
+      std::string mode("mode: 'CUDA', device_id : ");
+      internal::occaDevice.setup(mode.append(1,'0'+dev));
 #else
       MFEM_ABORT("the OCCA CUDA backend requires OCCA built with CUDA!");
 #endif
@@ -197,22 +190,10 @@ void Device::Setup(const int device)
                "the OpenMP and RAJA OpenMP backends require MFEM built with"
                " MFEM_USE_OPENMP=YES");
 #endif
-   // The check for MFEM_USE_OCCA is in the function OccaDeviceSetup().
-
-   // We initialize CUDA and/or RAJA_CUDA first so OccaDeviceSetup() can reuse
-   // the same initialized cuDevice and cuContext objects when OCCA_CUDA is
-   // enabled.
    if (Allows(Backend::CUDA)) { CudaDeviceSetup(dev, ngpu); }
    if (Allows(Backend::RAJA_CUDA)) { RajaDeviceSetup(dev, ngpu); }
-   if (Allows(Backend::OCCA_MASK))
-   {
-      OccaDeviceSetup(internal::cuDevice, internal::cuContext);
-   }
-}
-
-Device::~Device()
-{
-   delete internal::cuStream;
+   // The check for MFEM_USE_OCCA is in the function OccaDeviceSetup().
+   if (Allows(Backend::OCCA_MASK)) { OccaDeviceSetup(dev); }
 }
 
 } // mfem

general/device.hpp

@@ -181,8 +181,6 @@ public:
        Backend::*_MASK, or combinations of those. */
    static inline bool Allows(unsigned long b_mask)
    { return Get().allowed_backends & b_mask; }
-
-   ~Device();
 };
 
 } // mfem

general/forall.hpp

@@ -22,6 +22,9 @@
 
 #ifdef MFEM_USE_RAJA
 #include "RAJA/RAJA.hpp"
+#if defined(RAJA_ENABLE_CUDA) && !defined(MFEM_USE_CUDA)
+#error When RAJA is built with CUDA, MFEM_USE_CUDA=YES is required
+#endif
 #endif
 
 namespace mfem
@@ -106,8 +109,7 @@ void CuWrap(const int N, DBODY &&d_body)
    if (N==0) { return; }
    const int GRID = (N+BLOCKS-1)/BLOCKS;
    CuKernel<<<GRID,BLOCKS>>>(N,d_body);
-   const cudaError_t last = cudaGetLastError();
-   MFEM_VERIFY(last == cudaSuccess, cudaGetErrorString(last));
+   MFEM_CUDA_CHECK(cudaGetLastError());
 }
 
 #else  // MFEM_USE_CUDA

general/mem_manager.cpp

@@ -312,7 +312,6 @@ void MemoryManager::Pull(const void *ptr, const std::size_t bytes)
    { mfem_error("Unknown pointer to pull from!"); }
 }
 
-namespace internal { extern CUstream *cuStream; }
 void* MemoryManager::Memcpy(void *dst, const void *src,
                             const std::size_t bytes, const bool async)
 {
@@ -322,7 +321,7 @@ void* MemoryManager::Memcpy(void *dst, const void *src,
    const bool run_on_host = !Device::Allows(Backend::DEVICE_MASK);
    if (run_on_host) { return std::memcpy(dst, src, bytes); }
    if (!async) { return CuMemcpyDtoD(d_dst, d_src, bytes); }
-   return CuMemcpyDtoDAsync(d_dst, d_src, bytes, internal::cuStream);
+   return CuMemcpyDtoDAsync(d_dst, d_src, bytes);
 }
 
 void MemoryManager::RegisterCheck(void *ptr)

general/occa.hpp

@@ -13,7 +13,6 @@
 #define MFEM_OCCA_HPP
 
 #include "../config/config.hpp"
-#include "cuda.hpp" // for CUdevice, CUcontext
 
 #ifdef MFEM_USE_OCCA
 #include <occa.hpp>

linalg/complex_operator.cpp

@@ -10,8 +10,6 @@
 // Software Foundation) version 2.1 dated February 1999.
 
 #include "complex_operator.hpp"
-#include <set>
-#include <map>
 
 namespace mfem
 {
@@ -41,30 +39,6 @@ ComplexOperator::~ComplexOperator()
    delete v_;
 }
 
-Operator & ComplexOperator::real()
-{
-   MFEM_ASSERT(Op_Real_, "ComplexOperator has no real part!");
-   return *Op_Real_;
-}
-
-Operator & ComplexOperator::imag()
-{
-   MFEM_ASSERT(Op_Imag_, "ComplexOperator has no imaginary part!");
-   return *Op_Imag_;
-}
-
-const Operator & ComplexOperator::real() const
-{
-   MFEM_ASSERT(Op_Real_, "ComplexOperator has no real part!");
-   return *Op_Real_;
-}
-
-const Operator & ComplexOperator::imag() const
-{
-   MFEM_ASSERT(Op_Imag_, "ComplexOperator has no imaginary part!");
-   return *Op_Imag_;
-}
-
 void ComplexOperator::Mult(const Vector &x, Vector &y) const
 {
    double * x_data = x.GetData();
@@ -146,30 +120,6 @@ void ComplexOperator::MultTranspose(const Vector &x_r, const Vector &x_i,
 }
 
 
-SparseMatrix & ComplexSparseMatrix::real()
-{
-   MFEM_ASSERT(Op_Real_, "ComplexSparseMatrix has no real part!");
-   return dynamic_cast<SparseMatrix &>(*Op_Real_);
-}
-
-SparseMatrix & ComplexSparseMatrix::imag()
-{
-   MFEM_ASSERT(Op_Imag_, "ComplexSparseMatrix has no imaginary part!");
-   return dynamic_cast<SparseMatrix &>(*Op_Imag_);
-}
-
-const SparseMatrix & ComplexSparseMatrix::real() const
-{
-   MFEM_ASSERT(Op_Real_, "ComplexSparseMatrix has no real part!");
-   return dynamic_cast<const SparseMatrix &>(*Op_Real_);
-}
-
-const SparseMatrix & ComplexSparseMatrix::imag() const
-{
-   MFEM_ASSERT(Op_Imag_, "ComplexSparseMatrix has no imaginary part!");
-   return dynamic_cast<const SparseMatrix &>(*Op_Imag_);
-}
-
 SparseMatrix * ComplexSparseMatrix::GetSystemMatrix() const
 {
    SparseMatrix * A_r = dynamic_cast<SparseMatrix*>(Op_Real_);
@@ -234,358 +184,4 @@ SparseMatrix * ComplexSparseMatrix::GetSystemMatrix() const
    return new SparseMatrix(I, J, D, this->Height(), this->Width());
 }
 
-#ifdef MFEM_USE_MPI
-
-ComplexHypreParMatrix::ComplexHypreParMatrix(HypreParMatrix * A_Real,
-                                             HypreParMatrix * A_Imag,
-                                             bool ownReal, bool ownImag,
-                                             Convention convention)
-   : ComplexOperator(A_Real, A_Imag, ownReal, ownImag, convention)
-{
-   comm_ = (A_Real) ? A_Real->GetComm() :
-           ((A_Imag) ? A_Imag->GetComm() : MPI_COMM_WORLD);
-
-   MPI_Comm_rank(comm_, &myid_);
-   MPI_Comm_size(comm_, &nranks_);
-}
-
-HypreParMatrix & ComplexHypreParMatrix::real()
-{
-   MFEM_ASSERT(Op_Real_, "ComplexHypreParMatrix has no real part!");
-   return dynamic_cast<HypreParMatrix &>(*Op_Real_);
-}
-
-HypreParMatrix & ComplexHypreParMatrix::imag()
-{
-   MFEM_ASSERT(Op_Imag_, "ComplexHypreParMatrix has no imaginary part!");
-   return dynamic_cast<HypreParMatrix &>(*Op_Imag_);
-}
-
-const HypreParMatrix & ComplexHypreParMatrix::real() const
-{
-   MFEM_ASSERT(Op_Real_, "ComplexHypreParMatrix has no real part!");
-   return dynamic_cast<const HypreParMatrix &>(*Op_Real_);
-}
-
-const HypreParMatrix & ComplexHypreParMatrix::imag() const
-{
-   MFEM_ASSERT(Op_Imag_, "ComplexHypreParMatrix has no imaginary part!");
-   return dynamic_cast<const HypreParMatrix &>(*Op_Imag_);
-}
-
-HypreParMatrix * ComplexHypreParMatrix::GetSystemMatrix() const
-{
-   HypreParMatrix * A_r = dynamic_cast<HypreParMatrix*>(Op_Real_);
-   HypreParMatrix * A_i = dynamic_cast<HypreParMatrix*>(Op_Imag_);
-
-   if ( A_r == NULL && A_i == NULL ) { return NULL; }
-
-   HYPRE_Int global_num_rows_r = (A_r) ? A_r->GetGlobalNumRows() : 0;
-   HYPRE_Int global_num_rows_i = (A_i) ? A_i->GetGlobalNumRows() : 0;
-   HYPRE_Int global_num_rows = std::max(global_num_rows_r, global_num_rows_i);
-
-   HYPRE_Int global_num_cols_r = (A_r) ? A_r->GetGlobalNumCols() : 0;
-   HYPRE_Int global_num_cols_i = (A_i) ? A_i->GetGlobalNumCols() : 0;
-   HYPRE_Int global_num_cols = std::max(global_num_cols_r, global_num_cols_i);
-
-   int row_starts_size = (HYPRE_AssumedPartitionCheck()) ? 2 : nranks_ + 1;
-   HYPRE_Int * row_starts = hypre_CTAlloc(HYPRE_Int, row_starts_size);
-   HYPRE_Int * col_starts = hypre_CTAlloc(HYPRE_Int, row_starts_size);
-
-   const HYPRE_Int * row_starts_z = (A_r) ? A_r->RowPart() :
-                                    ((A_i) ? A_i->RowPart() : NULL);
-   const HYPRE_Int * col_starts_z = (A_r) ? A_r->ColPart() :
-                                    ((A_i) ? A_i->ColPart() : NULL);
-
-   for (int i = 0; i < row_starts_size; i++)
-   {
-      row_starts[i] = 2 * row_starts_z[i];
-      col_starts[i] = 2 * col_starts_z[i];
-   }
-
-   SparseMatrix diag_r, diag_i, offd_r, offd_i;
-   HYPRE_Int * cmap_r, * cmap_i;
-
-   int nrows_r = 0, nrows_i = 0, ncols_r = 0, ncols_i = 0;
-   int ncols_offd_r = 0, ncols_offd_i = 0;
-   if (A_r)
-   {
-      A_r->GetDiag(diag_r);
-      A_r->GetOffd(offd_r, cmap_r);
-      nrows_r = diag_r.Height();
-      ncols_r = diag_r.Width();
-      ncols_offd_r = offd_r.Width();
-   }
-   if (A_i)
-   {
-      A_i->GetDiag(diag_i);
-      A_i->GetOffd(offd_i, cmap_i);
-      nrows_i = diag_i.Height();
-      ncols_i = diag_i.Width();
-      ncols_offd_i = offd_i.Width();
-   }
-   int nrows = std::max(nrows_r, nrows_i);
-   int ncols = std::max(ncols_r, ncols_i);
-
-   // Determine the unique set of off-diagonal columns global indices
-   std::set<int> cset;
-   for (int i=0; i<ncols_offd_r; i++)
-   {
-      cset.insert(cmap_r[i]);
-   }
-   for (int i=0; i<ncols_offd_i; i++)
-   {
-      cset.insert(cmap_i[i]);
-   }
-   int num_cols_offd = (int)cset.size();
-
-   // Exatract pointers to the various CSR arrays of the diagonal blocks
-   const int * diag_r_I = (A_r) ? diag_r.GetI() : NULL;
-   const int * diag_i_I = (A_i) ? diag_i.GetI() : NULL;
-
-   const int * diag_r_J = (A_r) ? diag_r.GetJ() : NULL;
-   const int * diag_i_J = (A_i) ? diag_i.GetJ() : NULL;
-
-   const double * diag_r_D = (A_r) ? diag_r.GetData() : NULL;
-   const double * diag_i_D = (A_i) ? diag_i.GetData() : NULL;
-
-   int diag_r_nnz = (diag_r_I) ? diag_r_I[nrows] : 0;
-   int diag_i_nnz = (diag_i_I) ? diag_i_I[nrows] : 0;
-   int diag_nnz = 2 * (diag_r_nnz + diag_i_nnz);
-
-   // Exatract pointers to the various CSR arrays of the off-diagonal blocks
-   const int * offd_r_I = (A_r) ? offd_r.GetI() : NULL;
-   const int * offd_i_I = (A_i) ? offd_i.GetI() : NULL;
-
-   const int * offd_r_J = (A_r) ? offd_r.GetJ() : NULL;
-   const int * offd_i_J = (A_i) ? offd_i.GetJ() : NULL;
-
-   const double * offd_r_D = (A_r) ? offd_r.GetData() : NULL;
-   const double * offd_i_D = (A_i) ? offd_i.GetData() : NULL;
-
-   int offd_r_nnz = (offd_r_I) ? offd_r_I[nrows] : 0;
-   int offd_i_nnz = (offd_i_I) ? offd_i_I[nrows] : 0;
-   int offd_nnz = 2 * (offd_r_nnz + offd_i_nnz);
-
-   // Allocate CSR arrays for the combined matrix
-   HYPRE_Int * diag_I = hypre_CTAlloc(HYPRE_Int, 2 * nrows + 1);
-   HYPRE_Int * diag_J = hypre_CTAlloc(HYPRE_Int, diag_nnz);
-   double    * diag_D = hypre_CTAlloc(double, diag_nnz);
-
-   HYPRE_Int * offd_I = hypre_CTAlloc(HYPRE_Int, 2 * nrows + 1);
-   HYPRE_Int * offd_J = hypre_CTAlloc(HYPRE_Int, offd_nnz);
-   double    * offd_D = hypre_CTAlloc(double, offd_nnz);
-   HYPRE_Int * cmap   = hypre_CTAlloc(HYPRE_Int, 2 * num_cols_offd);
-
-   // Fill the CSR arrays for the diagonal portion of the matrix
-   const double factor = (convention_ == HERMITIAN) ? 1.0 : -1.0;
-
-   diag_I[0] = 0;
-   diag_I[nrows] = diag_r_nnz + diag_i_nnz;
-   for (int i=0; i<nrows; i++)
-   {
-      diag_I[i + 1]         = ((diag_r_I)?diag_r_I[i+1]:0) +
-                              ((diag_i_I)?diag_i_I[i+1]:0);
-      diag_I[i + nrows + 1] = diag_I[i+1] + diag_r_nnz + diag_i_nnz;
-
-      if (diag_r_I)
-      {
-         for (int j=0; j<diag_r_I[i+1] - diag_r_I[i]; j++)
-         {
-            diag_J[diag_I[i] + j] = diag_r_J[diag_r_I[i] + j];
-            diag_D[diag_I[i] + j] = diag_r_D[diag_r_I[i] + j];
-
-            diag_J[diag_I[i+nrows] + j] =
-               diag_r_J[diag_r_I[i] + j] + ncols;
-            diag_D[diag_I[i+nrows] + j] =
-               factor * diag_r_D[diag_r_I[i] + j];
-         }
-      }
-      if (diag_i_I)
-      {
-         const int off_r = (diag_r_I)?(diag_r_I[i+1] - diag_r_I[i]):0;
-         for (int j=0; j<diag_i_I[i+1] - diag_i_I[i]; j++)
-         {
-            diag_J[diag_I[i] + off_r + j] =  diag_i_J[diag_i_I[i] + j] + ncols;
-            diag_D[diag_I[i] + off_r + j] = -diag_i_D[diag_i_I[i] + j];
-
-            diag_J[diag_I[i+nrows] + off_r + j] = diag_i_J[diag_i_I[i] + j];
-            diag_D[diag_I[i+nrows] + off_r + j] =
-               factor * diag_i_D[diag_i_I[i] + j];
-         }
-      }
-   }
-
-   // Determine the mappings describing the layout of off-diagonal columns
-   int num_recv_procs = 0;
-   HYPRE_Int * offd_col_start_stop = NULL;
-   this->getColStartStop(A_r, A_i, num_recv_procs, offd_col_start_stop);
-
-   std::set<int>::iterator sit;
-   std::map<int,int> cmapa, cmapb, cinvmap;
-   for (sit=cset.begin(); sit!=cset.end(); sit++)
-   {
-      int col_orig = *sit;
-      int col_2x2  = -1;
-      int col_size = 0;
-      for (int i=0; i<num_recv_procs; i++)
-      {
-         if (offd_col_start_stop[2*i] <= col_orig &&
-             col_orig < offd_col_start_stop[2*i+1])
-         {
-            col_2x2 = offd_col_start_stop[2*i] + col_orig;
-            col_size = offd_col_start_stop[2*i+1] - offd_col_start_stop[2*i];
-            break;
-         }
-      }
-      cmapa[*sit] = col_2x2;
-      cmapb[*sit] = col_2x2 + col_size;
-      cinvmap[col_2x2] = -1;
-      cinvmap[col_2x2 + col_size] = -1;
-   }
-   delete [] offd_col_start_stop;
-
-   std::map<int, int>::iterator mit;
-   int i = 0;
-   for (mit=cinvmap.begin(); mit!=cinvmap.end(); mit++, i++)
-   {
-      mit->second = i;
-      cmap[i] = mit->first;
-   }
-
-   // Fill the CSR arrays for the off-diagonal portion of the matrix
-   offd_I[0] = 0;
-   offd_I[nrows] = offd_r_nnz + offd_i_nnz;
-   for (int i=0; i<nrows; i++)
-   {
-      offd_I[i + 1]         = ((offd_r_I)?offd_r_I[i+1]:0) +
-                              ((offd_i_I)?offd_i_I[i+1]:0);
-      offd_I[i + nrows + 1] = offd_I[i+1] + offd_r_nnz + offd_i_nnz;
-
-      if (offd_r_I)
-      {
-         const int off_i = (offd_i_I)?(offd_i_I[i+1] - offd_i_I[i]):0;
-         for (int j=0; j<offd_r_I[i+1] - offd_r_I[i]; j++)
-         {
-            offd_J[offd_I[i] + j] =
-               cinvmap[cmapa[cmap_r[offd_r_J[offd_r_I[i] + j]]]];
-            offd_D[offd_I[i] + j] = offd_r_D[offd_r_I[i] + j];
-
-            offd_J[offd_I[i+nrows] + off_i + j] =
-               cinvmap[cmapb[cmap_r[offd_r_J[offd_r_I[i] + j]]]];
-            offd_D[offd_I[i+nrows] + off_i + j] =
-               factor * offd_r_D[offd_r_I[i] + j];
-         }
-      }
-      if (offd_i_I)
-      {
-         const int off_r = (offd_r_I)?(offd_r_I[i+1] - offd_r_I[i]):0;
-         for (int j=0; j<offd_i_I[i+1] - offd_i_I[i]; j++)
-         {
-            offd_J[offd_I[i] + off_r + j] =
-               cinvmap[cmapb[cmap_i[offd_i_J[offd_i_I[i] + j]]]];
-            offd_D[offd_I[i] + off_r + j] = -offd_i_D[offd_i_I[i] + j];
-
-            offd_J[offd_I[i+nrows] + j] =
-               cinvmap[cmapa[cmap_i[offd_i_J[offd_i_I[i] + j]]]];
-            offd_D[offd_I[i+nrows] + j] = factor * offd_i_D[offd_i_I[i] + j];
-         }
-      }
-   }
-
-   // Construct the combined matrix
-   HypreParMatrix * A = new HypreParMatrix(comm_,
-                                           2 * global_num_rows,
-                                           2 * global_num_cols,
-                                           row_starts, col_starts,
-                                           diag_I, diag_J, diag_D,
-                                           offd_I, offd_J, offd_D,
-                                           2 * num_cols_offd, cmap);
-
-   // Give the new matrix ownership of its interanl arrays
-   A->SetOwnerFlags(-1,-1,-1);
-   hypre_CSRMatrixSetDataOwner(((hypre_ParCSRMatrix*)(*A))->diag,1);
-   hypre_CSRMatrixSetDataOwner(((hypre_ParCSRMatrix*)(*A))->offd,1);
-   hypre_ParCSRMatrixSetRowStartsOwner((hypre_ParCSRMatrix*)(*A),1);
-   hypre_ParCSRMatrixSetColStartsOwner((hypre_ParCSRMatrix*)(*A),1);
-
-   return A;
-}
-
-void
-ComplexHypreParMatrix::getColStartStop(const HypreParMatrix * A_r,
-                                       const HypreParMatrix * A_i,
-                                       int & num_recv_procs,
-                                       HYPRE_Int *& offd_col_start_stop) const
-{
-   hypre_ParCSRCommPkg * comm_pkg_r =
-      (A_r) ? hypre_ParCSRMatrixCommPkg((hypre_ParCSRMatrix*)(*A_r)) : NULL;
-   hypre_ParCSRCommPkg * comm_pkg_i =
-      (A_i) ? hypre_ParCSRMatrixCommPkg((hypre_ParCSRMatrix*)(*A_i)) : NULL;
-
-   std::set<HYPRE_Int> send_procs, recv_procs;
-   if ( comm_pkg_r )
-   {
-      for (HYPRE_Int i=0; i<comm_pkg_r->num_sends; i++)
-      {
-         send_procs.insert(comm_pkg_r->send_procs[i]);
-      }
-      for (HYPRE_Int i=0; i<comm_pkg_r->num_recvs; i++)
-      {
-         recv_procs.insert(comm_pkg_r->recv_procs[i]);
-      }
-   }
-   if ( comm_pkg_i )
-   {
-      for (HYPRE_Int i=0; i<comm_pkg_i->num_sends; i++)
-      {
-         send_procs.insert(comm_pkg_i->send_procs[i]);
-      }
-      for (HYPRE_Int i=0; i<comm_pkg_i->num_recvs; i++)
-      {
-         recv_procs.insert(comm_pkg_i->recv_procs[i]);
-      }
-   }
-
-   num_recv_procs = (int)recv_procs.size();
-
-   HYPRE_Int loc_start_stop[2];
-   offd_col_start_stop = new HYPRE_Int[2 * num_recv_procs];
-
-   const HYPRE_Int * row_part = (A_r) ? A_r->RowPart() :
-                                ((A_i) ? A_i->RowPart() : NULL);
-
-   int row_part_ind = (HYPRE_AssumedPartitionCheck()) ? 0 : myid_;
-   loc_start_stop[0] = row_part[row_part_ind];
-   loc_start_stop[1] = row_part[row_part_ind+1];
-
-   MPI_Request * req = new MPI_Request[send_procs.size()+recv_procs.size()];
-   MPI_Status * stat = new MPI_Status[send_procs.size()+recv_procs.size()];
-   int send_count = 0;
-   int recv_count = 0;
-   int tag = 0;
-
-   std::set<HYPRE_Int>::iterator sit;
-   for (sit=send_procs.begin(); sit!=send_procs.end(); sit++)
-   {
-      MPI_Isend(loc_start_stop, 2, HYPRE_MPI_INT,
-                *sit, tag, comm_, &req[send_count]);
-      send_count++;
-   }
-   for (sit=recv_procs.begin(); sit!=recv_procs.end(); sit++)
-   {
-      MPI_Irecv(&offd_col_start_stop[2*recv_count], 2, HYPRE_MPI_INT,
-                *sit, tag, comm_, &req[send_count+recv_count]);
-      recv_count++;
-   }
-
-   MPI_Waitall(send_count+recv_count, req, stat);
-
-   delete [] req;
-   delete [] stat;
-}
-
-
-#endif // MFEM_USE_MPI
-
 }

linalg/complex_operator.hpp

@@ -14,9 +14,6 @@
 
 #include "operator.hpp"
 #include "sparsemat.hpp"
-#ifdef MFEM_USE_MPI
-#include "hypre.hpp"
-#endif
 
 namespace mfem
 {
@@ -30,8 +27,7 @@ namespace mfem
     ComplexOperator allows one to choose a convention upon construction, which
     facilitates symmetry.
 
-    If we let (y_r + i y_i) = (Op_r + i Op_i)(x_r + i x_i) then Matrix-vector
-    products are then computed as:
+    Matrix-vector products are then computed as:
 
     1. When Convention::HERMITIAN is used (default)
     / y_r \   / Op_r -Op_i \ / x_r \
@@ -42,8 +38,6 @@ namespace mfem
     / y_r \   / Op_r -Op_i \ / x_r \
     |     | = |            | |     |
     \-y_i /   \-Op_i -Op_r / \ x_i /
-    In other words, Matrix-vector products with Convention::BLOCK_SYMMETRIC
-    compute the complex conjugate of Op*x.
 
     Either convention can be used with a given complex operator,
     however, each of them is best suited for certain classes of
@@ -88,30 +82,9 @@ public:
 
    virtual ~ComplexOperator();
 
-   /** @brief Check for existence of real or imaginary part of the operator
-
-       These methods do not check that the operators are non-zero but
-       only that the operators have been set.
-    */
-   bool hasRealPart() const { return Op_Real_ != NULL; }
-   bool hasImagPart() const { return Op_Imag_ != NULL; }
-
-   /** @brief Real or imaginary part accessor methods
-
-       The following accessor methods should only be called if the
-       requested part of the opertor is known to exist.  This
-       can be checked with hasRealPart() or hasImagPart().
-   */
-   virtual Operator & real();
-   virtual Operator & imag();
-   virtual const Operator & real() const;
-   virtual const Operator & imag() const;
-
    virtual void Mult(const Vector &x, Vector &y) const;
    virtual void MultTranspose(const Vector &x, Vector &y) const;
 
-   virtual Type GetType() const { return Complex_Operator; }
-
 protected:
    // Let this be hidden from the public interface since the implementation
    // depends on internal members
@@ -154,74 +127,9 @@ public:
       : ComplexOperator(A_Real, A_Imag, ownReal, ownImag, convention)
    {}
 
-   virtual SparseMatrix & real();
-   virtual SparseMatrix & imag();
-
-   virtual const SparseMatrix & real() const;
-   virtual const SparseMatrix & imag() const;
-
-   /** Combine the blocks making up this complex operator into a
-       single SparseMatrix.  The resulting matrix can be passed to
-       solvers which require access to the matrix entries themselves,
-       such as sparse direct solvers, rather than simply the action of
-       the opertor.  Note that this combined operator requires roughly
-       twice the memory of the block structured operator. */
    SparseMatrix * GetSystemMatrix() const;
-
-   virtual Type GetType() const { return MFEM_ComplexSparseMat; }
 };
 
-#ifdef MFEM_USE_MPI
-
-/** @brief Specialization of the ComplexOperator built from a pair of
-    HypreParMatrices.
-
-    The purpose of this specialization is to construct a single
-    HypreParMatrix object which is equivalent to the 2x2 block system
-    that the ComplexOperator mimics. The resulting HypreParMatrix can
-    then be passed along to solvers which require access to the CSR
-    matrix data such as SuperLU, STRUMPACK, or similar sparse linear
-    solvers.
-
-    See ComplexOperator documentation in operator.hpp for more information.
- */
-class ComplexHypreParMatrix : public ComplexOperator
-{
-public:
-   ComplexHypreParMatrix(HypreParMatrix * A_Real, HypreParMatrix * A_Imag,
-                         bool ownReal, bool ownImag,
-                         Convention convention = HERMITIAN);
-
-   virtual HypreParMatrix & real();
-   virtual HypreParMatrix & imag();
-
-   virtual const HypreParMatrix & real() const;
-   virtual const HypreParMatrix & imag() const;
-
-   /** Combine the blocks making up this complex operator into a
-       single HypreParMatrix.  The resulting matrix can be passed to
-       solvers which require access to the matrix entries themselves,
-       such as sparse direct solvers or Hypre preconditioners, rather
-       than simply the action of the opertor.  Note that this combined
-       operator requires roughly twice the memory of the block
-       structured operator. */
-   HypreParMatrix * GetSystemMatrix() const;
-
-   virtual Type GetType() const { return Complex_Hypre_ParCSR; }
-
-private:
-   void getColStartStop(const HypreParMatrix * A_r,
-                        const HypreParMatrix * A_i,
-                        int & num_recv_procs,
-                        HYPRE_Int *& offd_col_start_stop) const;
-
-   MPI_Comm comm_;
-   int myid_;
-   int nranks_;
-};
-
-#endif // MFEM_USE_MPI
-
 }
 
-#endif // MFEM_COMPLEX_OPERATOR
+#endif

linalg/hypre_parcsr.cpp

@@ -1520,10 +1520,14 @@ hypre_ParCSRMatrixSum(hypre_ParCSRMatrix *A,
    hypre_CSRMatrix *A_offd = hypre_ParCSRMatrixOffd(A);
    hypre_CSRMatrix *B_diag = hypre_ParCSRMatrixDiag(B);
    hypre_CSRMatrix *B_offd = hypre_ParCSRMatrixOffd(B);
+   HYPRE_Int ncols_B_offd  = hypre_CSRMatrixNumCols(B_offd);
    HYPRE_Int error;
 
    error = hypre_CSRMatrixSum(A_diag, beta, B_diag);
-   error = error ? error : hypre_CSRMatrixSum(A_offd, beta, B_offd);
+   if (ncols_B_offd > 0) /* treat B_offd as zero if it has no columns */
+   {
+      error = error ? error : hypre_CSRMatrixSum(A_offd, beta, B_offd);
+   }
 
    return error;
 }

linalg/ode.cpp

@@ -561,6 +561,102 @@ void GeneralizedAlphaSolver::Step(Vector &x, double &t, double &dt)
 }
 
 
+void IMEX_BE_FE::Init(TimeDependentOperator &_f)
+{
+   ODESolver::Init(_f);
+   k_imp.SetSize(f->Width());
+   y.SetSize(f->Width());
+   k_exp.SetSize(f->Width());
+}
+
+void IMEX_BE_FE::Step(Vector &x, double &t, double &dt)
+{
+   f->ExplicitMult(x, k_exp);
+   add(x, dt, k_exp, y);
+
+   f->SetTime(t + dt);
+   f->ImplicitSolve(dt, y, k_imp);
+
+   x.Add(dt, k_exp);
+   x.Add(dt, k_imp);
+   t += dt;
+}
+
+
+void IMEXRK2::Init(TimeDependentOperator &_f)
+{
+   ODESolver::Init(_f);
+   f = ODESolver::f;
+   k_imp.SetSize(f->Width());
+   k_exp.SetSize(f->Width());
+   y.SetSize(f->Width());
+   z.SetSize(f->Width());
+}
+
+void IMEXRK2::Step(Vector &x, double &t, double &dt)
+{
+   double gamma = 1 - sqrt(2)/2;
+   double delta = -2*sqrt(2)/3;
+
+   // The method is given by
+   // k1_exp = f(u)
+   // k1_imp = g(u + gamma*dt*k1_exp + gamma*dt*k1_imp)
+   // k2_exp = f(u + gamma*dt*k1_exp + gamma*dt*k1_imp)
+   // k2_imp = g(u + delta*dt*k1_exp + (1-gamma)*dt*k1_imp
+   //              + (1-delta)*dt*k2_exp + gamma*dt*k2_imp)
+   // k3_exp = f(u + delta*dt*k1_exp + (1-gamma)*dt*k1_imp
+   //              + (1-delta)*dt*k2_exp + gamma*dt*k2_imp)
+   // u_new = u + dt*((1-gamma)*k1_imp + (1-gamma)*k2_exp
+   //                 + gamma*k2_imp + gamma*k3_exp)
+
+   // Take first explicit step
+   // k1_exp = f(u)
+   f->ExplicitMult(x, k_exp);
+   // b corresponding to this stage is zero, so don't add to solution
+
+   // Solve first implicit step
+   // y = u + gamma*dt*k1_exp
+   add(x, gamma*dt, k_exp, y);
+   // Solve x1_imp = g(u + gamma*dt*k1_exp + gamma*dt*k1_imp)
+   f->SetTime(t + gamma*dt);
+   f->ImplicitSolve(gamma*dt, y, k_imp);
+   // x = u + (1-gamma)*dt*k1_imp
+   x.Add((1-gamma)*dt, k_imp);
+
+   // Begin setting up rhs for second solve
+   // z = u + (1-gamma)*dt*k_imp + delta*dt*k_exp
+   add(x, delta*dt, k_exp, z);
+
+   // Take second explicit step
+   // y = x + gamma*dt*k1_exp + gamma*dt*k1_imp
+   y.Add(gamma*dt, k_imp);
+   // k2_exp = f(x + gamma*dt*k1_exp + gamma*dt*k1_imp)
+   f->ExplicitMult(y, k_exp);
+   // x = u + (1-gamma)*dt*k1_imp + (1-gamma)*dt*k2_exp
+   x.Add((1-gamma)*dt, k_exp);
+
+   // Finish formoing rhs
+   // z = x + (1-gamma)*dt*k1_imp + delta*dt*k1_exp + (1-delta)*dt*k2_exp
+   z.Add((1-delta)*dt, k_exp);
+
+   // Solve second implicit step for k2_imp
+   f->SetTime(t + dt);
+   f->ImplicitSolve(gamma*dt, z, k_imp);
+   // x = u + (1-gamma)*dt*k1_imp + (1-gamma)*dt*k2_exp + gamma*dt*k2_imp
+   x.Add(gamma*dt, k_imp);
+
+   // Take final explicit step for k3_exp
+   z.Add(gamma*dt, k_imp);
+   f->ExplicitMult(z, k_exp);
+
+   // x = u + (1-gamma)*dt*k1_imp + (1-gamma)*dt*k2_exp + gamma*dt*k2_imp
+   //       + gamma*dt*k3_exp
+   x.Add(gamma*dt, k_exp);
+
+   t += dt;
+}
+
+
 void
 SIASolver::Init(Operator &P, TimeDependentOperator & F)
 {

linalg/ode.hpp

@@ -305,6 +305,33 @@ public:
 };
 
 
+/// IMEX Backward-Forward Euler ODE solver
+class IMEX_BE_FE : public ODESolver
+{
+protected:
+   Vector k_exp, k_imp, y;
+
+public:
+   virtual void Init(TimeDependentOperator &_f);
+
+   virtual void Step(Vector &x, double &t, double &dt);
+};
+
+/** Second-order IMEX (2,3,2) method, from "Implicit-explicit Runge-Kutta
+    methods for time-dependent partial differential equations" by Ascher, Ruuth
+    and Spiteri, Applied Numerical Mathematics (1997). */
+class IMEXRK2 : public ODESolver
+{
+protected:
+   Vector k_exp, k_imp, y, z;
+
+public:
+   virtual void Init(TimeDependentOperator &_f);
+
+   virtual void Step(Vector &x, double &t, double &dt);
+};
+
+
 /// The SIASolver class is based on the Symplectic Integration Algorithm
 /// described in "A Symplectic Integration Algorithm for Separable Hamiltonian
 /// Functions" by J. Candy and W. Rozmus, Journal of Computational Physics,

linalg/operator.hpp

@@ -124,17 +124,14 @@ public:
    enum Type
    {
       ANY_TYPE,         ///< ID for the base class Operator, i.e. any type.
-      MFEM_SPARSEMAT,   ///< ID for class SparseMatrix.
+      MFEM_SPARSEMAT,   ///< ID for class SparseMatrix
       Hypre_ParCSR,     ///< ID for class HypreParMatrix.
       PETSC_MATAIJ,     ///< ID for class PetscParMatrix, MATAIJ format.
       PETSC_MATIS,      ///< ID for class PetscParMatrix, MATIS format.
       PETSC_MATSHELL,   ///< ID for class PetscParMatrix, MATSHELL format.
       PETSC_MATNEST,    ///< ID for class PetscParMatrix, MATNEST format.
       PETSC_MATHYPRE,   ///< ID for class PetscParMatrix, MATHYPRE format.
-      PETSC_MATGENERIC, ///< ID for class PetscParMatrix, unspecified format.
-      Complex_Operator, ///< ID for class ComplexOperator.
-      MFEM_ComplexSparseMat, ///< ID for class ComplexSparseMatrix.
-      Complex_Hypre_ParCSR   ///< ID for class ComplexHypreParMatrix.
+      PETSC_MATGENERIC  ///< ID for class PetscParMatrix, unspecified format.
    };
 
    /// Return the type ID of the Operator class.
@@ -305,24 +302,6 @@ public:
 };
 
 
-/// Scaled Operator B: x -> a A(x).
-class ScaledOperator : public Operator
-{
-private:
-   const Operator &A_;
-   double a_;
-
-public:
-   /// Create a scalar product operator related to A.
-   explicit ScaledOperator(const Operator *A, double a)
-      : Operator(A->Width(), A->Height()), A_(*A), a_(a) { }
-
-   /// Operator application
-   virtual void Mult(const Vector &x, Vector &y) const
-   { A_.Mult(x, y); y *= a_; }
-};
-
-
 /** @brief The transpose of a given operator. Switches the roles of the methods
     Mult() and MultTranspose(). */
 class TransposeOperator : public Operator

linalg/sparsemat.cpp

@@ -38,6 +38,7 @@ SparseMatrix::SparseMatrix(int nrows, int ncols)
      current_row(-1),
      ColPtrJ(NULL),
      ColPtrNode(NULL),
+     At(NULL),
      ownGraph(true),
      ownData(true),
      isSorted(false)
@@ -60,6 +61,7 @@ SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n)
      Rows(NULL),
      ColPtrJ(NULL),
      ColPtrNode(NULL),
+     At(NULL),
      ownGraph(true),
      ownData(true),
      isSorted(false)
@@ -78,6 +80,7 @@ SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n,
      Rows(NULL),
      ColPtrJ(NULL),
      ColPtrNode(NULL),
+     At(NULL),
      ownGraph(ownij),
      ownData(owna),
      isSorted(issorted)
@@ -103,6 +106,7 @@ SparseMatrix::SparseMatrix(int nrows, int ncols, int rowsize)
    , Rows(NULL)
    , ColPtrJ(NULL)
    , ColPtrNode(NULL)
+   , At(NULL)
    , ownGraph(true)
    , ownData(true)
    , isSorted(false)
@@ -183,6 +187,7 @@ SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph)
    current_row = -1;
    ColPtrJ = NULL;
    ColPtrNode = NULL;
+   At = NULL;
    isSorted = mat.isSorted;
 }
 
@@ -191,6 +196,7 @@ SparseMatrix::SparseMatrix(const Vector &v)
    , Rows(NULL)
    , ColPtrJ(NULL)
    , ColPtrNode(NULL)
+   , At(NULL)
    , ownGraph(true)
    , ownData(true)
    , isSorted(true)
@@ -245,6 +251,7 @@ void SparseMatrix::SetEmpty()
    current_row = -1;
    ColPtrJ = NULL;
    ColPtrNode = NULL;
+   At = NULL;
 #ifdef MFEM_USE_MEMALLOC
    NodesMem = NULL;
 #endif
@@ -333,7 +340,7 @@ void SparseMatrix::SetWidth(int newWidth)
       // Nothing to be done here
       return;
    }
-   else if ( newWidth == -1)
+   else if (newWidth == -1)
    {
       // Compute the actual width
       width = ActualWidth();
@@ -552,12 +559,10 @@ void SparseMatrix::Mult(const Vector &x, Vector &y) const
 
 void SparseMatrix::AddMult(const Vector &x, Vector &y, const double a) const
 {
-   MFEM_ASSERT(width == x.Size(),
-               "Input vector size (" << x.Size() << ") must match matrix width (" << width
-               << ")");
-   MFEM_ASSERT(height == y.Size(),
-               "Output vector size (" << y.Size() << ") must match matrix height (" << height
-               << ")");
+   MFEM_ASSERT(width == x.Size(), "Input vector size (" << x.Size()
+               << ") must match matrix width (" << width << ")");
+   MFEM_ASSERT(height == y.Size(), "Output vector size (" << y.Size()
+               << ") must match matrix height (" << height << ")");
 
    int i, j, end;
    double *Ap = A, *yp = y.GetData();
@@ -636,12 +641,10 @@ void SparseMatrix::MultTranspose(const Vector &x, Vector &y) const
 void SparseMatrix::AddMultTranspose(const Vector &x, Vector &y,
                                     const double a) const
 {
-   MFEM_ASSERT(height == x.Size(),
-               "Input vector size (" << x.Size() << ") must match matrix height (" << height
-               << ")");
-   MFEM_ASSERT(width == y.Size(),
-               "Output vector size (" << y.Size() << ") must match matrix width (" << width
-               << ")");
+   MFEM_ASSERT(height == x.Size(), "Input vector size (" << x.Size()
+               << ") must match matrix height (" << height << ")");
+   MFEM_ASSERT(width == y.Size(), "Output vector size (" << y.Size()
+               << ") must match matrix width (" << width << ")");
 
    if (A == NULL)
    {
@@ -658,23 +661,40 @@ void SparseMatrix::AddMultTranspose(const Vector &x, Vector &y,
       }
       return;
    }
-   // Prepare the lambda capture and get our pointers from the memory manager
-   const int d_height = height;
-   const DeviceArray d_I(I);
-   const DeviceArray d_J(J);
-   const DeviceVector d_A(A);
-   const DeviceVector d_x(x, x.Size());
-   DeviceVector d_y(y, y.Size());
-   MFEM_FORALL(i, d_height,
+
+   if (At)
    {
-      const double xi = a * d_x[i];
-      const int end = d_I[i+1];
-      for (int j = d_I[i]; j < end; j++)
+      At->AddMult(x, y, a);
+   }
+   else
+   {
+      MFEM_VERIFY(Device::IsDisabled(), "transpose action on device is not "
+                  "enabled; see BuildTranspose() for details.");
+      for (int i = 0; i < height; i++)
       {
-         const int Jj = d_J[j];
-         AtomicAdd(&d_y[Jj], d_A[j] * xi);
+         const double xi = a * x[i];
+         const int end = I[i+1];
+         for (int j = I[i]; j < end; j++)
+         {
+            const int Jj = J[j];
+            y[Jj] += A[j] * xi;
+         }
       }
-   });
+   }
+}
+
+void SparseMatrix::BuildTranspose() const
+{
+   if (At == NULL)
+   {
+      At = Transpose(*this);
+   }
+}
+
+void SparseMatrix::ResetTranspose() const
+{
+   delete At;
+   At = NULL;
 }
 
 void SparseMatrix::PartMult(
@@ -2101,12 +2121,12 @@ void SparseMatrix::Set(const int i, const int j, const double A)
    if ((gi=i) < 0) { gi = -1-gi, s = -1; }
    else { s = 1; }
    MFEM_ASSERT(gi < height,
-               "Trying to insert a row " << gi << " outside the matrix height "
+               "Trying to set a row " << gi << " outside the matrix height "
                << height);
    if ((gj=j) < 0) { gj = -1-gj, t = -s; }
    else { t = s; }
    MFEM_ASSERT(gj < width,
-               "Trying to insert a column " << gj << " outside the matrix width "
+               "Trying to set a column " << gj << " outside the matrix width "
                << width);
    if (t < 0) { a = -a; }
    _Set_(gi, gj, a);
@@ -2142,7 +2162,7 @@ void SparseMatrix::SetSubMatrix(const Array<int> &rows, const Array<int> &cols,
       if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
       else { s = 1; }
       MFEM_ASSERT(gi < height,
-                  "Trying to insert a row " << gi << " outside the matrix height "
+                  "Trying to set a row " << gi << " outside the matrix height "
                   << height);
       SetColPtr(gi);
       for (j = 0; j < cols.Size(); j++)
@@ -2155,7 +2175,7 @@ void SparseMatrix::SetSubMatrix(const Array<int> &rows, const Array<int> &cols,
          if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
          else { t = s; }
          MFEM_ASSERT(gj < width,
-                     "Trying to insert a column " << gj << " outside the matrix width "
+                     "Trying to set a column " << gj << " outside the matrix width "
                      << width);
          if (t < 0) { a = -a; }
          _Set_(gj, a);
@@ -2177,7 +2197,7 @@ void SparseMatrix::SetSubMatrixTranspose(const Array<int> &rows,
       if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
       else { s = 1; }
       MFEM_ASSERT(gi < height,
-                  "Trying to insert a row " << gi << " outside the matrix height "
+                  "Trying to set a row " << gi << " outside the matrix height "
                   << height);
       SetColPtr(gi);
       for (j = 0; j < cols.Size(); j++)
@@ -2190,7 +2210,7 @@ void SparseMatrix::SetSubMatrixTranspose(const Array<int> &rows,
          if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
          else { t = s; }
          MFEM_ASSERT(gj < width,
-                     "Trying to insert a column " << gj << " outside the matrix width "
+                     "Trying to set a column " << gj << " outside the matrix width "
                      << width);
          if (t < 0) { a = -a; }
          _Set_(gj, a);
@@ -2210,7 +2230,7 @@ void SparseMatrix::GetSubMatrix(const Array<int> &rows, const Array<int> &cols,
       if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
       else { s = 1; }
       MFEM_ASSERT(gi < height,
-                  "Trying to insert a row " << gi << " outside the matrix height "
+                  "Trying to read a row " << gi << " outside the matrix height "
                   << height);
       SetColPtr(gi);
       for (j = 0; j < cols.Size(); j++)
@@ -2218,7 +2238,7 @@ void SparseMatrix::GetSubMatrix(const Array<int> &rows, const Array<int> &cols,
          if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
          else { t = s; }
          MFEM_ASSERT(gj < width,
-                     "Trying to insert a column " << gj << " outside the matrix width "
+                     "Trying to read a column " << gj << " outside the matrix width "
                      << width);
          a = _Get_(gj);
          subm(i, j) = (t < 0) ? (-a) : (a);
@@ -2236,7 +2256,7 @@ bool SparseMatrix::RowIsEmpty(const int row) const
       gi = -1-gi;
    }
    MFEM_ASSERT(gi < height,
-               "Trying to insert a row " << gi << " outside the matrix height "
+               "Trying to query a row " << gi << " outside the matrix height "
                << height);
    if (Rows)
    {
@@ -2255,7 +2275,7 @@ int SparseMatrix::GetRow(const int row, Array<int> &cols, Vector &srow) const
 
    if ((gi=row) < 0) { gi = -1-gi; }
    MFEM_ASSERT(gi < height,
-               "Trying to insert a row " << gi << " outside the matrix height "
+               "Trying to read a row " << gi << " outside the matrix height "
                << height);
    if (Rows)
    {
@@ -2282,7 +2302,7 @@ int SparseMatrix::GetRow(const int row, Array<int> &cols, Vector &srow) const
       j = I[gi];
       cols.MakeRef(J + j, I[gi+1]-j);
       srow.NewDataAndSize(A + j, cols.Size());
-      MFEM_ASSERT(row >= 0, "Row not valid: " << row );
+      MFEM_ASSERT(row >= 0, "Row not valid: " << row << ", height: " << height);
       return 1;
    }
 }
@@ -2296,7 +2316,7 @@ void SparseMatrix::SetRow(const int row, const Array<int> &cols,
    if ((gi=row) < 0) { gi = -1-gi, s = -1; }
    else { s = 1; }
    MFEM_ASSERT(gi < height,
-               "Trying to insert a row " << gi << " outside the matrix height "
+               "Trying to set a row " << gi << " outside the matrix height "
                << height);
 
    if (!Finalized())
@@ -2307,7 +2327,7 @@ void SparseMatrix::SetRow(const int row, const Array<int> &cols,
          if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
          else { t = s; }
          MFEM_ASSERT(gj < width,
-                     "Trying to insert a column " << gj << " outside the matrix"
+                     "Trying to set a column " << gj << " outside the matrix"
                      " width " << width);
          a = srow(j);
          if (t < 0) { a = -a; }
@@ -2325,7 +2345,7 @@ void SparseMatrix::SetRow(const int row, const Array<int> &cols,
          if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
          else { t = s; }
          MFEM_ASSERT(gj < width,
-                     "Trying to insert a column " << gj << " outside the matrix"
+                     "Trying to set a column " << gj << " outside the matrix"
                      " width " << width);
 
          J[i] = gj;
@@ -2771,9 +2791,10 @@ void SparseMatrix::Destroy()
       delete NodesMem;
    }
 #endif
+   delete At;
 }
 
-int SparseMatrix::ActualWidth()
+int SparseMatrix::ActualWidth() const
 {
    int awidth = 0;
    if (A)
@@ -2817,8 +2838,10 @@ SparseMatrix *Transpose (const SparseMatrix &A)
       "Finalize must be called before Transpose. Use TransposeRowMatrix instead");
 
    int i, j, end;
-   int m, n, nnz, *A_i, *A_j, *At_i, *At_j;
-   double *A_data, *At_data;
+   const int *A_i, *A_j;
+   int m, n, nnz, *At_i, *At_j;
+   const double *A_data;
+   double *At_data;
 
    m      = A.Height(); // number of rows of A
    n      = A.Width();  // number of columns of A
@@ -2945,8 +2968,10 @@ SparseMatrix *Mult (const SparseMatrix &A, const SparseMatrix &B,
                     SparseMatrix *OAB)
 {
    int nrowsA, ncolsA, nrowsB, ncolsB;
-   int *A_i, *A_j, *B_i, *B_j, *C_i, *C_j, *B_marker;
-   double *A_data, *B_data, *C_data;
+   const int *A_i, *A_j, *B_i, *B_j;
+   int *C_i, *C_j, *B_marker;
+   const double *A_data, *B_data;
+   double *C_data;
    int ia, ib, ic, ja, jb, num_nonzeros;
    int row_start, counter;
    double a_entry, b_entry;
@@ -3259,13 +3284,13 @@ SparseMatrix * Add(double a, const SparseMatrix & A, double b,
    int * C_j;
    double * C_data;
 
-   int * A_i = A.GetI();
-   int * A_j = A.GetJ();
-   double * A_data = A.GetData();
+   const int *A_i = A.GetI();
+   const int *A_j = A.GetJ();
+   const double *A_data = A.GetData();
 
-   int * B_i = B.GetI();
-   int * B_j = B.GetJ();
-   double * B_data = B.GetData();
+   const int *B_i = B.GetI();
+   const int *B_j = B.GetJ();
+   const double *B_data = B.GetData();
 
    int * marker = new int[ncols];
    std::fill(marker, marker+ncols, -1);
@@ -3493,6 +3518,7 @@ void SparseMatrix::Swap(SparseMatrix &other)
    mfem::Swap(current_row, other.current_row);
    mfem::Swap(ColPtrJ, other.ColPtrJ);
    mfem::Swap(ColPtrNode, other.ColPtrNode);
+   mfem::Swap(At, other.At);
 
 #ifdef MFEM_USE_MEMALLOC
    mfem::Swap(NodesMem, other.NodesMem);

linalg/sparsemat.hpp

@@ -64,6 +64,9 @@ protected:
    mutable int* ColPtrJ;
    mutable RowNode ** ColPtrNode;
 
+   /// Transpose of A. Owned. Used to perform MultTranspose() on devices.
+   mutable SparseMatrix *At;
+
 #ifdef MFEM_USE_MEMALLOC
    typedef MemAlloc <RowNode, 1024> RowNodeAlloc;
    RowNodeAlloc * NodesMem;
@@ -97,6 +100,9 @@ public:
 
    /** @brief Create a sparse matrix in CSR format. Ownership of @a i, @a j, and
        @a data is optionally transferred to the SparseMatrix. */
+   /** If the parameter @a data is NULL, then the internal #A array is allocated
+       by this constructor (initializing it with zeros and taking ownership,
+       regardless of the parameter @a owna). */
    SparseMatrix(int *i, int *j, double *data, int m, int n, bool ownij,
                 bool owna, bool issorted);
 
@@ -112,7 +118,7 @@ public:
        ownership. */
    SparseMatrix(const SparseMatrix &mat, bool copy_graph = true);
 
-   /// Create a SparseMatrix with diagonal v, i.e. A = Diag(v)
+   /// Create a SparseMatrix with diagonal @a v, i.e. A = Diag(v)
    SparseMatrix(const Vector & v);
 
 
@@ -134,21 +140,35 @@ public:
    /// Check if the SparseMatrix is empty.
    bool Empty() const { return (A == NULL) && (Rows == NULL); }
 
-   /// Return the array #I
-   inline int *GetI() const { return I; }
-   /// Return the array #J
-   inline int *GetJ() const { return J; }
-   /// Return element data, i.e. array #A
-   inline double *GetData() const { return A; }
-   /// Returns the number of elements in row @a i
+   /// Return the array #I.
+   inline int *GetI() { return I; }
+   /// Return the array #I, const version.
+   inline const int *GetI() const { return I; }
+
+   /// Return the array #J.
+   inline int *GetJ() { return J; }
+   /// Return the array #J, const version.
+   inline const int *GetJ() const { return J; }
+
+   /// Return the element data, i.e. the array #A.
+   inline double *GetData() { return A; }
+   /// Return the element data, i.e. the array #A, const version.
+   inline const double *GetData() const { return A; }
+
+   /// Returns the number of elements in row @a i.
    int RowSize(const int i) const;
-   /// Returns the maximum number of elements among all rows
+
+   /// Returns the maximum number of elements among all rows.
    int MaxRowSize() const;
-   /// Return a pointer to the column indices in a row
+
+   /// Return a pointer to the column indices in a row.
    int *GetRowColumns(const int row);
+   /// Return a pointer to the column indices in a row, const version.
    const int *GetRowColumns(const int row) const;
-   /// Return a pointer to the entries in a row
+
+   /// Return a pointer to the entries in a row.
    double *GetRowEntries(const int row);
+   /// Return a pointer to the entries in a row, const version.
    const double *GetRowEntries(const int row) const;
 
    /// Change the width of a SparseMatrix.
@@ -163,7 +183,7 @@ public:
 
    /// Returns the actual Width of the matrix.
    /*! This method can be called for matrices finalized or not. */
-   int ActualWidth();
+   int ActualWidth() const;
 
    /// Sort the column indices corresponding to each row.
    void SortColumnIndices();
@@ -206,13 +226,45 @@ public:
    void AddMultTranspose(const Vector &x, Vector &y,
                          const double a = 1.0) const;
 
+   /** @brief Build and store internally the transpose of this matrix which will
+       be used in the methods AddMultTranspose() and MultTranspose(). */
+   /** If this method has been called, the internal transpose matrix will be
+       used to perform the action of the transpose matrix in AddMultTranspose(),
+       and MultTranspose().
+
+       Warning: any changes in this matrix will invalidate the internal
+       transpose. To rebuild the transpose, call ResetTranspose() followed by a
+       call to this method. If the internal transpose is already built, this
+       method has no effect.
+
+       When any non-default backend is enabled, i.e. Device::IsEnabled() is
+       true, the methods AddMultTranspose(), and MultTranspose(), require the
+       internal transpose to be built. If that is not the case (i.e. the
+       internal transpose is not built), these methods will raise an error with
+       an appropriate message pointing to this method. When using the default
+       backend, calling this method is optional.
+
+       This method can only be used when the sparse matrix is finalized. */
+   void BuildTranspose() const;
+
+   /** Reset (destroy) the internal transpose matrix. See BuildTranspose() for
+       more details. */
+   void ResetTranspose() const;
+
    void PartMult(const Array<int> &rows, const Vector &x, Vector &y) const;
    void PartAddMult(const Array<int> &rows, const Vector &x, Vector &y,
                     const double a=1.0) const;
 
-   /// y = A * x, but treat all elements as booleans (zero=false, nonzero=true).
+   /// y = A * x, treating all entries as booleans (zero=false, nonzero=true).
+   /** The actual values stored in the data array, #A, are not used - this means
+       and that all entries in the sparsity pattern are considered to be true by
+       this method. */
    void BooleanMult(const Array<int> &x, Array<int> &y) const;
-   /// y = At * x, but treat all elements as booleans (zero=false, nonzero=true).
+
+   /// y = At * x, treating all entries as booleans (zero=false, nonzero=true).
+   /** The actual values stored in the data array, #A, are not used - this means
+       and that all entries in the sparsity pattern are considered to be true by
+       this method. */
    void BooleanMultTranspose(const Array<int> &x, Array<int> &y) const;
 
    /// Compute y^t A x
@@ -452,11 +504,11 @@ SparseMatrix *Transpose(const SparseMatrix &A);
 SparseMatrix *TransposeAbstractSparseMatrix (const AbstractSparseMatrix &A,
                                              int useActualWidth);
 
-/** Matrix product A.B.
-    If OAB is not NULL, we assume it has the structure
-    of A.B and store the result in OAB.
-    If OAB is NULL, we create a new SparseMatrix to store
+/// Matrix product A.B.
+/** If @a OAB is not NULL, we assume it has the structure of A.B and store the
+    result in @a OAB. If @a OAB is NULL, we create a new SparseMatrix to store
     the result and return a pointer to it.
+
     All matrices must be finalized. */
 SparseMatrix *Mult(const SparseMatrix &A, const SparseMatrix &B,
                    SparseMatrix *OAB = NULL);
@@ -555,16 +607,20 @@ inline void SparseMatrix::SetColPtr(const int row) const
 inline void SparseMatrix::ClearColPtr() const
 {
    if (Rows)
+   {
       for (RowNode *node_p = Rows[current_row]; node_p != NULL;
            node_p = node_p->Prev)
       {
          ColPtrNode[node_p->Column] = NULL;
       }
+   }
    else
+   {
       for (int j = I[current_row], end = I[current_row+1]; j < end; j++)
       {
          ColPtrJ[J[j]] = -1;
       }
+   }
 }
 
 inline double &SparseMatrix::SearchRow(const int col)

linalg/vector.cpp

@@ -857,11 +857,11 @@ static double cuVectorMin(const int N, const double *X)
    const int bytes = min_sz*sizeof(double);
    static double *h_min = NULL;
    if (!h_min) { h_min = (double*)calloc(min_sz,sizeof(double)); }
-   static CUdeviceptr gdsr = (CUdeviceptr) NULL;
-   if (!gdsr) { ::cuMemAlloc(&gdsr,bytes); }
+   static void *gdsr = NULL;
+   if (!gdsr) { MFEM_CUDA_CHECK(cudaMalloc(&gdsr, bytes)); }
    cuKernelMin<<<gridSize,blockSize>>>(N, (double*)gdsr, x);
-   MFEM_CUDA_CHECK_RT(cudaGetLastError());
-   ::cuMemcpy((CUdeviceptr)h_min,(CUdeviceptr)gdsr,bytes);
+   MFEM_CUDA_CHECK(cudaGetLastError());
+   MFEM_CUDA_CHECK(cudaMemcpy(h_min, gdsr, bytes, cudaMemcpyDeviceToHost));
    double min = std::numeric_limits<double>::infinity();
    for (int i = 0; i < min_sz; i++) { min = fmin(min, h_min[i]); }
    return min;
@@ -909,19 +909,19 @@ static double cuVectorDot(const int N, const double *X, const double *Y)
       if (h_dot) { free(h_dot); }
       h_dot = (double*)calloc(dot_sz,sizeof(double));
    }
-   static CUdeviceptr gdsr = (CUdeviceptr) NULL;
+   static void *gdsr = NULL;
    if (!gdsr or dot_block_sz!=dot_sz)
    {
-      if (gdsr) { MFEM_CUDA_CHECK_DRV(::cuMemFree(gdsr)); }
-      MFEM_CUDA_CHECK_DRV(::cuMemAlloc(&gdsr,bytes));
+      if (gdsr) { MFEM_CUDA_CHECK(cudaFree(gdsr)); }
+      MFEM_CUDA_CHECK(cudaMalloc(&gdsr,bytes));
    }
    if (dot_block_sz!=dot_sz)
    {
       dot_block_sz = dot_sz;
    }
    cuKernelDot<<<gridSize,blockSize>>>(N, (double*)gdsr, x, y);
-   MFEM_CUDA_CHECK_RT(cudaGetLastError());
-   MFEM_CUDA_CHECK_DRV(::cuMemcpy((CUdeviceptr)h_dot,(CUdeviceptr)gdsr,bytes));
+   MFEM_CUDA_CHECK(cudaGetLastError());
+   MFEM_CUDA_CHECK(cudaMemcpy(h_dot, gdsr, bytes, cudaMemcpyDeviceToHost));
    double dot = 0.0;
    for (int i = 0; i < dot_sz; i++) { dot += h_dot[i]; }
    return dot;

makefile

@@ -211,6 +211,9 @@ ifneq ($(MFEM_USE_CUDA),YES)
    XCOMPILER = $(CXX_XCOMPILER)
    XLINKER   = $(CXX_XLINKER)
 else
+   ifneq ($(MFEM_USE_MM),YES)
+      $(error MFEM_USE_CUDA=YES requires MFEM_USE_MM=YES)
+   endif
    MFEM_CXX ?= $(CUDA_CXX)
    CXXFLAGS += $(CUDA_FLAGS) -ccbin $(CXX_OR_MPICXX)
    XCOMPILER = $(CUDA_XCOMPILER)
@@ -230,7 +233,7 @@ endif
 
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
 MFEM_REQ_LIB_DEPS = SUPERLU METIS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
- SUITESPARSE STRUMPACK GECKO GNUTLS NETCDF PETSC MPFR PUMI CUDA OCCA RAJA
+ SUITESPARSE STRUMPACK GECKO GNUTLS NETCDF PETSC MPFR PUMI OCCA RAJA
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
 
 define mfem_check_dependency
@@ -246,7 +249,7 @@ ifeq ($(MAKECMDGOALS),config)
 endif
 
 # List of MFEM dependencies, processed below
-MFEM_DEPENDENCIES = $(MFEM_REQ_LIB_DEPS) LIBUNWIND OPENMP
+MFEM_DEPENDENCIES = $(MFEM_REQ_LIB_DEPS) LIBUNWIND OPENMP CUDA
 
 # List of deprecated MFEM dependencies, processed below
 MFEM_LEGACY_DEPENDENCIES = OPENMP
@@ -319,8 +322,8 @@ MFEM_TEST_MK   ?= @MFEM_DIR@/config/test.mk
 # Use "\n" (interpreted by sed) to add a newline.
 MFEM_CONFIG_EXTRA ?= $(if $(BUILD_DIR_DEF),MFEM_BUILD_DIR ?= @MFEM_DIR@,)
 
-MFEM_SOURCE_DIR  := $(MFEM_REAL_DIR)
-MFEM_INSTALL_DIR := $(BUILD_REAL_DIR)
+MFEM_SOURCE_DIR  = $(MFEM_REAL_DIR)
+MFEM_INSTALL_DIR = $(abspath $(MFEM_PREFIX))
 
 # If we have 'config' target, export variables used by config/makefile
 ifneq (,$(filter config,$(MAKECMDGOALS)))
@@ -344,6 +347,11 @@ ifneq (,$(filter install,$(MAKECMDGOALS)))
    MFEM_LIBS     = $(if $(shared),$(INSTALL_RPATH)) -L@MFEM_LIB_DIR@ -lmfem\
       @MFEM_EXT_LIBS@
    MFEM_LIB_FILE = @MFEM_LIB_DIR@/libmfem.$(if $(shared),$(SO_VER),a)
+   ifeq ($(MFEM_USE_OCCA),YES)
+      ifneq ($(MFEM_INSTALL_DIR),$(abspath $(PREFIX))
+         $(error OCCA is enabled: PREFIX must be set during configuration!)
+      endif
+   endif
    MFEM_PREFIX := $(abspath $(PREFIX))
    MFEM_INC_DIR = $(abspath $(PREFIX_INC))
    MFEM_LIB_DIR = $(abspath $(PREFIX_LIB))
@@ -358,6 +366,7 @@ DIRS = general linalg mesh fem
 SOURCE_FILES = $(foreach dir,$(DIRS),$(wildcard $(SRC)$(dir)/*.cpp))
 RELSRC_FILES = $(patsubst $(SRC)%,%,$(SOURCE_FILES))
 OBJECT_FILES = $(patsubst $(SRC)%,$(BLD)%,$(SOURCE_FILES:.cpp=.o))
+OKL_DIRS = fem
 
 .PHONY: lib all clean distclean install config status info deps serial parallel\
 	debug pdebug cuda pcuda cudebug pcudebug style check test unittest\
@@ -498,7 +507,12 @@ install: $(if $(static),$(BLD)libmfem.a) $(if $(shared),$(BLD)libmfem.$(SO_EXT))
 # install remaining includes in each subdirectory
 	for dir in $(DIRS); do \
 	   mkdir -p $(PREFIX_INC)/mfem/$$dir && \
-	   $(INSTALL) -m 640 $(SRC)$$dir/*.hpp $(SRC)$$dir/*.okl $(PREFIX_INC)/mfem/$$dir; \
+	   $(INSTALL) -m 640 $(SRC)$$dir/*.hpp $(PREFIX_INC)/mfem/$$dir; \
+	done
+# install *.okl files
+	for dir in $(OKL_DIRS); do \
+	   mkdir -p $(PREFIX_INC)/mfem/$$dir && \
+	   $(INSTALL) -m 640 $(SRC)$$dir/*.okl $(PREFIX_INC)/mfem/$$dir; \
 	done
 # install config.mk in $(PREFIX_SHARE)
 	mkdir -p $(PREFIX_SHARE)

mesh/pmesh.cpp

@@ -3090,6 +3090,23 @@ void ParMesh::Rebalance()
                  " meshes.");
    }
 
+   // Make sure the Nodes use a ParFiniteElementSpace
+   if (Nodes && dynamic_cast<ParFiniteElementSpace*>(Nodes->FESpace()) == NULL)
+   {
+      ParFiniteElementSpace *pfes =
+         new ParFiniteElementSpace(*Nodes->FESpace(), *this);
+      ParGridFunction *new_nodes = new ParGridFunction(pfes);
+      *new_nodes = *Nodes;
+      if (Nodes->OwnFEC())
+      {
+         new_nodes->MakeOwner(Nodes->OwnFEC());
+         Nodes->MakeOwner(NULL); // takes away ownership of 'fec' and 'fes'
+         delete Nodes->FESpace();
+      }
+      delete Nodes;
+      Nodes = new_nodes;
+   }
+
    DeleteFaceNbrData();
 
    pncmesh->Rebalance();
@@ -3110,22 +3127,6 @@ void ParMesh::Rebalance()
    last_operation = Mesh::REBALANCE;
    sequence++;
 
-   // Make sure the Nodes use a ParFiniteElementSpace
-   if (Nodes && dynamic_cast<ParFiniteElementSpace*>(Nodes->FESpace()) == NULL)
-   {
-      ParFiniteElementSpace *pfes =
-         new ParFiniteElementSpace(*Nodes->FESpace(), *this);
-      ParGridFunction *new_nodes = new ParGridFunction(pfes);
-      *new_nodes = *Nodes;
-      if (Nodes->OwnFEC())
-      {
-         new_nodes->MakeOwner(Nodes->OwnFEC());
-         Nodes->MakeOwner(NULL); // takes away ownership of 'fec' and 'fes'
-         delete Nodes->FESpace();
-      }
-      delete Nodes;
-      Nodes = new_nodes;
-   }
    UpdateNodes();
 }
 

tests/unit/CMakeLists.txt

@@ -9,11 +9,11 @@
 # terms of the GNU Lesser General Public License (as published by the Free
 # Software Foundation) version 2.1 dated February 1999.
 
-# Include the build directory where mfem.hpp and mfem-performance.hpp are.
-include_directories(BEFORE ${PROJECT_BINARY_DIR})
 # Include the top mfem source directory - needed by some tests, e.g. to
 # #include "general/text.hpp".
 include_directories(BEFORE ${PROJECT_SOURCE_DIR})
+# Include the build directory where mfem.hpp and mfem-performance.hpp are.
+include_directories(BEFORE ${PROJECT_BINARY_DIR})
 # Include the source directory for the unit tests - catch.hpp is there.
 include_directories(BEFORE ${CMAKE_CURRENT_SOURCE_DIR})
 

tests/unit/fem/test_datacollection.cpp

@@ -12,7 +12,13 @@
 #include "mfem.hpp"
 #include "catch.hpp"
 #include <stdio.h>
-#include <unistd.h>  // rmdir
+
+#ifndef _WIN32
+#include <unistd.h> // rmdir
+#else
+#include <direct.h> // _rmdir
+#define rmdir(dir) _rmdir(dir)
+#endif
 
 using namespace mfem;
 

tests/unit/fem/test_inversetransform.cpp

@@ -125,7 +125,7 @@ TEST_CASE("InverseElementTransformation",
       REQUIRE( mesh_file.good() );
 
       const int npts = 100; // number of random points to test
-      const int min_found_pts = 94;
+      const int min_found_pts = 93;
       const int rand_seed = 189548;
       srand(rand_seed);