/*
 * PolyGas.cc
 *
 *  Created on: Mar 26, 2012
 *      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 "PolyGas.hh"

#include "Memory.hh"
#include "InputFile.hh"
#include "Hydro.hh"
#include "Mesh.hh"

using namespace std;


PolyGas::PolyGas(const InputFile* inp, Hydro* h) : hydro(h) {
    gamma = inp->getDouble("gamma", 5. / 3.);
    ssmin = inp->getDouble("ssmin", 0.);

}


void PolyGas::calcStateAtHalf(
        const double* zr0,
        const double* zvolp,
        const double* zvol0,
        const double* ze,
        const double* zwrate,
        const double* zm,
        const double dt,
        double* zp,
        double* zss,
        const int zfirst,
        const int zlast) {

    double* z0per = Memory::alloc<double>(zlast - zfirst);

    const double dth = 0.5 * dt;

    // compute EOS at beginning of time step
    calcEOS(zr0, ze, zp, z0per, zss, zfirst, zlast);

    // now advance pressure to the half-step
    #pragma ivdep
    for (int z = zfirst; z < zlast; ++z) {
        int z0 = z - zfirst;
        double zminv = 1. / zm[z];
        double dv = (zvolp[z] - zvol0[z]) * zminv;
        double bulk = zr0[z] * zss[z] * zss[z];
        double denom = 1. + 0.5 * z0per[z0] * dv;
        double src = zwrate[z] * dth * zminv;
        zp[z] += (z0per[z0] * src - zr0[z] * bulk * dv) / denom;
    }

    Memory::free(z0per);
}


void PolyGas::calcEOS(
        const double* zr,
        const double* ze,
        double* zp,
        double* z0per,
        double* zss,
        const int zfirst,
        const int zlast) {

    const double gm1 = gamma - 1.;
    const double ss2 = max(ssmin * ssmin, 1.e-99);

    #pragma ivdep
    for (int z = zfirst; z < zlast; ++z) {
        int z0 = z - zfirst;
        double rx = zr[z];
        double ex = max(ze[z], 0.0);
        double px = gm1 * rx * ex;
        double prex = gm1 * ex;
        double perx = gm1 * rx;
        double csqd = max(ss2, prex + perx * px / (rx * rx));
        zp[z] = px;
        z0per[z0] = perx;
        zss[z] = sqrt(csqd);
    }

}


void PolyGas::calcForce(
        const double* zp,
        const double2* ssurfp,
        double2* sf,
        const int sfirst,
        const int slast) {

    const Mesh* mesh = hydro->mesh;

    #pragma ivdep
    for (int s = sfirst; s < slast; ++s) {
        int z = mesh->mapsz[s];
        double2 sfx = -zp[z] * ssurfp[s];
        sf[s] = sfx;

    }
}