main.edp

include "getARGV.idp"
include "solvers.edp"
load "medit";

// Config parameters
int config = getARGV("--config", 1);                        // 1: Navier-Stokes, 2: Natural convection
int testcase = getARGV("--testcase", 1);                    // the case of test (1, 2 and so on)
string meshname = getARGV("--meshname", "liddriven19");     // the name of mesh file
int dc = getARGV("--dc", 0);                                // the config of DC methods (1 for DC or 0 for non-DC)
real nu = getARGV("--nu", 1e-3);                            // the constant kinematic viscosity of the fluid (Navier-Stokes only) 
real Pr = getARGV("--Pr", 1.0);							    // the Prandtl number (natural convection only)
real Ra = getARGV("--Ra", 100000.0);	                    // the Rayleigh number (natural convection only)
real tf = getARGV("--tf", 100.0);                           // the final time (second)
real dt = getARGV("--dt", 0.01);;                           // the time step (second)
int save = getARGV("--save", 0);                            // save results to .sol files configuration (1 for save, 0 for non-save)
int ns = getARGV("--ns", 100);                              // number of steps to save results once

if(config == 2) 
    nu = Pr;

// Make output
string outputFolder = getARGV("--resu", "results\Lid_driven_cavity");
cout << "Results and figures will be saved in " << outputFolder << endl;
system("mkdir "+outputFolder);


// Save the command 
ofstream cmd(outputFolder+"/command.sh");
for (int ii = 0; ii < ARGV.n; ii++)
  cmd << ARGV[ii] << " ";
cmd << endl;
cmd.flush;


// Load computational mesh
mesh Th;
string mname = "meshes/" + meshname + ".mesh";
cout << "Loading mesh " << meshname << "...";
Th = readmesh(mname);
cout << "done." << endl;
cout.flush;


real sigma = getARGV("--sigma", 0.4 * Th.hmax);                             // the stabilizing factor of defect-correction method (DC only)


// Finite element spaces
fespace Xh(Th,P2);                                                          // definition of the velocity component space
fespace Mh(Th,P1);                                                          // definition of the pressure space


// Declare variables
Xh ux, uy, vx, vy, upx, upy, u0x, u0y, up0x, up0y;
Xh T, Tp, Tau;
Xh dcx1, dcx2, dcx3, dcx4, dcy1, dcy2, dcy3, dcy4, w, psi, phi;
Xh tbc1, tbc2, tbc3, tbc4;
Mh p, q, p0;
Xh fx = 0.0, fy = 0.0, Q = 0.0;


// Set boundary conditions
if(config == 1)
{
    if(testcase == 1)                                 // Lid - driven cavity
    {
        dcx1 = 1.;
        dcy1 = 0.;
        dcx2 = 0.;
        dcy2 = 0.;
    }
    else if(testcase == 2)                            // Backward facing step
    {
        dcx1 = 0.5 * y * (2 - y);
        dcx2 = 0.;
        dcx2 = 0.;
        dcy2 = 0.;
    }
    else if (testcase == 3)                           // Flow around cylinder
    {                     
        real D = 0.1, H = 0.41;
        real Um = getARGV("--Um", 1.5);
        dcx1 = 4. * Um * y * (H - y) / (H^2);
        dcy1 = 0.;
        dcx2 = 0.;
        dcy2 = 0.;
        int nn = getARGV("--nn", 30);               // number of meshing points                  
        dt = D / nn / Um;
    }
}
else if(config == 2)
{
    ux = 0.0;
    uy = 0.0;
    T = 0.0; 

    if(testcase == 1)                               // Thermal driven cavity
    {
        dcx2 = 0.0;
        dcx4 = 0.0;
        dcy2 = 0.0;
        dcy4 = 0.0;

        tbc2 = -0.5;
        tbc4 = 0.5;
    }
    else if(testcase == 2)                          // Benard test case
    {
        dcx1 = 0.0;
        dcx3 = 0.0;
        dcy1 = 0.0;
        dcy3 = 0.0;

        tbc1 = 1.0;
        tbc3 = 0.0;
    }
}


real alpha = 1.0 / dt;
int i = 0;
int imax = tf / dt;

// Main code

for (i = 1; i <= imax; i++)
{
    upx = ux;
    upy = uy;
    if(dc == 1)
    {
        up0x = u0x;
        up0y = u0y;
    }

    if(config == 1)
    {
        NavierStokes;
    }
    else if(config == 2)
    {
        Tp = T;
        NaturalConvection;

        // Only for Thermal driven cavity test
        real currtime = i * dt;
        if(testcase == 1 && (currtime == 0.003 || currtime == 0.01 || currtime == 0.025 || currtime == 0.1 || currtime == 0.2))
        {
            ofstream re(outputFolder + "/maxvelocity." + currtime + ".txt");
            Xh unorm = sqrt(ux^2 + uy^2);
            re << unorm[].max << endl;
            re.flush;
        }
    }
        
    Streamlines;
    w = -dy(ux) + dx(uy);

    if(i % 10 == 0)
    {
        plot([ux, uy], cmm = "Velocity: iteration " + i + "/" + imax);
        plot(p, cmm = "Pressure: iteration " + i + "/" + imax);	
        plot(w, cmm = "Vorticity: iteration " + i + "/" + imax);		
        plot(psi, cmm = "Streamlines: iteration " + i + "/" + imax);	

        if(config == 2)
        {
            plot(T, cmm = "Thermal: iteration " + i + "/" + imax);
        }		
    }

    if(save == 1 && i % ns == 0)
    {
        savesol(outputFolder + "/velocity." + i/ns + ".sol", Th, [ux, uy]);	
        savemesh(Th, outputFolder + "/velocity." + i/ns + ".mesh");	

        savesol(outputFolder + "/pressure." + i/ns + ".sol", Th, p);
        savemesh(Th, outputFolder + "/pressure." + i/ns + ".mesh");

        savesol(outputFolder + "/vorticity." + i/ns + ".sol", Th, w);
        savemesh(Th, outputFolder + "/vorticity." + i/ns +".mesh");

        savesol(outputFolder + "/streamlines." + i/ns + ".sol", Th, psi);
        savemesh(Th, outputFolder + "/streamlines." + i/ns + ".mesh");

        if(config == 2)
        {
            savesol(outputFolder + "/thermal." + i/ns + ".sol", Th, T);
            savemesh(Th, outputFolder + "/thermal." + i/ns + ".mesh");
        }
    }
}


// Only for Lid-driven cavity test
int nbp = getARGV("--nbp", 100);                            // the number of points in vertical and horizontal lines

if (config == 1 && testcase == 1)
{
    // Get velocity of flow at points in vertical and horizontal lines that pass center of cavity
    real dd = 1.0 / nbp;

    ofstream fxprof(outputFolder + "/xprof");
    for(i = 0; i <= nbp; i++)
        fxprof << i * dd << " " << ux(0.5, i * dd) << endl;
    fxprof.flush;

    ofstream fyprof(outputFolder + "/yprof");
    for(i = 0; i <= nbp; i++)
        fyprof << i * dd << " " << uy(i * dd, 0.5)  << endl;
    fyprof.flush;


    // Find center of primary vortex
    real minx = 0.5;
    real maxx = 0.8;
    real miny = 0.5;
    real maxy = 0.8;

    real centerx = minx;
    real centery = miny;

    int ni = getARGV("--ni", 300);
    real slx = (maxx - minx)/ni;
    real sly = (maxy - miny)/ni;

    real minu = ux(minx, miny)^2 + uy(minx, miny)^2;

    for (int j = 1; j <= ni; j++)
        for(int k = 1; k <= ni; k++)
        {
            real tempx = minx + j * slx;
            real tempy = miny + k * sly;
            real tempu = ux(tempx, tempy)^2 + uy(tempx, tempy)^2;
            if(tempu < minu)
                {
                    minu = tempu;
                    centerx = tempx;
                    centery = tempy;
                }
        }

    ofstream ct(outputFolder + "/center.txt");
    ct << centerx << " " << centery  << endl;
    ct.flush;
}