Commit 9f5a390e authored by Kilian Lackhove's avatar Kilian Lackhove
Browse files

APE: Back to the constant c formulation

works, but the received fields are _very_ discontinuous
parent 81d49533
......@@ -213,25 +213,20 @@ void APE::GetFluxVector(
ASSERTL1(flux[0].num_elements() == m_spacedim,
"Dimension of flux array and velocity array do not match");
// F_{adv,p',j} = \rho_0 u'_j + p' \bar{u}_j / c^2
// F_{adv,p',j} = \gamma p_0 u'_j + p' \bar{u}_j
for (int j = 0; j < m_spacedim; ++j)
{
Vmath::Zero(nq, flux[0][j], 1);
// construct rho_0 u'_j term
Vmath::Vmul(nq, m_basefield[1], 1, physfield[j + 1], 1, flux[0][j], 1);
// construct \gamma p_0 u'_j term
Vmath::Smul(nq, m_gamma, m_basefield[0], 1, tmp1, 1);
Vmath::Vmul(nq, tmp1, 1, physfield[j+1], 1, tmp1, 1);
// construct p' \bar{u}_j / c^2 term
// c^2
Vmath::Vdiv(nq, m_basefield[0], 1, m_basefield[1], 1, tmp1, 1);
Vmath::Smul(nq, m_gamma, tmp1, 1, tmp1, 1);
// construct p' \bar{u}_j term
Vmath::Vmul(nq, physfield[0], 1, m_basefield[j+2], 1, tmp2, 1);
// p' \bar{u}_j / c^2 term
Vmath::Vmul(nq, physfield[0], 1, m_basefield[j + 2], 1, tmp2, 1);
Vmath::Vdiv(nq, tmp2, 1, tmp1, 1, tmp2, 1);
// \rho_0 u'_j + p' \bar{u}_j / c^2
Vmath::Vadd(nq, flux[0][j], 1, tmp2, 1, flux[0][j], 1);
// add both terms
Vmath::Vadd(nq, tmp1, 1, tmp2, 1, flux[0][j], 1);
}
for (int i = 1; i < flux.num_elements(); ++i)
......@@ -306,22 +301,14 @@ void APE::DoOdeRhs(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
{
int nVariables = inarray.num_elements();
int nq = GetTotPoints();
Array<OneD, NekDouble> tmp1(nq);
// WeakDG does not use advVel, so we only provide a dummy array
Array<OneD, Array<OneD, NekDouble> > advVel(m_spacedim);
m_advection->Advect(nVariables, m_fields, advVel, inarray, outarray, time);
// Negate the LHS terms
for (int i = 0; i < nVariables; ++i)
{
if (i == 0)
{
// c^2 = gamma*p0/rho0
Vmath::Vdiv(nq, m_basefield[0], 1, m_basefield[1], 1, tmp1, 1);
Vmath::Smul(nq, m_gamma, tmp1, 1, tmp1, 1);
Vmath::Vmul(nq, tmp1, 1, outarray[i], 1, outarray[i], 1);
}
Vmath::Neg(nq, outarray[i], 1);
}
......
......@@ -95,18 +95,18 @@ void LaxFriedrichsSolver::v_PointSolve(
// max absolute eigenvalue of the jacobian of F_n1
NekDouble a_1_max = std::max(std::abs(u0 - c), std::abs(u0 + c));
NekDouble pFL = rho0*uL + u0*pL/(c*c);
NekDouble pFL = gamma*p0*uL + pL*u0;
NekDouble uFL = pL/rho0 + u0*uL + v0*vL + w0*wL;
NekDouble vFL = 0;
NekDouble wFL = 0;
NekDouble pFR = rho0*uR + u0*pR/(c*c);
NekDouble pFR = gamma*p0*uR + pR*u0;
NekDouble uFR = pR/rho0 + u0*uR + v0*vR + w0*wR;
NekDouble vFR = 0;
NekDouble wFR = 0;
// assemble the face-normal fluxes
pF = 0.5*(pFL + pFR - a_1_max*(pR/(c*c) - pL/(c*c)));
pF = 0.5*(pFL + pFR - a_1_max*(pR - pL));
uF = 0.5*(uFL + uFR - a_1_max*(uR - uL));
vF = 0.5*(vFL + vFR - a_1_max*(vR - vL));
wF = 0.5*(wFL + wFR - a_1_max*(wR - wL));
......
......@@ -123,7 +123,7 @@ void UpwindSolver::v_PointSolve(
NekDouble u = (W[0] - W[1])/(c*rho0);
// assemble the fluxes
pF = rho0*u + u0*p/(c*c);
pF = gamma*p0*u + u0*p;
uF = p/rho0 + u0*u + v0*vL + w0*wL;
vF = 0.0;
wF = 0.0;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment