diff --git a/src/PDE/MultiMat/RiemannChoice.hpp b/src/PDE/MultiMat/RiemannChoice.hpp
index 14ff5f2a26..bc2342dfe5 100644
--- a/src/PDE/MultiMat/RiemannChoice.hpp
+++ b/src/PDE/MultiMat/RiemannChoice.hpp
@@ -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 {
 
@@ -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.");
diff --git a/src/PDE/Riemann/HLLCMultiMat.hpp b/src/PDE/Riemann/HLLCMultiMat.hpp
new file mode 100644
index 0000000000..2ffd425ec8
--- /dev/null
+++ b/src/PDE/Riemann/HLLCMultiMat.hpp
@@ -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
diff --git a/src/PDE/Riemann/HLLCSolids.hpp b/src/PDE/Riemann/HLLCSolids.hpp
index c7d33eccb9..7424a12b3e 100644
--- a/src/PDE/Riemann/HLLCSolids.hpp
+++ b/src/PDE/Riemann/HLLCSolids.hpp
@@ -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])
@@ -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)]);
@@ -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] + (
@@ -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 + (
@@ -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] + (
@@ -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 + (
@@ -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
@@ -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) {
@@ -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) {
@@ -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];
     }
 
     // -------------------------------------------------------------------------
@@ -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) {
@@ -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) {