src/pennant/src/Hydro.cc

/*
 * Hydro.cc
 *
 *  Created on: Dec 22, 2011
 *      Author: cferenba
 *
 * Copyright (c) 2012, Los Alamos National Security, LLC.
 * All rights reserved.
 * Use of this source code is governed by a BSD-style open-source
 * license; see top-level LICENSE file for full license text.
 */

#include "Hydro.hh"

#include <string>
#include <vector>
#include <cmath>
#include <algorithm>
#include <sstream>
#include <iomanip>

#include "Memory.hh"
#include "InputFile.hh"
#include "Mesh.hh"
#include "PolyGas.hh"
#include "TTS.hh"
#include "QCS.hh"
#include "HydroBC.hh"
#include "HydroGPU.hh"

using namespace std;


Hydro::Hydro(const InputFile* inp, Mesh* m) : mesh(m) {
    cfl = inp->getDouble("cfl", 0.6);
    cflv = inp->getDouble("cflv", 0.1);
    rinit = inp->getDouble("rinit", 1.);
    einit = inp->getDouble("einit", 0.);
    rinitsub = inp->getDouble("rinitsub", 1.);
    einitsub = inp->getDouble("einitsub", 0.);
    uinitradial = inp->getDouble("uinitradial", 0.);
    bcx = inp->getDoubleList("bcx", vector<double>());
    bcy = inp->getDoubleList("bcy", vector<double>());

    pgas = new PolyGas(inp, this);
    tts = new TTS(inp, this);
    qcs = new QCS(inp, this);

    init();

    hydroInitGPU();
    hydroInit(
            mesh->nump, mesh->numz, mesh->nums, mesh->numc, mesh->nume,
            pgas->gamma, pgas->ssmin,
            tts->alfa, tts->ssmin,
            qcs->qgamma, qcs->q1, qcs->q2,
            cfl, cflv,
            bcx.size(), &bcx[0], bcy.size(), &bcy[0],
            mesh->px,
            pu,
            zm,
            zr,
            mesh->zvol,
            ze, zetot,
            zwrate,
            mesh->smf,
            mesh->mapsp1,
            mesh->mapsp2,
            mesh->mapsz,
            mesh->mapss4,
            mesh->mapse,
            mesh->znump);

}


Hydro::~Hydro() {

    hydroFinalGPU();

    delete tts;
    delete qcs;
    for (int i = 0; i < bcs.size(); ++i) {
        delete bcs[i];
    }
}


void Hydro::init() {

    dtrec = 1.e99;
    msgdtrec = "Hydro default";

    const int nump = mesh->nump;
    const int numz = mesh->numz;

    const double2* zx = mesh->zx;
    const double* zvol = mesh->zvol;

    // allocate arrays
    pu = Memory::alloc<double2>(nump);
    zm = Memory::alloc<double>(numz);
    zr = Memory::alloc<double>(numz);
    ze = Memory::alloc<double>(numz);
    zetot = Memory::alloc<double>(numz);
    zwrate = Memory::alloc<double>(numz);
    zp = Memory::alloc<double>(numz);

    // initialize hydro vars
    fill(&zr[0], &zr[numz], rinit);
    fill(&ze[0], &ze[numz], einit);
    fill(&zwrate[0], &zwrate[numz], 0.);

    const vector<double>& subrgn = mesh->subregion;
    if (!subrgn.empty()) {
        const double eps = 1.e-12;
        for (int z = 0; z < numz; ++z) {
            if (zx[z].x > (subrgn[0] - eps) &&
                zx[z].x < (subrgn[1] + eps) &&
                zx[z].y > (subrgn[2] - eps) &&
                zx[z].y < (subrgn[3] + eps)) {
                zr[z] = rinitsub;
                ze[z] = einitsub;
            }
        }
    }

    for (int z = 0; z < numz; ++z) {
        zm[z] = zr[z] * zvol[z];
        zetot[z] = ze[z] * zm[z];
    }

    if (uinitradial != 0.)
        initRadialVel(uinitradial);
    else
        fill(&pu[0], &pu[nump], double2(0., 0.));

}


void Hydro::getData() {

    hydroGetData(
            mesh->nump, mesh->numz,
            mesh->px,
            zr, ze, zp);

}


void Hydro::initRadialVel(const double vel) {
    const int nump = mesh->nump;
    const double2* px = mesh->px;
    const double eps = 1.e-12;

    for (int p = 0; p < nump; ++p) {
        double pmag = length(px[p]);
        if (pmag > eps)
            pu[p] = vel * px[p] / pmag;
        else
            pu[p] = double2(0., 0.);
    }
}


void Hydro::doCycle(
            const double dt) {

    int idtrec;

    hydroDoCycle(dt, dtrec, idtrec);

    int z = idtrec >> 1;
    bool dtfromvol = idtrec & 1;
    ostringstream oss;
    if (dtfromvol)
        oss << "Hydro dV/V limit for z = " << setw(6) << z;
    else
        oss << "Hydro Courant limit for z = " << setw(6) << z;
    msgdtrec = oss.str();

}


void Hydro::getDtHydro(
        double& dtnew,
        string& msgdtnew) {

    if (dtrec < dtnew) {
        dtnew = dtrec;
        msgdtnew = msgdtrec;
    }

}