/***********************************************************************************
                 OpenCS Project: www.daetools.com
                 Copyright (C) Dragan Nikolic
************************************************************************************
OpenCS is free software; you can redistribute it and/or modify it under the
terms of the GNU Lesser General Public License version 3 as published by the Free Software
Foundation. OpenCS 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 the
OpenCS software; if not, see <http://www.gnu.org/licenses/>.
***********************************************************************************/
#include <mpi.h>
#include "brusselator_kernels_2d.h"
#include <OpenCS/opencs.h>
#include <OpenCS/models/partitioner_metis.h>
using namespace cs;

/* Reimplementation of IDAS idasBruss_kry_bbd_p example.
 * The PDE system is a two-species time-dependent PDE known as
 * Brusselator PDE and models a chemically reacting system:
 *
 *  du/dt = eps1(d2u/dx2  + d2u/dy2) + u^2 v - (B+1)u + A
 *  dv/dt = eps2(d2v/dx2  + d2v/dy2) - u^2 v + Bu
 *
 *  BC: Homogenous Neumann
 *  IC: u(x,y,t0) = u0(x,y) =  1  - 0.5*cos(pi*y/L)
 *      v(x,y,t0) = v0(x,y) = 3.5 - 2.5*cos(pi*x/L)
 *
 * The PDEs are discretized by central differencing on a uniform (Nx, Ny) grid.
 * The original results are in dae_example_3.csv file.
 *
 * The model is described in:
 *  - R. Serban and A. C. Hindmarsh. CVODES, the sensitivity-enabled ODE solver in SUNDIALS.
 *    In Proceedings of the 5th International Conference on Multibody Systems,
 *    Nonlinear Dynamics and Control, Long Beach, CA, 2005. ASME.
 *  - M. R. Wittman. Testing of PVODE, a Parallel ODE Solver.
 *    Technical Report UCRL-ID-125562, LLNL, August 1996. */
int main(int argc, char *argv[])
{
    uint32_t Nx = 80;
    uint32_t Ny = 80;

    if(argc == 3)
    {
        Nx = atoi(argv[1]);
        Ny = atoi(argv[2]);
    }
    else if(argc == 1)
    {
        // Use the default grid size.
    }
    else
    {
        printf("Usage:\n");
        printf("  %s (using the default grid: %u x %u)\n", argv[0], Nx, Ny);
        printf("  %s Nx Ny\n", argv[0]);
        return -1;
    }

    printf("########################################################\n");
    printf(" Brusselator kernels model in 2D: %u x %u grid\n", Nx, Ny);
    printf("########################################################\n");

    // Homogenous Neumann BCs at all four edges.
    csNumber_t bc_u_flux(0.0);
    csNumber_t bc_v_flux(0.0);
    BrusselatorKernels_2D bruss(Nx, Ny, bc_u_flux, bc_v_flux);

    uint32_t    Ndofs                         = 0;
    uint32_t    Nvariables                    = bruss.Nequations;
    real_t      defaultVariableValue          = 0.0;
    real_t      defaultVariableTimeDerivative = 0.0;
    real_t      defaultAbsoluteTolerance      = 1e-5;
    std::string defaultVariableName           = "x";

    /* 1. Initialise model builder with the number of variables/equations. */
    csModelBuilder_t mb;
    mb.Initialize_DAE_System(Nvariables,
                             Ndofs,
                             defaultVariableValue,
                             defaultVariableTimeDerivative,
                             defaultAbsoluteTolerance,
                             defaultVariableName);
    printf("Model builder initialised\n");

    /* 2. Create and set model equations. */
    csNumber_t                     time    = mb.GetTime();
    const std::vector<csNumber_t>& x_vars  = mb.GetVariables();
    const std::vector<csNumber_t>& xt_vars = mb.GetTimeDerivatives();

    std::vector<csEquation_t> equations;
    std::vector<csKernelPtr>  kernels;
    bruss.CreateEquations(x_vars, xt_vars, time, equations, kernels);
    printf("Model equations generated\n");

    mb.SetModelEquations(equations, kernels);
    printf("Model equations set\n");

    // Set variable names
    std::vector<std::string> names;
    bruss.GetVariableNames(names);
    mb.SetVariableNames(names);

    /* 3. Export model(s) to a specified directory for sequential and parallel simulations. */
    // Set initial conditions.
    std::vector<real_t> uv0(Nvariables, 0.0);
    bruss.SetInitialConditions(uv0);
    mb.SetVariableValues(uv0);

    // Set the simulation options.
    csSimulationOptionsPtr options = mb.GetSimulationOptions();
    options->LoadString(simulation_options);
    options->SetDouble("Simulation.TimeHorizon",              20.0);
    options->SetDouble("Simulation.ReportingInterval",         0.5);
    options->SetDouble("Solver.Parameters.RelativeTolerance", 1e-5);

    std::vector<std::string> grps  = {"BoundaryConditions"};
    std::vector<std::string> krnls = {"Brusselator_u", "Brusselator_v"};
    /* (1) The boundary condition equations are placed into the BoundaryConditions. */
    //options->SetString     ("Model.Evaluators.Device_0.Library", "Sequential");
    //options->SetString     ("Model.Evaluators.Device_0.Name",    "seq");
    //options->SetStringArray("Model.Evaluators.Device_0.Groups",  grps);
    /* (2a) Kernels are evaluated in parallel using OpenMP, one at the time. */
    //options->SetString     ("Model.Evaluators.Device_1.Library",               "Kernels_OpenMP");
    //options->SetString     ("Model.Evaluators.Device_1.Name",                  "kopenomp");
    //options->SetInteger    ("Model.Evaluators.Device_1.Parameters.numThreads", 0);
    //options->SetStringArray("Model.Evaluators.Device_1.Kernels",               krnls);
    /* (2b) Kernels are evaluated in parallel using OpenCL, one at the time. */
    //options->SetString("Model.Evaluators.Device_1.Library", "Kernels_OpenCL");
    //options->SetString("Model.Evaluators.Device_1.Name",    "kopencl");
    //options->SetInteger("Model.Evaluators.Device_1.Parameters.platformID", 0);
    //options->SetInteger("Model.Evaluators.Device_1.Parameters.deviceID",   0);
    //options->SetString("Model.Evaluators.Device_1.Parameters.API",        "C99");
    //options->SetStringArray("Model.Evaluators.Device_1.Kernels", krnls);

    // For Ncpu = 1: k = 3, rho = 1.0, alpha = 1e-1, w = 0.5
    options->SetInteger("LinearSolver.Preconditioner.Parameters.fact: level-of-fill",         3);
    options->SetDouble ("LinearSolver.Preconditioner.Parameters.fact: relax value",         0.5);
    options->SetDouble ("LinearSolver.Preconditioner.Parameters.fact: absolute threshold", 1e-1);
    options->SetDouble ("LinearSolver.Preconditioner.Parameters.fact: relative threshold",  1.0);
    std::string simulationOptions_seq = options->ToString();

    // Generate a single model (no graph partitioner required).
    std::string inputFilesDirectory = "dae_example_3_kernels-sequential";
    inputFilesDirectory += "-" + std::to_string(Nx) + "x" + std::to_string(Ny);
    {
        uint32_t Npe = 1;
        std::vector<std::string> balancingConstraints;
        csGraphPartitioner_t* gp = NULL;
        std::vector<csModelPtr> models_sequential = mb.PartitionSystem(Npe, gp, balancingConstraints, true);
        csModelBuilder_t::ExportModels(models_sequential, inputFilesDirectory, simulationOptions_seq);
        printf("Generated %s model for Npe = 1\n", inputFilesDirectory.c_str());
    }

    /* 4. Simulate the sequential model using the libOpenCS_Simulators. */
    printf("Simulation of '%s' (using libOpenCS_Simulators)\n\n", inputFilesDirectory.c_str());
    MPI_Init(&argc, &argv);
    csSimulate(inputFilesDirectory);
    MPI_Finalize();

    return 0;
}