From 3ed1176ca74646d3f6b678e8872d62c2ba588bd1 Mon Sep 17 00:00:00 2001
From: Peter Hill <zed.three@gmail.com>
Date: Tue, 23 Feb 2021 11:04:59 +0000
Subject: [PATCH] Remove the Karniadakis time solver

---
 CHANGELOG.md                                  |   6 +-
 CMakeLists.txt                                |   2 -
 .../IMEX/advection-reaction/data/BOUT.inp     |   2 +-
 .../constraints/laplace-dae/data/BOUT.inp     |   2 +-
 examples/eigen-box/data/BOUT.inp              |   2 +-
 examples/elm-pb/eigen/BOUT.inp                |   2 +-
 examples/gravity_reduced/data/BOUT.inp        |   2 +-
 examples/gyro-gem/data/BOUT.inp               |   2 +-
 examples/hasegawa-wakatani/data/BOUT.inp      |   2 +-
 examples/jorek-compare/data/BOUT.inp          |   2 +-
 include/bout/solver.hxx                       |   1 -
 manual/sphinx/user_docs/time_integration.rst  |   4 +-
 src/solver/impls/karniadakis/karniadakis.cxx  | 196 ------------------
 src/solver/impls/karniadakis/karniadakis.hxx  |  79 -------
 src/solver/impls/karniadakis/makefile         |   8 -
 src/solver/impls/makefile                     |   2 +-
 src/solver/solver.cxx                         |   1 -
 tests/integrated/test-delp2/data/BOUT.inp     |   2 +-
 tests/integrated/test-solver/test_solver.cxx  |   3 -
 19 files changed, 16 insertions(+), 304 deletions(-)
 delete mode 100644 src/solver/impls/karniadakis/karniadakis.cxx
 delete mode 100644 src/solver/impls/karniadakis/karniadakis.hxx
 delete mode 100644 src/solver/impls/karniadakis/makefile

diff --git a/CHANGELOG.md b/CHANGELOG.md
index b93363159..7fbda3edc 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -47,7 +47,11 @@
   have been set; this is only a problem if you don't call
   `BoutInitialise`. [\#2062](https://github.com/boutproject/BOUT-dev/pull/2062)
 - Support for reading/writing HDF5 files has been removed ahead of completely
-  refactoring the I/O systems. [\#2208](https://github.com/boutproject/BOUT-dev/pull/2208)
+  refactoring the I/O
+  systems. [\#2208](https://github.com/boutproject/BOUT-dev/pull/2208)
+- Removed the Karniadakis time solver. Other choices for split-operator schemes
+  are: `splitrk` (built-in), `imexbdf2` (requires PETSc), and `arkode` (requires
+  SUNDIALS) [\#2241](https://github.com/boutproject/BOUT-dev/pull/2241)
 
 
 ## [v4.3.2](https://github.com/boutproject/BOUT-dev/tree/v4.3.2) (2020-10-19)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index e0f741825..ec2ccb1ff 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -251,8 +251,6 @@ set(BOUT_SOURCES
   ./src/solver/impls/ida/ida.hxx
   ./src/solver/impls/imex-bdf2/imex-bdf2.cxx
   ./src/solver/impls/imex-bdf2/imex-bdf2.hxx
-  ./src/solver/impls/karniadakis/karniadakis.cxx
-  ./src/solver/impls/karniadakis/karniadakis.hxx
   ./src/solver/impls/petsc/petsc.cxx
   ./src/solver/impls/petsc/petsc.hxx
   ./src/solver/impls/power/power.cxx
diff --git a/examples/IMEX/advection-reaction/data/BOUT.inp b/examples/IMEX/advection-reaction/data/BOUT.inp
index 83326ea14..7b2a74275 100644
--- a/examples/IMEX/advection-reaction/data/BOUT.inp
+++ b/examples/IMEX/advection-reaction/data/BOUT.inp
@@ -58,7 +58,7 @@ ATOL = 1e-10 # absolute tolerance
 RTOL = 1e-05  # relative tolerance
 mxstep = 50000
 
-type = karniadakis
+type = splitrk
 timestep = 0.0001
 
 ##################################################
diff --git a/examples/constraints/laplace-dae/data/BOUT.inp b/examples/constraints/laplace-dae/data/BOUT.inp
index 31d05ba87..5e0afe134 100644
--- a/examples/constraints/laplace-dae/data/BOUT.inp
+++ b/examples/constraints/laplace-dae/data/BOUT.inp
@@ -59,7 +59,7 @@ ATOL = 1e-10 # absolute tolerance
 RTOL = 1e-05  # relative tolerance
 mxstep = 50000
 
-#type=karniadakis
+#type=splitrk
 #timestep=1e-4
 use_precon = true
 
diff --git a/examples/eigen-box/data/BOUT.inp b/examples/eigen-box/data/BOUT.inp
index 79f21c76f..d92d62cfa 100644
--- a/examples/eigen-box/data/BOUT.inp
+++ b/examples/eigen-box/data/BOUT.inp
@@ -37,7 +37,7 @@ selfSolve = false
 timestep = 0.01
 
 [solver:advance]
-type = karniadakis
+type = splitrk
 timestep = 0.001
 
 [f]
diff --git a/examples/elm-pb/eigen/BOUT.inp b/examples/elm-pb/eigen/BOUT.inp
index ba3adadba..cafddd603 100644
--- a/examples/elm-pb/eigen/BOUT.inp
+++ b/examples/elm-pb/eigen/BOUT.inp
@@ -79,7 +79,7 @@ tol = 0.0001
 selfSolve = false
 
 [solver:advance]
-type = karniadakis
+type = splitrk
 timestep = 0.01
 
 
diff --git a/examples/gravity_reduced/data/BOUT.inp b/examples/gravity_reduced/data/BOUT.inp
index 93f1d0be3..825c13449 100644
--- a/examples/gravity_reduced/data/BOUT.inp
+++ b/examples/gravity_reduced/data/BOUT.inp
@@ -69,7 +69,7 @@ upwind = W3
 ATOL = 1e-10 # absolute tolerance
 RTOL = 1e-05  # relative tolerance
 
-#type=karniadakis
+#type=splitrk
 #timestep = 1e-3
 
 ##################################################
diff --git a/examples/gyro-gem/data/BOUT.inp b/examples/gyro-gem/data/BOUT.inp
index 4f7fb7920..043f3947b 100644
--- a/examples/gyro-gem/data/BOUT.inp
+++ b/examples/gyro-gem/data/BOUT.inp
@@ -59,7 +59,7 @@ upwind = W3
 ATOL = 1e-12 # absolute tolerance
 RTOL = 1e-05  # relative tolerance
 mxstep = 50000
-type = karniadakis
+type = splitrk
 timestep = 1e-06
 
 ##################################################
diff --git a/examples/hasegawa-wakatani/data/BOUT.inp b/examples/hasegawa-wakatani/data/BOUT.inp
index 3f2bf47cd..139846622 100644
--- a/examples/hasegawa-wakatani/data/BOUT.inp
+++ b/examples/hasegawa-wakatani/data/BOUT.inp
@@ -22,7 +22,7 @@ dz = 0.2
 flags = 0   # Flags for Laplacian inversion
 
 [solver]
-#type=karniadakis
+#type=splitrk
 timestep = 0.001
 
 [hw]
diff --git a/examples/jorek-compare/data/BOUT.inp b/examples/jorek-compare/data/BOUT.inp
index 132f4fc8f..40cf90a68 100644
--- a/examples/jorek-compare/data/BOUT.inp
+++ b/examples/jorek-compare/data/BOUT.inp
@@ -71,7 +71,7 @@ ATOL = 1e-10 # absolute tolerance
 RTOL = 1e-05  # relative tolerance
 mxstep = 50000
 
-#type=karniadakis
+#type=splitrk
 #timestep=5e-5
 
 ##################################################
diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx
index 9718311a0..493186161 100644
--- a/include/bout/solver.hxx
+++ b/include/bout/solver.hxx
@@ -84,7 +84,6 @@ constexpr auto SOLVERPVODE = "pvode";
 constexpr auto SOLVERIDA = "ida";
 constexpr auto SOLVERPETSC = "petsc";
 constexpr auto SOLVERSLEPC = "slepc";
-constexpr auto SOLVERKARNIADAKIS = "karniadakis";
 constexpr auto SOLVERRK4 = "rk4";
 constexpr auto SOLVEREULER = "euler";
 constexpr auto SOLVERRK3SSP = "rk3ssp";
diff --git a/manual/sphinx/user_docs/time_integration.rst b/manual/sphinx/user_docs/time_integration.rst
index 6e84dd716..315e7f8ff 100644
--- a/manual/sphinx/user_docs/time_integration.rst
+++ b/manual/sphinx/user_docs/time_integration.rst
@@ -43,8 +43,6 @@ needed to make the solver available.
    +---------------+-----------------------------------------+--------------------+
    | rkgeneric     | Generic Runge Kutta explicit methods    | Always available   |
    +---------------+-----------------------------------------+--------------------+
-   | karniadakis   | Karniadakis explicit method             | Always available   |
-   +---------------+-----------------------------------------+--------------------+
    | rk3ssp        | 3rd-order Strong Stability Preserving   | Always available   |
    +---------------+-----------------------------------------+--------------------+
    | splitrk       | Split RK3-SSP and RK-Legendre           | Always available   |
@@ -83,7 +81,7 @@ given in table :numref:`tab-solveropts`.
    +------------------+--------------------------------------------+-------------------------------------+
    | max\_timestep    | Maximum timestep                           | rk4, cvode                          |
    +------------------+--------------------------------------------+-------------------------------------+
-   | timestep         | Starting timestep                          | rk4, karniadakis, euler, imexbdf2   |
+   | timestep         | Starting timestep                          | rk4, euler, imexbdf2                |
    +------------------+--------------------------------------------+-------------------------------------+
    | adaptive         | Adapt timestep? (Y/N)                      | rk4, imexbdf2                       |
    +------------------+--------------------------------------------+-------------------------------------+
diff --git a/src/solver/impls/karniadakis/karniadakis.cxx b/src/solver/impls/karniadakis/karniadakis.cxx
deleted file mode 100644
index aa9bab3ff..000000000
--- a/src/solver/impls/karniadakis/karniadakis.cxx
+++ /dev/null
@@ -1,196 +0,0 @@
-/**************************************************************************
- * Karniadakis split-operator solver
- * 
- * Formulation from:
- * "GEM - An Energy Conserving Electromagnetic Gyrofluid Model"
- *  by Bruce D Scott. arXiv:physics/0501124v1 23 Jan 2005 
- *
- * Original paper:
- *   J. Comput. Phys. 97 (1991) p414-443
- * 
- * Always available, since doesn't depend on external library
- * 
- * Solves a system df/dt = S(f) + D(f)
- * 
- * where S is the RHS of each equation, and D is the diffusion terms
- * 
- **************************************************************************
- * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu
- *
- * Contact: Ben Dudson, bd512@york.ac.uk
- * 
- * This file is part of BOUT++.
- *
- * BOUT++ 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, either version 3 of the License, or
- * (at your option) any later version.
- *
- * BOUT++ is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with BOUT++.  If not, see <http://www.gnu.org/licenses/>.
- *
- **************************************************************************/
-
-#include "karniadakis.hxx"
-
-#include <boutcomm.hxx>
-#include <msg_stack.hxx>
-#include <output.hxx>
-#include <bout/openmpwrap.hxx>
-
-KarniadakisSolver::KarniadakisSolver(Options *options) : Solver(options) {
-  canReset = true;  
-}
-
-int KarniadakisSolver::init(int nout, BoutReal tstep) {
-  TRACE("Initialising Karniadakis solver");
-  
-  output_error << "\nWARNING:\n"
-    "        The Karniadakis solver is now deprecated and will be removed in BOUT++ 5.0!\n"
-    "        Try the \"splitrk\", \"imexbdf2\" (requires PETSc) or \"arkode\" (requires SUNDIALS)\n"
-    "        solvers for other split-schemes\n\n";
-
-  /// Call the generic initialisation first
-  if (Solver::init(nout, tstep))
-    return 1;
-  
-  output << "\n\tKarniadakis solver\n";
- 
-  nsteps = nout; // Save number of output steps
-  out_timestep = tstep;
-  
-  // Calculate number of variables
-  nlocal = getLocalN();
-  
-  // Get total problem size
-  int neq;
-  if (bout::globals::mpi->MPI_Allreduce(&nlocal, &neq, 1, MPI_INT, MPI_SUM,
-                                        BoutComm::get())) {
-    output_error.write("\tERROR: MPI_Allreduce failed!\n");
-    return 1;
-  }
-  
-  output.write("\t3d fields = {:d}, 2d fields = {:d} neq={:d}, local_N={:d}\n",
-	       n3Dvars(), n2Dvars(), neq, nlocal);
-  
-  // Allocate memory
-
-  f1.reallocate(nlocal);
-  f0.reallocate(nlocal);
-  fm1.reallocate(nlocal);
-  fm2.reallocate(nlocal);
-
-  S0.reallocate(nlocal);
-  Sm1.reallocate(nlocal);
-  Sm2.reallocate(nlocal);
-
-  D0.reallocate(nlocal);
-
-  first_time = true;
-
-  // Put starting values into f0
-  save_vars(std::begin(f0));
-
-  // Get options
-  OPTION(options, timestep, tstep);
-  
-  // Make sure timestep divides into tstep
-  
-  // Number of sub-steps, rounded up
-  nsubsteps = static_cast<int>(std::round(tstep / timestep));
-
-  output.write("\tNumber of substeps: {:e} / {:e} -> {:d}\n", tstep, timestep, nsubsteps);
-
-  timestep = tstep / static_cast<BoutReal>(nsubsteps);
-
-  return 0;
-}
-
-int KarniadakisSolver::run() {
-  TRACE("KarniadakisSolver::run()");
-  
-  for(int i=0;i<nsteps;i++) {
-    // Run through a fixed number of steps
-    for(int j=0; j<nsubsteps; j++) {
-      // Advance f0 -> f1
-      take_step(timestep);
-      
-      // Cycle buffers
-      auto tmp = fm2;
-      fm2 = fm1;
-      fm1 = f0;
-      f0 = f1;
-      f1 = tmp;
-      
-      tmp = Sm2;
-      Sm2 = Sm1;
-      Sm1 = S0;
-      S0 = tmp;
-      
-      simtime += timestep;
-      
-      call_timestep_monitors(simtime, timestep);
-    }
-    iteration++;
-    
-    // Call RHS to communicate and get auxilliary variables
-    load_vars(std::begin(f0));
-    run_rhs(simtime);
-    
-    /// Call the monitor function
-    
-    if(call_monitors(simtime, i, nsteps)) {
-      // User signalled to quit
-      break;
-    }
-  }
-  
-  return 0;
-}
-
-void KarniadakisSolver::resetInternalFields(){
-  //Make sure all other fields get reset
-  first_time=true;
-
-  //Copy fields into vector
-  save_vars(std::begin(f0));
-}
-
-void KarniadakisSolver::take_step(BoutReal dt) {
-  // S0 = S(f0)
-
-  load_vars(std::begin(f0));
-  run_convective(simtime);
-  save_derivs(std::begin(S0));
-
-  if(first_time) {
-    // Initialise values
-    BOUT_OMP(parallel for)
-    for(int i=0;i<nlocal;i++) {
-    //fm1[i] = fm2[i] = f0[i];
-      fm1[i] = f0[i] - dt*S0[i];
-      fm2[i] = fm1[i] - dt*S0[i];
-      Sm1[i] = Sm2[i] = S0[i];
-    }
-    first_time = false;
-  }
-
-  BOUT_OMP(parallel for)
-  for(int i=0;i<nlocal;i++)
-    f1[i] = (6./11.) * (3.*f0[i] - 1.5*fm1[i] + (1./3.)*fm2[i] + dt*(3.*S0[i] - 3.*Sm1[i] + Sm2[i]));
-  
-  // D0 = S(f0)
-  load_vars(std::begin(f0));
-  run_diffusive(simtime);
-  save_derivs(std::begin(D0));
-
-  // f1 = f1 + dt*D0
-  BOUT_OMP(parallel for)
-  for(int i=0;i<nlocal;i++)
-    f1[i] += (6./11.) * dt*D0[i];
-}
diff --git a/src/solver/impls/karniadakis/karniadakis.hxx b/src/solver/impls/karniadakis/karniadakis.hxx
deleted file mode 100644
index 33a82e3d9..000000000
--- a/src/solver/impls/karniadakis/karniadakis.hxx
+++ /dev/null
@@ -1,79 +0,0 @@
-/**************************************************************************
- * Karniadakis split-operator solver
- *   J. Comput. Phys. 97 (1991) p414-443
- * 
- * Always available, since doesn't depend on external library
- * 
- * Solves a system df/dt = S(f) + D(f)
- * 
- * where S is the RHS of each equation, and D is the diffusion terms
- * 
- **************************************************************************
- * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu
- *
- * Contact: Ben Dudson, bd512@york.ac.uk
- * 
- * This file is part of BOUT++.
- *
- * BOUT++ 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, either version 3 of the License, or
- * (at your option) any later version.
- *
- * BOUT++ is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with BOUT++.  If not, see <http://www.gnu.org/licenses/>.
- *
- **************************************************************************/
-
-class KarniadakisSolver;
-
-#ifndef __KARNIADAKIS_SOLVER_H__
-#define __KARNIADAKIS_SOLVER_H__
-
-#include "mpi.h"
-
-#include <bout_types.hxx>
-#include <bout/solver.hxx>
-
-namespace {
-RegisterSolver<KarniadakisSolver> registersolverkarniadakis("karniadakis");
-}
-
-class KarniadakisSolver : public Solver {
- public:
-  KarniadakisSolver(Options *options);
-  ~KarniadakisSolver(){};
-
-  BoutReal getCurrentTimestep() override {return timestep; }
-
-  int init(int nout, BoutReal tstep) override;
-  
-  int run() override;
-  void resetInternalFields() override;
-
- private:
-  
-  Array<BoutReal> f1, f0, fm1, fm2; // System state at current, and two previous time points
-  Array<BoutReal> S0, Sm1, Sm2; // Convective part of the RHS equations
-  Array<BoutReal> D0;             // Dissipative part of the RHS
-  
-  bool first_time; // Need to initialise values
-
-  BoutReal out_timestep; // The output timestep
-  int nsteps; // Number of output steps
-  
-  BoutReal timestep; // The internal timestep
-  int nsubsteps; // Number of sub steps
-  
-  int nlocal; // Number of variables on local processor
-  
-  void take_step(BoutReal dt); // Take a single step to calculate f1
-};
-
-#endif // __KARNIADAKIS_SOLVER_H__
-
diff --git a/src/solver/impls/karniadakis/makefile b/src/solver/impls/karniadakis/makefile
deleted file mode 100644
index bf529307a..000000000
--- a/src/solver/impls/karniadakis/makefile
+++ /dev/null
@@ -1,8 +0,0 @@
-
-BOUT_TOP = ../../../..
-
-SOURCEC		= karniadakis.cxx
-SOURCEH		= $(SOURCEC:%.cxx=%.hxx)
-TARGET		= lib
-
-include $(BOUT_TOP)/make.config
diff --git a/src/solver/impls/makefile b/src/solver/impls/makefile
index ef5fa158d..25729c164 100644
--- a/src/solver/impls/makefile
+++ b/src/solver/impls/makefile
@@ -6,7 +6,7 @@ DIRS		= arkode \
 	petsc \
 	snes imex-bdf2 \
 	power slepc adams_bashforth \
-	karniadakis rk4 euler rk3-ssp rkgeneric split-rk
+	rk4 euler rk3-ssp rkgeneric split-rk
 
 TARGET		= lib
 
diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx
index de5b5bff2..b6233d815 100644
--- a/src/solver/solver.cxx
+++ b/src/solver/solver.cxx
@@ -48,7 +48,6 @@
 #include "impls/euler/euler.hxx"
 #include "impls/ida/ida.hxx"
 #include "impls/imex-bdf2/imex-bdf2.hxx"
-#include "impls/karniadakis/karniadakis.hxx"
 #include "impls/petsc/petsc.hxx"
 #include "impls/power/power.hxx"
 #include "impls/pvode/pvode.hxx"
diff --git a/tests/integrated/test-delp2/data/BOUT.inp b/tests/integrated/test-delp2/data/BOUT.inp
index 927478f16..3d6e93677 100644
--- a/tests/integrated/test-delp2/data/BOUT.inp
+++ b/tests/integrated/test-delp2/data/BOUT.inp
@@ -17,7 +17,7 @@ dx = 0.2
 dy = 1.0
 
 [solver]
-#type=karniadakis
+#type=splitrk
 timestep = 0.001
 
 [diffusion]
diff --git a/tests/integrated/test-solver/test_solver.cxx b/tests/integrated/test-solver/test_solver.cxx
index 8ad79add1..6f5fb36bb 100644
--- a/tests/integrated/test-solver/test_solver.cxx
+++ b/tests/integrated/test-solver/test_solver.cxx
@@ -101,9 +101,6 @@ int main(int argc, char** argv) {
   root["imexbdf2"]["adaptive"] = true;
   root["imexbdf2"]["adaptRtol"] = 1.e-5;
 
-  root["karniadakis"]["nout"] = 100;
-  root["karniadakis"]["timestep"] = end / (NOUT * 10000);
-
   root["petsc"]["nout"] = 10000;
   root["petsc"]["output_step"] = end / 10000;
 
-- 
GitLab