// This is a sandbox example to get a feeling with the stabilization methods
// proposed in sparselizard for advection-dominated advection-diffusion problems.
//
// The geometry is a square of side 'a' with 'n' mesh points in the x and y direction.


#include "sparselizardbase.h"


using namespace mathop;

mesh createmesh(double a, int n);

void sparselizard()
{	
    // Square surface and its four sides:
    int sur = 1, left = 2, right = 3, up = 4, down = 5;

    mesh mymesh = createmesh(1.0, 101);

    // Field c can be a chemical concentration:
    field c("h1"), x("x"), y("y");

    // Use interpolation order 2:
    c.setorder(sur, 2);

    // Steady supply of chemical concentration from the left side: 
    c.setconstraint(left, 1);

    // A 0.1 m/s fluid flow in the x direction is considered:
    expression v = array2x1(0.1,0);

    // Initial all zero concentration except in the square 
    // delimited by [0,0.25] for x and [0.4,0.6] for y:
    expression cinit = ifpositive(x-0.25, 0.0, ifpositive(abs(y-0.5)-0.1, 0.0, 1.0) );

    c.setvalue(sur, cinit);
    c.write(sur, "cinit.vtu", 2);

    // Tuning factor for the stabilization:
    double delta = 0.1;
    
    // Timestep and end time for the simulation:
    double ts = 0.1, tend = 5.0;	
    
    // This will hold the solution vectors for every timestep:
    std::vector<vec> sol;


    std::cout << std::endl << "No stabilization:" << std::endl;
    formulation adnostab;

    adnostab += integral(sur, predefinedadvectiondiffusion(dof(c), tf(c), v, 1e-6, 1.0, 1.0, true));
    
    c.setvalue(sur, cinit);
    vec init(adnostab);
    init.setdata(sur, c);
    genalpha genanostab(adnostab, init, vec(adnostab), vec(adnostab), {true,true,true,true});
    genanostab.setparameter(0.75);
    sol = genanostab.runlinear(0, ts, tend, 2)[0];
    c.setdata(sur, sol[sol.size()-1]);
    c.write(sur, "cnostab.vtu", 2);


    std::cout << std::endl << "Isotropic:" << std::endl;
    formulation adiso;

    adiso += integral(sur, predefinedadvectiondiffusion(dof(c), tf(c), v, 1e-6, 1.0, 1.0, true));
    adiso += integral(sur, predefinedstabilization("iso", delta, c, v, 0.0, 0.0));

    c.setvalue(sur, cinit);
    vec initiso(adiso);
    initiso.setdata(sur, c);
    genalpha genaiso(adiso, initiso, vec(adiso), vec(adiso), {true,true,true,true});
    genaiso.setparameter(0.75);
    sol = genaiso.runlinear(0, ts, tend, 2)[0];
    c.setdata(sur, sol[sol.size()-1]);
    c.write(sur, "ciso.vtu", 2);


    std::cout << std::endl << "Streamline anisotropic:" << std::endl;
    formulation adaniso;

    adaniso += integral(sur, predefinedadvectiondiffusion(dof(c), tf(c), v, 1e-6, 1.0, 1.0, true));
    adaniso += integral(sur, predefinedstabilization("aniso", delta, c, v, 0.0, 0.0));

    c.setvalue(sur, cinit);
    vec initaniso(adaniso);
    initaniso.setdata(sur, c);
    genalpha genaaniso(adaniso, initaniso, vec(adaniso), vec(adaniso), {true,true,true,true});
    genaaniso.setparameter(0.75);
    sol = genaaniso.runlinear(0, ts, tend, 2)[0];
    c.setdata(sur, sol[sol.size()-1]);
    c.write(sur, "caniso.vtu", 2);


    // Define the strong form residual to be used by the stabilization methods below:
    expression dofres = v*grad(dof(c));
    expression res = v*grad(c);

    // Tuning parameter for streamline and for crosswind:
    double deltas = 0.1, deltac = 0.001;

    std::cout << std::endl << "Streamline and crosswind combined:" << std::endl;
    formulation adcomb;

    adcomb += integral(sur, predefinedadvectiondiffusion(dof(c), tf(c), v, 1e-6, 1.0, 1.0, true));
    adcomb += integral(sur, predefinedstabilization("supg", deltas, c, v, 1e-6, dofres));
    adcomb += integral(sur, predefinedstabilization("cws", deltac, c, v, 1e-6, res));
    
    c.setvalue(sur, cinit);
    vec initcomb(adcomb);
    initcomb.setdata(sur, c);

    genalpha genacomb(adcomb, initcomb, vec(adcomb), vec(adcomb), {true,true,true,true});
    genacomb.setparameter(0.75);
    sol = genacomb.runlinear(0, ts, tend, 2)[0];
    c.setdata(sur, sol[sol.size()-1]);
    c.write(sur, "ccombined.vtu", 2);
    

    // Code validation line. Can be removed.
    double cinterpolated = c.interpolate(sur, {0.76,0.399,0.0})[0];
    std::cout << (cinterpolated < 0.128046 && cinterpolated > 0.128043);
}


mesh createmesh(double a, int n)
{
    // Give names to the physical region numbers:
    int sur = 1, left = 2, right = 3, up = 4, down = 5;

    shape square("quadrangle", sur , {0,0,0, a,0,0, a,a,0, 0,a,0}, {n,n,n,n});

    shape line1 = square.getsons()[0];
    line1.setphysicalregion(down);
    shape line2 = square.getsons()[1];
    line2.setphysicalregion(right);
    shape line3 = square.getsons()[2];
    line3.setphysicalregion(up);
    shape line4 = square.getsons()[3];
    line4.setphysicalregion(left);

    // Provide to the mesh all shapes of interest:
    mesh mymesh({square,line1,line2,line3,line4});

    return mymesh;
}

int main(void)
{	
    SlepcInitialize(0,{},0,0);

    sparselizard();

    SlepcFinalize();

    return 0;
}