Tutorial Step-67 Extension HLLC

[Yeji Wi]

Hi all!

I am trying to implement HLLC flux in the code for tutorial step-67.

I have implemented the flux just as described in the textbook.

However, this code results 'dt' to 'inf' and not a number errors...

It looks like I am missing something but not sure what it could be as I am not insightful enough to see the issue by myself..

The only modification I have done to step-67.cc is adding this hllc case to the switch in euler_numerical_flux and adding hllc option to enum EulerNumericalFlux. Any suggestion will be helpful for me and if I am missing any crucial info that would be helpful, please let me know. I have attached the code as well. Thank you in advance!

case hllc:

{

const auto avg_vel_normal = 0.5 * ((velocity_m + velocity_p) * normal);

const auto avg_c = std::sqrt(std::abs(0.5 * gamma * (pressure_p * (1. / u_p[0]) + pressure_m * (1. / u_m[0]))));

const auto avg_p = 0.5 * (pressure_m + pressure_p);

const auto rho_l = u_m[0];

const auto rho_r = u_p[0];

const auto p_l = pressure_m;

const auto p_r = pressure_p;

const auto v_l = velocity_m * normal;

const auto v_r = velocity_p * normal;

const auto avg_rho = 0.5 * (rho_l + rho_r);

const auto avg_a = avg_c;

const auto p_pvrs = 0.5 * (p_l + p_r) - 0.5 * (v_r - v_l) * avg_rho * avg_a ;

const auto p_star = std::max(Number(), p_pvrs);

double q_L = 0.;

if (p_star.data <= p_l.data)

q_L = 1.;

else

q_L = std::sqrt((1+ (gamma + 1)*(p_star.data/p_l.data - 1)/(2 * gamma)));

double q_R = 0.;

if (p_star.data <= p_r.data)

q_R = 1.;

else

q_R = std::sqrt((1+ (gamma + 1)*(p_star.data/p_r.data - 1)/(2 * gamma)));

const auto S_L = v_l - avg_a * q_L ;

const auto S_R = v_r + avg_a * q_R ;

const auto S_num = p_r - p_l + rho_l * v_l * (S_L - v_l) - rho_r * v_r * (S_R - v_r);

const auto S_den = rho_l * (S_L - v_l) - rho_r * (S_R - v_r);

const auto S_star = S_num / S_den;

const auto P_LR = 0.5 * (p_l + p_r + rho_l*(S_L - v_l)*(S_star - v_l) + rho_r * (S_R - v_r)*(S_star - v_r));

// Compute fluxes

const auto F_L = flux_m * normal;

const auto F_R = flux_p * normal;

const double zero = 0.;

// Compute D* tensor

Tensor<1, dim + 2, Number> D_star;

D_star[0] = 0;

D_star[1] = 1;

for (unsigned int d = 0; d < dim-1; ++d) {

D_star[d + 2] = 0;

}

D_star[dim + 1] = S_star;

// Compute the numerical flux

Tensor<1, dim + 2, Number> numerical_flux;

if (zero <= S_L.data)

numerical_flux = F_L;

else if (S_L.data <= zero && zero <= S_star.data)

for (unsigned int d = 0; d < dim + 2; ++d)

numerical_flux[d] = (S_star * (S_L * u_m[d] - F_L[d]) + S_L * P_LR * D_star[d])/(S_L - S_star);

else if (S_star.data <= zero && zero <= S_R.data)

for (unsigned int d = 0; d < dim + 2; ++d)

numerical_flux[d] = (S_star * (S_R * u_p[d] - F_R[d]) + S_R * P_LR * D_star[d])/(S_R - S_star);

else if (zero >= S_R.data)

numerical_flux = F_R;

return numerical_flux;

}

enum EulerNumericalFlux

{

lax_friedrichs_modified,

harten_lax_vanleer,

hllc,

};

constexpr EulerNumericalFlux numerical_flux_type = hllc;

[Wolfgang Bangerth]

Yeji:

> I am trying to implement HLLC flux in the code for tutorial step-67.
> I have implemented the flux just as described in the textbook.
> However, this code results 'dt' to 'inf' and not a number errors...

> Screenshot 2024-04-14 at 2.06.49 PM.png
>
> It looks like I am missing something but not sure what it could be as I
> am not insightful enough to see the issue by myself..

Rather than trying to find the problem in your code, let me outline how
you find out what is wrong. Finding out why a function computes inf as
an answer is actually one of the easier tasks. That's because once you
have an infinity, it doesn't tend to go away: If a variable is infinity,
all the things you compute from it will also be infinities of some sort.
So all you have to figure out is where the infinity is first created.

Let's say you have a function that perhaps reads like this:

double hllc(double a, double b, double c) {
d = sqrt(a+b+c)
e = sin(a*d)
f = acos(d+e)
return f;
}

In a place where you call this function, the function returns inf when
you think that it should not. So for given arguments a,b,c, you now know
that f=inf. This is wrong, so we may as well abort the program at this
point. Make it so:

double hllc(double a, double b, double c) {
d = sqrt(a+b+c)
e = sin(a*d)
f = acos(d+e)
Assert (numbers::is_finite(f), ExcMessage ("f is infinite!"));
return f;
}

When you run your program, it should not abort in the Assert line
because you already know that the function returns inf.

Since you know that f was computed from d and e, using the acos
function, there are only two possibilities for d to be inf: either the
inputs to acos were already inf, or the argument to acos was out of
range (it needs to be between -1 and +1). Test this:

double hllc(double a, double b, double c) {
d = sqrt(a+b+c)
e = sin(a*d)
Assert (numbers::is_finite(d), ExcMessage ("d is infinite!"));
Assert (numbers::is_finite(e), ExcMessage ("e is infinite!"));
Assert ((d+e >= -1) && (d+e <= +1), ExcMessage ("d+e out of range!"));
f = acos(d+e)
Assert (numbers::is_finite(f), ExcMessage ("f is infinite!"));
return f;
}

Your program should now abort in one of the three new assertions.
Depending on which one, you have now gained knowledge which variable
seems to be wrong. Use this to add more assertions further up the chain.
For example

double hllc(double a, double b, double c) {
Assert (numbers::is_finite(a), ExcMessage ("a is infinite!"));
Assert (numbers::is_finite(b), ExcMessage ("b is infinite!"));
Assert (numbers::is_finite(c), ExcMessage ("c is infinite!"));
Assert (a+b+c>=0, ExcMessage ("a+b+c out of range!"));
d = sqrt(a+b+c)
Assert (numbers::is_finite(a), ExcMessage ("a is infinite!"));
Assert (numbers::is_finite(d), ExcMessage ("d is infinite!"));
e = sin(a*d)
Assert (numbers::is_finite(d), ExcMessage ("d is infinite!"));
Assert (numbers::is_finite(e), ExcMessage ("e is infinite!"));
Assert ((d+e >= -1) && (d+e <= +1), ExcMessage ("d+e out of range!"));
f = acos(d+e)
Assert (numbers::is_finite(f), ExcMessage ("f is infinite!"));
return f;
}

The point is that you work yourself bottom-up until you figure out in
which place the infinity is being created. When you know which line
creates the infinity, you can start thinking about what it is that's wrong.

Best
W.