diff --git a/CHANGELOG.md b/CHANGELOG.md
index b93363159e001edb7469d04b8ba4d4e08e2358b6..7fbda3edc2687139b517ac7647058cc0a4d2ef4a 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 e0f741825334c1889b302f8bf0c06a56460e4dd7..ec2ccb1ff9951b714bf49a0263ebd7455949f2d8 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 83326ea14c46c8d94efd8eb20f3f6b07b07792d7..7b2a7427584f46c49424151b3ad27bf5cdfd288f 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 31d05ba870145eafe1325c004d2b650ad8ce82c3..5e0afe134ecc1bb5f7177afb05cdd5c43889a94f 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 79f21c76f0ecea387dd4a161e09917dc33ce239a..d92d62cfac563d2490736f22b7e2cfd4a9990406 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 ba3adadbae7cf23ed89eb5c728f71d44f45e0547..cafddd60352f90c5cc66045e6f06f45c8bbf36b9 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 93f1d0be3535780ff9b65888dcf95df0f866c9e9..825c1344958ce2e0ec839cbc9a888d0cf5f10a88 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 4f7fb7920dcdf0430793dfd95283e7b03ab34e46..043f3947b3163928d318b0f71b89df30acbec5d3 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 3f2bf47cdad2277875999c1c92a70a03bc21052e..1398466223aeb6715dc88ccfe5c34c6cf2c463c6 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 132f4fc8f6e2bc4a04aa46b8d5ff505054ed5322..40cf90a6831293f53b6a79fc68084133ab3c9df1 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 9718311a036d2a1174739851c1efe70c1b6ff4a6..493186161cef788b3f4c56c130878567803c79d0 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 6e84dd716ed0b3c927e39a1474b2552bf15b7023..315e7f8ff8ec8e9536f4d2daef38d7dfcadba7ab 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 aa9bab3ff667e823cecc5311f0b61e5d7656f77f..0000000000000000000000000000000000000000
--- 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 33a82e3d9acb509774b2d0f9af25a17f767ab7ea..0000000000000000000000000000000000000000
--- 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 bf529307af04b93d4bd81229ce4e50d4097c391b..0000000000000000000000000000000000000000
--- 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 ef5fa158d1a0c2c9063d616fd763051ca273d783..25729c164131f838f9cc62aed6cea1de76715089 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 de5b5bff27f805dbfeb945bf95e2feb8eb37ed95..b6233d8153dc8061c7e0d7e63b1a45147ebc0f0d 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 927478f16621f3adf9e79687e62681b47cbc9919..3d6e936774b2ba492ecb1ff2509f655d4d09dd65 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 8ad79add17a8ef1cabe50920cc5485b1c8a645e4..6f5fb36bbde3ba49884444a196d5435f889c8335 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;