Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Mmhllc #629

Open
wants to merge 9 commits into
base: develop
Choose a base branch
from
+266 −21
Open

Mmhllc #629

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
a0abc90
simplified code for mm hllc
adityakpandare Oct 21, 2024
0638f99
HLLCSolids works on single material fluids in the same way as HLLC
airaudofacundo Oct 22, 2024
436474f
Implemented multi-fluid HLLCSolids. Does not work
airaudofacundo Oct 23, 2024
be9c662
Implemented HLLCMultiMat, not working yet
airaudofacundo Oct 23, 2024
b136206
Merge branch 'develop' into mmhllc_fluids_test
adityakpandare Oct 23, 2024
e14b08b
Partially fixed multi-fluid. Using averaged pStar
airaudofacundo Oct 23, 2024
5026661
Merge branch 'mmhllc_fluids_test' of https://github.com/quinoacomputi…
airaudofacundo Oct 23, 2024
d668990
indexing bug fix
adityakpandare Oct 24, 2024
ddde629
use independently limited material pressures, this fixes the big spik…
adityakpandare Oct 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/PDE/MultiMat/RiemannChoice.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
#include "Riemann/HLL.hpp"
#include "Riemann/AUSM.hpp"
#include "Riemann/LaxFriedrichsSolids.hpp"
#include "Riemann/HLLCSolids.hpp"
#include "Riemann/HLLCMultiMat.hpp"

namespace inciter {

Expand All @@ -41,7 +41,7 @@ namespace inciter {
fluxfn = LaxFriedrichsSolids::flux;
}
else if (flux == ctr::FluxType::HLLC) {
fluxfn = HLLCSolids::flux;
fluxfn = HLLCMultiMat::flux;
}
else {
Throw("Riemann solver not set up for multi-material PDEs.");
Expand Down
237 changes: 237 additions & 0 deletions src/PDE/Riemann/HLLCMultiMat.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,237 @@
// *****************************************************************************
/*!
\file src/PDE/Riemann/HLLCMultiMat.hpp
\copyright 2012-2015 J. Bakosi,
2016-2018 Los Alamos National Security, LLC.,
2019-2023 Triad National Security, LLC.
All rights reserved. See the LICENSE file for details.
\brief HLLC Riemann flux function for solids
\details This file implements the HLLC Riemann solver for solids.
Ref. Ndanou, S., Favrie, N., & Gavrilyuk, S. (2015). Multi-solid
and multi-fluid diffuse interface model: Applications to dynamic
fracture and fragmentation. Journal of Computational Physics, 295,
523-555.
*/
// *****************************************************************************
#ifndef HLLCMultiMat_h
#define HLLCMultiMat_h

#include <vector>

#include "Fields.hpp"
#include "FunctionPrototypes.hpp"
#include "Inciter/Options/Flux.hpp"
#include "EoS/EOS.hpp"
#include "MultiMat/MultiMatIndexing.hpp"
#include "MultiMat/MiscMultiMatFns.hpp"

namespace inciter {

//! HLLC approximate Riemann solver for solids
struct HLLCMultiMat {

//! HLLC approximate Riemann solver flux function
//! \param[in] fn Face/Surface normal
//! \param[in] u Left and right unknown/state vector
//! \return Riemann solution according to Harten-Lax-van Leer-Contact
//! \note The function signature must follow tk::RiemannFluxFn
static tk::RiemannFluxFn::result_type
flux( const std::vector< EOS >& mat_blk,
const std::array< tk::real, 3 >& fn,
const std::array< std::vector< tk::real >, 2 >& u,
const std::vector< std::array< tk::real, 3 > >& = {} )
{
auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();

auto nsld = numSolids(nmat, solidx);
auto ncomp = u[0].size()-(3+nmat+nsld*6);
std::vector< tk::real > flx(ncomp, 0), fl(ncomp, 0), fr(ncomp, 0),
ftl(ncomp, 0), ftr(ncomp, 0);

// Primitive variables
// -------------------------------------------------------------------------
tk::real rhol(0.0), rhor(0.0);
for (size_t k=0; k<nmat; ++k) {
rhol += u[0][densityIdx(nmat, k)];
rhor += u[1][densityIdx(nmat, k)];
}

auto ul = u[0][ncomp+velocityIdx(nmat, 0)];
auto vl = u[0][ncomp+velocityIdx(nmat, 1)];
auto wl = u[0][ncomp+velocityIdx(nmat, 2)];
auto ur = u[1][ncomp+velocityIdx(nmat, 0)];
auto vr = u[1][ncomp+velocityIdx(nmat, 1)];
auto wr = u[1][ncomp+velocityIdx(nmat, 2)];

// Outer states
// -------------------------------------------------------------------------
tk::real pl(0.0), pr(0.0);
std::vector< tk::real > apl(nmat, 0.0), apr(nmat, 0.0);
tk::real acl(0.0), acr(0.0);

for (std::size_t k=0; k<nmat; ++k) {
// Left state
apl[k] = u[0][ncomp+pressureIdx(nmat, k)];
pl += apl[k];
auto amatl = mat_blk[k].compute< EOS::soundspeed >(
u[0][densityIdx(nmat, k)], apl[k], u[0][volfracIdx(nmat, k)], k );

// Right state
apr[k] = u[1][ncomp+pressureIdx(nmat, k)];
pr += apr[k];
auto amatr = mat_blk[k].compute< EOS::soundspeed >(
u[1][densityIdx(nmat, k)], apr[k], u[1][volfracIdx(nmat, k)], k );

// Mixture speed of sound
acl += u[0][densityIdx(nmat, k)] * amatl * amatl;
acr += u[1][densityIdx(nmat, k)] * amatr * amatr;
}
acl = std::sqrt(acl/rhol);
acr = std::sqrt(acr/rhor);

// Face-normal velocities
tk::real vnl = ul*fn[0] + vl*fn[1] + wl*fn[2];
tk::real vnr = ur*fn[0] + vr*fn[1] + wr*fn[2];

// Signal velocities
auto Sl = std::min((vnl-acl), (vnr-acr));
auto Sr = std::max((vnl+acl), (vnr+acr));
auto Sm = ( rhor*vnr*(Sr-vnr) - rhol*vnl*(Sl-vnl) + pl-pr )
/( rhor*(Sr-vnr) - rhol*(Sl-vnl) );

// Middle-zone (star) variables
// -------------------------------------------------------------------------
tk::real pStar(0.0);
std::vector< tk::real > apStar(nmat, 0.0);
for (std::size_t k=0; k<nmat; ++k) {
apStar[k] = 0.5 * (apl[k] + u[0][densityIdx(nmat, k)]*(vnl-Sl)*(vnl-Sm)
+ apr[k] + u[1][densityIdx(nmat, k)]*(vnr-Sr)*(vnr-Sm));
pStar += apStar[k];
}
auto uStar = u;

for (std::size_t i=0; i<3; ++i) {
uStar[0][momentumIdx(nmat, i)] =
((Sl-vnl)*u[0][momentumIdx(nmat, i)] + (pStar-pl)*fn[i])/(Sl-Sm);
uStar[1][momentumIdx(nmat, i)] =
((Sr-vnr)*u[1][momentumIdx(nmat, i)] + (pStar-pr)*fn[i])/(Sr-Sm);
}

for (std::size_t k=0; k<nmat; ++k) {
// Left
uStar[0][volfracIdx(nmat, k)] = u[0][volfracIdx(nmat, k)];
uStar[0][densityIdx(nmat, k)] =
(Sl-vnl) * u[0][densityIdx(nmat, k)] / (Sl-Sm);
uStar[0][energyIdx(nmat, k)] =
((Sl-vnl)*u[0][energyIdx(nmat,k)]
+ u[0][densityIdx(nmat,k)]*(Sl-vnl)*(Sm-vnl)*Sm
+ (Sm-vnl)*apl[k]) / (Sl-Sm);

// Right
uStar[1][volfracIdx(nmat, k)] = u[1][volfracIdx(nmat, k)];
uStar[1][densityIdx(nmat, k)] =
(Sr-vnr) * u[1][densityIdx(nmat, k)] / (Sr-Sm);
uStar[1][energyIdx(nmat, k)] =
((Sr-vnr)*u[1][energyIdx(nmat,k)]
+ u[1][densityIdx(nmat,k)]*(Sr-vnr)*(Sm-vnr)*Sm
+ (Sm-vnr)*apr[k]) / (Sr-Sm);
}

// Numerical fluxes
// -------------------------------------------------------------------------
if (Sl > 0.0) {

for (std::size_t idir=0; idir<3; ++idir)
flx[momentumIdx(nmat, idir)] =
u[0][momentumIdx(nmat, idir)] * vnl + pl*fn[idir];

for (std::size_t k=0; k<nmat; ++k) {
flx[volfracIdx(nmat, k)] = u[0][volfracIdx(nmat, k)] * vnl;
flx[densityIdx(nmat, k)] = u[0][densityIdx(nmat, k)] * vnl;
flx[energyIdx(nmat, k)] = (u[0][energyIdx(nmat, k)] + apl[k]) * vnl;
}

// Quantities for non-conservative terms
// Store Riemann-advected partial pressures
for (std::size_t k=0; k<nmat; ++k)
flx.push_back(apl[k]);
// Store Riemann velocity
flx.push_back(vnl);
}

else if (Sl <= 0.0 && Sm > 0.0) {

for (std::size_t idir=0; idir<3; ++idir)
flx[momentumIdx(nmat, idir)] =
uStar[0][momentumIdx(nmat, idir)] * Sm + pStar*fn[idir];

for (std::size_t k=0; k<nmat; ++k) {
flx[volfracIdx(nmat, k)] = uStar[0][volfracIdx(nmat, k)] * Sm;
flx[densityIdx(nmat, k)] = uStar[0][densityIdx(nmat, k)] * Sm;
flx[energyIdx(nmat, k)] = (uStar[0][energyIdx(nmat, k)] + apStar[k]) * Sm;
}

// Quantities for non-conservative terms
// Store Riemann-advected partial pressures
for (std::size_t k=0; k<nmat; ++k)
flx.push_back(apStar[k]);
// Store Riemann velocity
flx.push_back(Sm);
}

else if (Sm <= 0.0 && Sr >= 0.0) {

for (std::size_t idir=0; idir<3; ++idir)
flx[momentumIdx(nmat, idir)] =
uStar[1][momentumIdx(nmat, idir)] * Sm + pStar*fn[idir];

for (std::size_t k=0; k<nmat; ++k) {
flx[volfracIdx(nmat, k)] = uStar[1][volfracIdx(nmat, k)] * Sm;
flx[densityIdx(nmat, k)] = uStar[1][densityIdx(nmat, k)] * Sm;
flx[energyIdx(nmat, k)] = (uStar[1][energyIdx(nmat, k)] + apStar[k]) * Sm;
}

// Quantities for non-conservative terms
// Store Riemann-advected partial pressures
for (std::size_t k=0; k<nmat; ++k)
flx.push_back(apStar[k]);
// Store Riemann velocity
flx.push_back(Sm);
}

else {

for (std::size_t idir=0; idir<3; ++idir)
flx[momentumIdx(nmat, idir)] =
u[1][momentumIdx(nmat, idir)] * vnr + pr*fn[idir];

for (std::size_t k=0; k<nmat; ++k) {
flx[volfracIdx(nmat, k)] = u[1][volfracIdx(nmat, k)] * vnr;
flx[densityIdx(nmat, k)] = u[1][densityIdx(nmat, k)] * vnr;
flx[energyIdx(nmat, k)] = (u[1][energyIdx(nmat, k)] + apr[k]) * vnr;
}

// Quantities for non-conservative terms
// Store Riemann-advected partial pressures
for (std::size_t k=0; k<nmat; ++k)
flx.push_back(apr[k]);
// Store Riemann velocity
flx.push_back(vnr);
}

Assert( flx.size() == (ncomp+nmat+1+3*nsld), "Size of "
"multi-material flux vector incorrect" );

return flx;
}

////! Flux type accessor
////! \return Flux type
//static ctr::FluxType type() noexcept {
// return ctr::FluxType::HLLCMultiMat; }
};

} // inciter::

#endif // HLLCMultiMat_h
46 changes: 27 additions & 19 deletions src/PDE/Riemann/HLLCSolids.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -158,8 +158,8 @@ struct HLLCSolids {
ac_r = std::sqrt(ac_r/rhor);

// Signal velocities
auto Sl = std::min((vn_l[0]-ac_l), (0.5*(vn_r[0]+vn_l[0]-ac_r-ac_l)));
auto Sr = std::max((vn_l[0]+ac_l), (0.5*(vn_r[0]+vn_l[0]+ac_r+ac_l)));
auto Sl = std::min((vn_l[0]-ac_l), (vn_r[0]-ac_r));
auto Sr = std::max((vn_l[0]+ac_l), (vn_r[0]+ac_r));
auto Si = (rhol*vn_l[0]*(Sl-vn_l[0])
-rhor*vn_r[0]*(Sr-vn_r[0])
+signn_l[0][0]-signn_r[0][0])
Expand Down Expand Up @@ -200,7 +200,7 @@ struct HLLCSolids {
asig_t[k][j] = ( (vn_r[0]-Sr)*u[1][densityIdx(nmat,k)]*asignn_l[k][0][j]
- (vn_l[0]-Sl)*u[0][densityIdx(nmat,k)]*asignn_r[k][0][j]
+ (vn_l[0]-Sl)*u[0][densityIdx(nmat,k)]
*(vn_r[0]-Sr)*u[1][densityIdx(nmat,k)]
* (vn_r[0]-Sr)*u[1][densityIdx(nmat,k)]
* (vn_r[j]-vn_l[j]) ) /
( (vn_r[0]-Sr)*u[1][densityIdx(nmat,k)] -
(vn_l[0]-Sl)*u[0][densityIdx(nmat,k)]);
Expand All @@ -213,6 +213,7 @@ struct HLLCSolids {
rho_t_l += arho_t_l[k];

// inv deformation gradient
if (solidx[k] > 0) {
gn_t_l = gn_l[k];
for (std::size_t i=0; i<3; ++i)
gn_t_l[i][0] = w_l*gn_l[k][i][0] + (
Expand All @@ -221,6 +222,7 @@ struct HLLCSolids {

// rotate g back to original frame of reference
g_t_l.push_back(tk::unrotateTensor(gn_t_l, fn));
}

// energy
arhoe_t_l[k] = u[0][energyIdx(nmat, k)] * w_l + (
Expand All @@ -240,6 +242,7 @@ struct HLLCSolids {
rho_t_r += arho_t_r[k];

// inv deformation gradient
if (solidx[k] > 0) {
gn_t_r = gn_r[k];
for (std::size_t i=0; i<3; ++i)
gn_t_r[i][0] = w_r*gn_r[k][i][0] + (
Expand All @@ -248,6 +251,7 @@ struct HLLCSolids {

// rotate g back to original frame of reference
g_t_r.push_back(tk::unrotateTensor(gn_t_r, fn));
}

// energy
arhoe_t_r[k] = u[1][energyIdx(nmat, k)] * w_r + (
Expand All @@ -264,6 +268,8 @@ struct HLLCSolids {
v_t_l = tk::unrotateVector(vn_t_l, fn);
v_t_r = tk::unrotateVector(vn_t_r, fn);

sig_t = tk::unrotateVector(sig_t, fn);

// -------------------------------------------------------------------------

// Flux functions based on HLLC
Expand Down Expand Up @@ -320,15 +326,15 @@ struct HLLCSolids {
}

// Star state fluxes
ftl = fl;
ftr = fr;
for (std::size_t k=0; k<nmat; ++k)
{
// Left fluxes
ftl[volfracIdx(nmat, k)] += Sl * (al_t_l[k] - al_l[k]);
ftl[densityIdx(nmat, k)] +=
Sl * (arho_t_l[k] - u[0][densityIdx(nmat, k)]);
ftl[energyIdx(nmat, k)] += Sl * (arhoe_t_l[k] - u[0][energyIdx(nmat, k)]);
ftl[volfracIdx(nmat, k)] = vn_t_l[0] * al_t_l[k];
ftl[densityIdx(nmat, k)] = vn_t_l[0] * arho_t_l[k];
ftl[energyIdx(nmat, k)] = vn_t_l[0] * arhoe_t_l[k] - (//(sig_11*u + sig_12*v + sig_13*w);
+ asig_t[k][0]*vn_t_l[0]
+ asig_t[k][1]*vn_t_l[1]
+ asig_t[k][2]*vn_t_l[2] );

// inv deformation gradient tensor
if (solidx[k] > 0) {
Expand All @@ -339,10 +345,12 @@ struct HLLCSolids {
}

// Right fluxes
ftr[volfracIdx(nmat, k)] += Sr * (al_t_r[k] - al_r[k]);
ftr[densityIdx(nmat, k)] +=
Sr * (arho_t_r[k] - u[1][densityIdx(nmat, k)]);
ftr[energyIdx(nmat, k)] += Sr * (arhoe_t_r[k] - u[1][energyIdx(nmat, k)]);
ftr[volfracIdx(nmat, k)] = vn_t_r[0] * al_t_r[k];
ftr[densityIdx(nmat, k)] = vn_t_r[0] * arho_t_r[k];
ftr[energyIdx(nmat, k)] = vn_t_r[0] * arhoe_t_r[k] - (//(sig_11*u + sig_12*v + sig_13*w);
+ asig_t[k][0]*vn_t_r[0]
+ asig_t[k][1]*vn_t_r[1]
+ asig_t[k][2]*vn_t_r[2] );

// inv deformation gradient tensor
if (solidx[k] > 0) {
Expand All @@ -355,10 +363,10 @@ struct HLLCSolids {
// bulk momentum
for (std::size_t idir=0; idir<3; ++idir)
{
ftl[momentumIdx(nmat, idir)] +=
Sl * (rho_t_l*v_t_l[idir] - u[0][momentumIdx(nmat, idir)]);
ftr[momentumIdx(nmat, idir)] +=
Sr * (rho_t_r*v_t_r[idir] - u[1][momentumIdx(nmat, idir)]);
ftl[momentumIdx(nmat, idir)] =
rho_t_l*v_t_l[idir]*Si - sig_t[idir];
ftr[momentumIdx(nmat, idir)] =
rho_t_r*v_t_r[idir]*Si - sig_t[idir];
}

// -------------------------------------------------------------------------
Expand Down Expand Up @@ -394,7 +402,7 @@ struct HLLCSolids {
for (std::size_t k=0; k<nmat; ++k)
flx.push_back(-(aTn_l[k][0]+Sl/(Si-Sl)*(aTn_l[k][0]-asig_t[k][0])));
// Store Riemann velocity
flx.push_back( vn_l[0] );
flx.push_back( vn_t_l[0] );
// Store Riemann aTn_ij (3*nsld)
for (std::size_t k=0; k<nmat; ++k) {
if (solidx[k] > 0) {
Expand All @@ -413,7 +421,7 @@ struct HLLCSolids {
for (std::size_t k=0; k<nmat; ++k)
flx.push_back(-(aTn_r[k][0]+Sr/(Si-Sr)*(aTn_r[k][0]-asig_t[k][0])));
// Store Riemann velocity
flx.push_back( vn_r[0] );
flx.push_back( vn_t_r[0] );
// Store Riemann aTn_ij (3*nsld)
for (std::size_t k=0; k<nmat; ++k) {
if (solidx[k] > 0) {
Expand Down