//========================================================================================
// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved.
//
// This program was produced under U.S. Government contract 89233218CNA000001 for Los
// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
// in the program are reserved by Triad National Security, LLC, and the U.S. Department
// of Energy/National Nuclear Security Administration. The Government is granted for
// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
// license in this material to reproduce, prepare derivative works, distribute copies to
// the public, perform publicly and display publicly, and to permit others to do so.
//========================================================================================

#include <string>
#include <vector>

#include "basic_types.hpp"
#include "parameter_input.hpp"
#include "staged_integrator.hpp"

namespace parthenon {

/*
 *classic Butcher Tableau integrators
 *
 * Notation:
 * c | a
 * --+--
 *   | b
 *
 * Note that a good resource is NASA-TR-R-287
 * https://ntrs.nasa.gov/citations/19680027281
 * Note that in the report, Fehlberg used the following notation,
 * as translated into standard Butcher tableau notation:
 * beta_{k lambda} = a_{ij}
 * alpha_k = c_k
 * c_k = b_k
 */
ButcherIntegrator::ButcherIntegrator(ParameterInput *pin)
    : StagedIntegrator(pin->GetOrAddString("parthenon/time", "integrator", "rk2")) {
  if (name_ == "rk1") {
    nstages = nbuffers = 1;
    Resize_(nstages);

    a[0][0] = 0;
    b[0] = 1;
    c[0] = 0;
  } else if (name_ == "rk2") {
    // Heun's method. Should match minimal storage solution
    nstages = nbuffers = 2;
    Resize_(nstages);

    a[0] = {0, 0};
    a[1] = {1, 0};
    b = {0, 1};
    c = {0, 1. / 3., 2. / 3.};
  } else if (name_ == "rk4") {
    // Classic RK4 because why not
    nstages = nbuffers = 4;
    Resize_(nstages);

    /* clang-format off */
    a[0] = {0,   0,   0, 0};
    a[1] = {0.5, 0,   0, 0};
    a[2] = {0,   0.5, 0, 0};
    a[3] = {0,   0,   1, 0};
    /* clang-format on */
    b = {1. / 6., 1. / 3., 1. / 3., 1. / 6.};
    c = {0, 0.5, 0.5, 1};
  } else if (name_ == "rk10") {
    // Feagin's family of high-order embedded methods as introduced in
    // Feagin, Neural, Parallel, and Scientific Computations 20 (2012)
    // 437-458
    // Note that Feagin uses the following notation, which we have
    // translated:
    //
    // a | beta
    // --+-----
    //   | c
    //
    // Feagin's 10th-order method has the same number of stages
    // as Fehlberg's classic RK8(9) adaptive method
    //
    // Feagin's coefficients are vailable on his website:
    // https://sce.uhcl.edu/rungekutta/

    nstages = nbuffers = 17;
    Resize_(nstages);

    // computed up to 60 digits
    c[0] = 0.000000000000000000000000000000000000000000000000000000000000;
    c[1] = 0.100000000000000000000000000000000000000000000000000000000000;
    c[2] = 0.539357840802981787532485197881302436857273449701009015505500;
    c[3] = 0.809036761204472681298727796821953655285910174551513523258250;
    c[4] = 0.309036761204472681298727796821953655285910174551513523258250;
    c[5] = 0.981074190219795268254879548310562080489056746118724882027805;
    c[6] = 0.833333333333333333333333333333333333333333333333333333333333;
    c[7] = 0.354017365856802376329264185948796742115824053807373968324184;
    c[8] = 0.882527661964732346425501486979669075182867844268052119663791;
    c[9] = 0.642615758240322548157075497020439535959501736363212695909875;
    c[10] = 0.357384241759677451842924502979560464040498263636787304090125;
    c[11] = 0.117472338035267653574498513020330924817132155731947880336209;
    c[12] = 0.833333333333333333333333333333333333333333333333333333333333;
    c[13] = 0.309036761204472681298727796821953655285910174551513523258250;
    c[14] = 0.539357840802981787532485197881302436857273449701009015505500;
    c[15] = 0.100000000000000000000000000000000000000000000000000000000000;
    c[16] = 1.00000000000000000000000000000000000000000000000000000000000;

    b[0] = 0.0333333333333333333333333333333333333333333333333333333333333;
    b[1] = 0.0250000000000000000000000000000000000000000000000000000000000;
    b[2] = 0.0333333333333333333333333333333333333333333333333333333333333;
    b[3] = 0.000000000000000000000000000000000000000000000000000000000000;
    b[4] = 0.0500000000000000000000000000000000000000000000000000000000000;
    b[5] = 0.000000000000000000000000000000000000000000000000000000000000;
    b[6] = 0.0400000000000000000000000000000000000000000000000000000000000;
    b[7] = 0.000000000000000000000000000000000000000000000000000000000000;
    b[8] = 0.189237478148923490158306404106012326238162346948625830327194;
    b[9] = 0.277429188517743176508360262560654340428504319718040836339472;
    b[10] = 0.277429188517743176508360262560654340428504319718040836339472;
    b[11] = 0.189237478148923490158306404106012326238162346948625830327194;
    b[12] = -0.0400000000000000000000000000000000000000000000000000000000000;
    b[13] = -0.0500000000000000000000000000000000000000000000000000000000000;
    b[14] = -0.0333333333333333333333333333333333333333333333333333333333333;
    b[15] = -0.0250000000000000000000000000000000000000000000000000000000000;
    b[16] = 0.0333333333333333333333333333333333333333333333333333333333333;

    // Set matrix to zero first, to only set non-zero coeffs
    for (int i = 0; i < nstages; ++i) {
      for (int j = 0; j < nstages; ++j) {
        a[i][j] = 0;
      }
    }
    a[1][0] = 0.100000000000000000000000000000000000000000000000000000000000;
    a[2][0] = -0.915176561375291440520015019275342154318951387664369720564660;
    a[2][1] = 1.45453440217827322805250021715664459117622483736537873607016;
    a[3][0] = 0.202259190301118170324681949205488413821477543637878380814562;
    a[3][2] = 0.606777570903354510974045847616465241464432630913635142443687;
    a[4][0] = 0.184024714708643575149100693471120664216774047979591417844635;
    a[4][2] = 0.197966831227192369068141770510388793370637287463360401555746;
    a[4][3] = -0.0729547847313632629185146671595558023015011608914382961421311;
    a[5][0] = 0.0879007340206681337319777094132125475918886824944548534041378;
    a[5][3] = 0.410459702520260645318174895920453426088035325902848695210406;
    a[5][4] = 0.482713753678866489204726942976896106809132737721421333413261;
    a[6][0] = 0.0859700504902460302188480225945808401411132615636600222593880;
    a[6][3] = 0.330885963040722183948884057658753173648240154838402033448632;
    a[6][4] = 0.489662957309450192844507011135898201178015478433790097210790;
    a[6][5] = -0.0731856375070850736789057580558988816340355615025188195854775;
    a[7][0] = 0.120930449125333720660378854927668953958938996999703678812621;
    a[7][4] = 0.260124675758295622809007617838335174368108756484693361887839;
    a[7][5] = 0.0325402621549091330158899334391231259332716675992700000776101;
    a[7][6] = -0.0595780211817361001560122202563305121444953672762930724538856;
    a[8][0] = 0.110854379580391483508936171010218441909425780168656559807038;
    a[8][5] = -0.0605761488255005587620924953655516875526344415354339234619466;
    a[8][6] = 0.321763705601778390100898799049878904081404368603077129251110;
    a[8][7] = 0.510485725608063031577759012285123416744672137031752354067590;
    a[9][0] = 0.112054414752879004829715002761802363003717611158172229329393;
    a[9][5] = -0.144942775902865915672349828340980777181668499748506838876185;
    a[9][6] = -0.333269719096256706589705211415746871709467423992115497968724;
    a[9][7] = 0.499269229556880061353316843969978567860276816592673201240332;
    a[9][8] = 0.509504608929686104236098690045386253986643232352989602185060;
    a[10][0] = 0.113976783964185986138004186736901163890724752541486831640341;
    a[10][5] = -0.0768813364203356938586214289120895270821349023390922987406384;
    a[10][6] = 0.239527360324390649107711455271882373019741311201004119339563;
    a[10][7] = 0.397774662368094639047830462488952104564716416343454639902613;
    a[10][8] = 0.0107558956873607455550609147441477450257136782823280838547024;
    a[10][9] = -0.327769124164018874147061087350233395378262992392394071906457;
    a[11][0] = 0.0798314528280196046351426864486400322758737630423413945356284;
    a[11][5] = -0.0520329686800603076514949887612959068721311443881683526937298;
    a[11][6] = -0.0576954146168548881732784355283433509066159287152968723021864;
    a[11][7] = 0.194781915712104164976306262147382871156142921354409364738090;
    a[11][8] = 0.145384923188325069727524825977071194859203467568236523866582;
    a[11][9] = -0.0782942710351670777553986729725692447252077047239160551335016;
    a[11][10] = -0.114503299361098912184303164290554670970133218405658122674674;
    a[12][0] = 0.985115610164857280120041500306517278413646677314195559520529;
    a[12][3] = 0.330885963040722183948884057658753173648240154838402033448632;
    a[12][4] = 0.489662957309450192844507011135898201178015478433790097210790;
    a[12][5] = -1.37896486574843567582112720930751902353904327148559471526397;
    a[12][6] = -0.861164195027635666673916999665534573351026060987427093314412;
    a[12][7] = 5.78428813637537220022999785486578436006872789689499172601856;
    a[12][8] = 3.28807761985103566890460615937314805477268252903342356581925;
    a[12][9] = -2.38633905093136384013422325215527866148401465975954104585807;
    a[12][10] = -3.25479342483643918654589367587788726747711504674780680269911;
    a[12][11] = -2.16343541686422982353954211300054820889678036420109999154887;
    a[13][0] = 0.895080295771632891049613132336585138148156279241561345991710;
    a[13][2] = 0.197966831227192369068141770510388793370637287463360401555746;
    a[13][3] = -0.0729547847313632629185146671595558023015011608914382961421311;
    a[13][5] = -0.851236239662007619739049371445966793289359722875702227166105;
    a[13][6] = 0.398320112318533301719718614174373643336480918103773904231856;
    a[13][7] = 3.63937263181035606029412920047090044132027387893977804176229;
    a[13][8] = 1.54822877039830322365301663075174564919981736348973496313065;
    a[13][9] = -2.12221714704053716026062427460427261025318461146260124401561;
    a[13][10] = -1.58350398545326172713384349625753212757269188934434237975291;
    a[13][11] = -1.71561608285936264922031819751349098912615880827551992973034;
    a[13][12] = -0.0244036405750127452135415444412216875465593598370910566069132;
    a[14][0] = -0.915176561375291440520015019275342154318951387664369720564660;
    a[14][1] = 1.45453440217827322805250021715664459117622483736537873607016;
    a[14][4] = -0.777333643644968233538931228575302137803351053629547286334469;
    a[14][6] = -0.0910895662155176069593203555807484200111889091770101799647985;
    a[14][12] = 0.0910895662155176069593203555807484200111889091770101799647985;
    a[14][13] = 0.777333643644968233538931228575302137803351053629547286334469;
    a[15][0] = 0.100000000000000000000000000000000000000000000000000000000000;
    a[15][2] = -0.157178665799771163367058998273128921867183754126709419409654;
    a[15][14] = 0.157178665799771163367058998273128921867183754126709419409654;
    a[16][0] = 0.181781300700095283888472062582262379650443831463199521664945;
    a[16][1] = 0.675000000000000000000000000000000000000000000000000000000000;
    a[16][2] = 0.342758159847189839942220553413850871742338734703958919937260;
    a[16][4] = 0.259111214548322744512977076191767379267783684543182428778156;
    a[16][5] = -0.358278966717952089048961276721979397739750634673268802484271;
    a[16][6] = -1.04594895940883306095050068756409905131588123172378489286080;
    a[16][7] = 0.930327845415626983292300564432428777137601651182965794680397;
    a[16][8] = 1.77950959431708102446142106794824453926275743243327790536000;
    a[16][9] = 0.100000000000000000000000000000000000000000000000000000000000;
    a[16][10] = -0.282547569539044081612477785222287276408489375976211189952877;
    a[16][11] = -0.159327350119972549169261984373485859278031542127551931461821;
    a[16][12] = -0.145515894647001510860991961081084111308650130578626404945571;
    a[16][13] = -0.259111214548322744512977076191767379267783684543182428778156;
    a[16][14] = -0.342758159847189839942220553413850871742338734703958919937260;
    a[16][15] = -0.675000000000000000000000000000000000000000000000000000000000;
  }
}

void ButcherIntegrator::Resize_(int nstages) {
  a.resize(nstages);
  for (int i = 0; i < a.size(); ++i) {
    a[i].resize(nstages);
  }
  b.resize(nstages);
  c.resize(nstages);
}

} // namespace parthenon
