Convergence study of the inviscid part of MFC with multicomponent species

[RadouanBoukharfane]

Hello everyone,

Has anyone performed a convergence study of the inviscid part of the MFC solver with multicomponent species?

I ran the 1D_titarevtorro_analytical case using the default WENO5 reconstruction. Based on the scheme, I was expecting 5th-order convergence, but my results show only about 1st order:

     N         h≈dx                 L2
-----------------------------------------------------------
    27   7.407407e-02   1.0759333255e-03
    52   3.846154e-02   5.5889558236e-04
   102   1.960784e-02   2.8496073423e-04
   202   9.900990e-03   1.4389545890e-04
   402   4.975124e-03   7.2311041434e-05

Has anyone encountered this before, or can confirm the expected convergence rate in this setup?
Any hints on whether this is a known limitation of the current implementation, or if additional options/settings are needed to recover the theoretical WENO5 order would be very helpful.

Thanks in advance!

Best,
Radouan

---

Here are my case.py parameters:

#!/usr/bin/env python3
import math
import json
import cantera as ct
import argparse

parser = argparse.ArgumentParser(prog="1D_titarevtorro_shocktube", formatter_class=argparse.ArgumentDefaultsHelpFormatter)

parser.add_argument(
    "--mfc",
    type=json.loads,
    default="{}",
    metavar="DICT",
    help="MFC's toolchain's internal state.",
)
parser.add_argument(
    "--no-chem",
    dest="chemistry",
    default=True,
    action="store_false",
    help="Disable chemistry.",
)

args = parser.parse_args()

ctfile = "h2o2.yaml"
solT = ct.Solution(ctfile)
solT.TPX = 300, 8000, "O2:0.21, N2:0.79"

# Numerical setup
xmin = 0.0
xmax = 2.0
Nx = 25
Lx = xmax - xmin
dx = Lx / Nx
Tend =  2

# Configuring case dictionary
case = {
    "run_time_info": "T",
    "x_domain%beg": xmin,
    "x_domain%end": xmax,
    "m": Nx+1,
    "n": 0,
    "p": 0,
    "cfl_adap_dt": "T",
    "cfl_target": 0.5,
    "n_start": 0,
    "t_stop": Tend,
    "t_save": Tend/10.0,
    "num_patches": 1,
    "model_eqns": 2,
    "num_fluids": 1,
    "mpp_lim": "F",
    "mixture_err": "F",
    "time_stepper": 3,
    "weno_order": 5,
    "weno_eps": 1.0e-16,
    "mapped_weno": "T",
    "null_weights": "F",
    "mp_weno": "F",
    "riemann_solver": 2,
    "wave_speeds": 1,
    "avg_state": 2,
    "bc_x%beg": -1,
    "bc_x%end": -1,
    "format": 2,
    "precision": 2,
    "prim_vars_wrt": "T",
    "parallel_io": "F",
    "gamma_wrt" : "T", 
    "chemistry": "F" if not args.chemistry else "T",
    "chem_params%diffusion": "F",
    "chem_params%reactions": "F",
    "patch_icpp(1)%geometry": 1,
    "patch_icpp(1)%x_centroid": Lx/2,
    "patch_icpp(1)%length_x": Lx,
    "patch_icpp(1)%vel(1)": 0.05,
    "patch_icpp(1)%pres": solT.P,
    "patch_icpp(1)%alpha_rho(1)": "0.09253158446106331*(1 + 0.1*sin(x*pi))",
    "patch_icpp(1)%alpha(1)": 1.0,
    "fluid_pp(1)%gamma": 1.0e00 / (1.55e00 - 1.0e00),
    "fluid_pp(1)%pi_inf": 0.0,
    "cantera_file": ctfile,
}

if args.chemistry:
    for i in range(len(solT.Y)):
        case[f"patch_icpp(1)%Y({i+1})"] = solT.Y[i]

if __name__ == "__main__":
    print(json.dumps(case))

---

[sbryngelson]

Hi @RadouanBoukharfane -- thanks for the question and investigating MFC!

Would you not expect the seemingly seen first-order accuracy for a problem with a shock? You can confirm the design order of accuracy on a smooth problem, of course.

We have multicomponent chemistry cases that have discontinuities, though they have not seen a grid study yet. But I would not expect them to converge at order 5, for example.

Useful ref: https://www3.nd.edu/~powers/paper.list/jcp.weno5m.05.pdf

[RadouanBoukharfane]

Hi @sbryngelson and thanks for the quick reply.

The test case I’m running is smooth -it’s purely sinusoidal (no shocks or discontinuities). That’s why I was expecting to observe the design-order convergence rather than the first-order behavior I’m seeing.

[sbryngelson]

Ah yes I should have seen this. It's not the usual TT problem. Ok, do you see 5th order convergence (or something >1) if you take out chemistry? What about single-species?

[RadouanBoukharfane]

That’s exactly the “magic” I can’t quite explain yet. For purely aero cases (no chemistry), I do recover the expected high-order convergence of WENO. But as soon as I couple with Cantera and the thermodynamic/transport quantities depend on temperature, the observed order collapses back to first order.

[sbryngelson]

This is a good finding. I don't have a definitive answer at the moment. I'm curious what your cantera file looks like, and what happens if you drop to single-species chemistry? cc'ing @DimAdam-01, @henryleberre, @ecisneros8

[RadouanBoukharfane]

I’m currently using ``h2o2.yaml'' with 9 species, and for the setup I’ve run so far I used a mixture of 0.21 O₂ and 0.79 N₂. I will see as sugggested the single-species case you mention, to see how the behavior compares.

[sbryngelson]

Got it. Keep us posted and I think we will hear from some of the cc'ed above as well.

[ecisneros8]

Thanks for looking into this @RadouanBoukharfane.

Can you clarify whether your case has temperature variations? I vaguely remember we went back and forth on how to compute the specific heat ratio for multi-component mixtures with variable T in the Riemann solver. I can see that part of the code as a potential source of error. @DimAdam-01 should be able to tell us more about what he did there.

The 1 species case is definitely the place to start. (There are other things you can test too, like modifying the NASA polynomials in your Cantera input file to yield constant thermodynamic properties. But, if you do continue to run multi-species, start with only O2 and N2 in your mechanism).

Keep us posted!

[DimAdam-01]

Thanks for your question @RadouanBoukharfane.

Are you comparing against an analytical solution or a solution from a very dense grid?

That's the key question. The method has to change from using the analytical solution to using a dense grid.

The analytical solution is only valid for a simplified problem with constant thermodynamic properties. Once chemistry is enabled with Cantera, those properties become temperature-dependent, which fundamentally changes the problem. This makes the analytical solution the wrong benchmark.

The proper procedure is a grid refinement study:

A benchmark solution is generated by running the full chemistry simulation on the finest possible grid (e.g., 4096 cells).

The L2 error for all coarser grids is then calculated by comparing them against this single benchmark solution.