[martinsshorts]

I am writing a FVM solver in C++ for elastic waves. I've got most of it working, but I cannot seem to get the Riemann solver to work. It is so unstable the flux values will blow up in just a few time steps.



The Riemann solver I am implementing can be seen here:

https://nbviewer.org/github/maojrs/i..._riemann.ipynb

It is from a Randy Leveque book.


For my own sanity, I need to switch to a simpler version of the Riemann solver so I can get it working. This project is just a proof of concept, so it doesn't need to be the most cutting edge solver ever. It should just be very simple, and work for linear elastic waves on an unstructured mesh with heterogeneous material. This is for learning purposes so I want to code it myself. I don't want to use a library that already has a Riemann solver in it.

Does anyone have any resources that could help me find a simple Riemann solver for elastic waves? Thank you very much.


================================================== ================================================== ==


For reference, here is my current Riemann solver:

static void ElasticRiemann(
    PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n,
    const PetscScalar *uL, const PetscScalar *uR, PetscInt numConstants,
    const PetscScalar constants[], PetscScalar *flux, void *ctx)
{
  /* Input Parameters:
       dim          - The spatial dimension
       Nf           - The number of fields
       x            - The coordinates at a point on the interface
       n            - The normal vector to the interface
       uL           - The state vector to the left of the interface
       uR           - The state vector to the right of the interface
       numConstants - number of constant parameters
       constants    - constant parameters
       ctx          - optional user context
       flux         - output array of flux through the interface
  */

  Physics phys = (Physics) ctx;  // Cast context to phys struct
  
  PetscReal du[5]; // difference in stress and velocity values between the sides
  PetscReal detP, detS;
  PetscReal nn[dim];
  PetscReal a1, a2, a3, a4; // alpha values
  
  // Compute jumps in the states (uL is left, uR is right)
  // and store in du

  /* Normalize the normal vector */
  nn[0] = n[0];
  nn[1] = n[1];
  Normalize2Real(nn);

  for (PetscInt i = 0; i < 5; ++i) {
      du[i] = uR[i] - uL[i];
  }

  // obtain material paramter values for L and R sides
  // these are currently part of the solution field because PETSc doesn't 
  // support auxiliary fields for FVM. I realize this is not ideal.

  PetscReal lambdaL = uL[5];
  PetscReal muL = uL[6];
  PetscReal densityL = uL[7];
  PetscReal lambdaR = uR[5];
  PetscReal muR = uR[6];
  PetscReal densityR = uR[7];

  // Calculate cp (P-wave speed) and cs (S-wave speed)
  PetscReal bulkL = lambdaL + 2.0 * muL;
  PetscReal bulkR = lambdaR + 2.0 * muR;
  PetscReal cpL = PetscSqrtReal(bulkL / densityL);
  PetscReal csL = PetscSqrtReal(muL / densityL);
  PetscReal cpR = PetscSqrtReal(bulkR / densityR);
  PetscReal csR = PetscSqrtReal(muR / densityR);

  // Calculate useful multiples of the norm vector components
  PetscReal nx = nn[0];
  PetscReal ny = nn[1];
  PetscReal nx2 = nx * nx;
  PetscReal ny2 = ny * ny;
  PetscReal nxy = nx * ny;
  
  // Define Eigenvalues (wave speeds)
  PetscReal s[4] = {-cpL, cpR, -csL, csR};

  // Define the 4 eigenvectors (from columns of Matrix R)
  PetscReal r1[5] = {lambdaL + 2 * muL * nx2, lambdaL + 2 * muL * ny2,
                     2 * muL * nxy, nx * cpL, ny * cpL};
  PetscReal r2[5] = {lambdaR + 2 * muR * nx2, lambdaR + 2 * muR * ny2,
                     2 * muR * nxy, -nx * cpR, -ny * cpR};
  PetscReal r3[5] = {-2 * muL * nxy, 2 * muL * nxy, muL * (nx2 - ny2),
                     -ny * csL, nx * csL};
  PetscReal r4[5] = {-2 * muR * nxy, 2 * muR * nxy, muR * (nx2 - ny2), ny * csR,
                     -nx * csR};

  // Compute the 4 alphas
  detP = cpR * bulkL + cpL * bulkR;
  detS = csR * muL + csL * muR;

  // P wave strengths
  a1 = (cpR * (du[0] * nx2 + du[1] * ny2 + 2 * nxy * du[2]) +
         bulkR * (nx * du[3] + ny * du[4])) / detP;
  a2 = (cpL * (du[0] * nx2 + du[1] * ny2 + 2 * nxy * du[2]) -
         bulkL * (nx * du[3] + ny * du[4])) / detP;
  // S wave strengths
  a3 = (csR * (du[2] * (nx2 - ny2) + nxy * (du[1] - du[0])) +
         muR * (nx * du[4] - ny * du[3])) / detS;
  a4 = (csL * (du[2] * (nx2 - ny2) + nxy * (du[1] - du[0])) -
         muL * (nx * du[4] - ny * du[3])) / detS;

// Compute the waves
  PetscReal W1[5], W2[5], W3[5], W4[5];
  for (int i = 0; i < 5; ++i) {
      W1[i] = a1 * r1[i];
      W2[i] = a2 * r2[i];
      W3[i] = a3 * r3[i];
      W4[i] = a4 * r4[i];
  }
  // First wave (P-wave, left-going)
  // Second wave (P-wave, right-going)
  // Third wave (S-wave, left-going)
  // Fourth wave (S-wave, right-going)
  // Use the waves to update the flux
  flux[0] = s[1] * W2[0] + s[3] * W4[0] + s[0] * W1[0] + s[2] * W3[0];
  flux[1] = s[1] * W2[1] + s[3] * W4[1] + s[0] * W1[1] + s[2] * W3[1];
  flux[2] = s[1] * W2[2] + s[3] * W4[2] + s[0] * W1[2] + s[2] * W3[2];
  flux[3] = s[1] * W2[3] + s[3] * W4[3] + s[0] * W1[3] + s[2] * W3[3];
  flux[4] = s[1] * W2[4] + s[3] * W4[4] + s[0] * W1[4] + s[2] * W3[4];
  flux[5] =0;
  flux[6] =0;
  flux[7] =0;
}