Conversation
fahnab666
changed the title
42854d4 to
a492966
Compare
a492966 to
65fb100
Compare
65fb100 to
010c3c7
Compare
Resolve material files at the case-loading boundary using the EOS family registry. Preserve the existing runtime coefficients, provenance checks, Q metadata, and qv namelist input. Stage referenced files for convergence runs and demonstrate all three analytic families with synthetic examples.
010c3c7 to
3ed3031
Compare
A reactive burn now takes Q from the products' material file instead of rejecting it. The loader sets the product qv to zero and the reactant qv to Q - Pi_r(rho0)/rho0, placing the unreacted explosive at energy Q on the products' JWL scale, so the CJ state follows from the fit. The validator and the loader share one registry helper for a family's coefficients. The new 1D_jwl_detonation example runs at 7422 m/s against the fit's CJ speed of 7404 m/s. The public LLNL handbook fits, loaded from outside the repository, reproduce TNT at 6944 m/s (6930 published) and PETN at 8363 m/s (8300 published).
Leave unknown and missing EOS parameters and a nonpositive Q to the schema and case validator, and merge the loader's parsing error paths. Refuse a reactant material file's qv alongside JWL Q instead of overwriting it, and evaluate an integer-coded reactant eos for its Pi. Keep the loader importable on Python 3.9. Trim the case.md section, keep the converted Mie-Gruneisen examples equivalent to master, and add the golden for the 1D JWL detonation example.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #1923 +/- ##
=======================================
Coverage 61.07% 61.07%
=======================================
Files 86 86
Lines 22662 22662
Branches 3342 3342
=======================================
Hits 13841 13841
Misses 6351 6351
Partials 2470 2470 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
sbryngelson
left a comment
sbryngelson
left a comment
There was a problem hiding this comment.
Thanks, this is a well-shaped feature. Expanding the file into ordinary parameters before validation keeps the solver untouched and costs nothing at run time, and the physics is right. I checked out the branch, ran test_materials.py (14 pass) and the validator/EOS tests (156 pass), and probed the loader with about 17 malformed or edge-case inputs. Inline comments carry the specifics; this is the summary.
Physics check. I recomputed the CJ state from the JWL Hugoniot using the energy form MFC implements, e = e_ref + (p - p_ref)/(omega rho), with the unreacted state (rho0, p = 0) at e0 = Q:
| Fit | CJ speed (mine) | PR | Handbook | derived C vs handbook |
|---|---|---|---|---|
| synthetic | 7403.9 m/s | 7404 | ||
| TNT | 6929.4 m/s | 6929 | 6930 | 0.010452 vs 0.01045 |
| PETN | 8300.8 m/s | 8301 | 8300 | 0.006987 vs 0.00699 |
So the reactant qv = Q - Pi_r(rho0)/rho0, product qv = 0 mapping is correct and C is indeed not an input. The golden's initial energy is also consistent (0.55e9 + 9e8 + 1600(4e6 - 9e8/1600) = 8.9e9). One request: the commit message quotes simulated speeds (6944, 8363) and the PR body quotes analytic ones (6929, 8301); please say which is which.
Must fix
- A reactive-burn products file with
Qbut norho0crashes with a bareKeyError(inline,materials.py:81). - Error messages lose the file: an unknown or missing coefficient is reported under a
fluid_pp(i)%jwl_*name the user never wrote, and seven different layout failures share one message (inline). - The reactant density is silently assumed to equal the products'
jwl_rho0(inline).
Should fix
- The layout check rejects any extra metadata (
doi,notes,units), which is too strict for a provenance format;release_statusis checked for nonemptiness and otherwise unused. - The converted convergence examples now hold
rho0/c0/gruneisen/sin bothcase.pyand the YAML (inline). - The convergence harness copies every
*.yamlin the example directory (inline). - The
.nantest does not exerciseisfinite(inline).
What would make this excellent
- Generate a per-family JSON schema for the material file from
EOS_FAMILIES(required, optional,cv/qv/qvp,Qfor JWL) and validate withfastjsonschema, which is already a dependency. That replaces_LAYOUTand the catch-alltry, gives precise file-located errors, and can be published for editor completion. - Units. A handbook JWL file is in Mbar and g/cc; MFC cases are usually SI. A provenance format without a
unitsfield is a trap. Require one and convert or fail on mismatch, or at minimum state incase.mdthat no conversion happens. - Keep provenance past expansion: write the resolved path, sha256 and citation next to the
.inpfiles (or print them on load). Today the citation is checked and then dropped, so a run cannot say which file produced its coefficients. - Turn the CJ claim into a test: a ~15-line
jwl_cj_state(...)ineos.pywith a unit test pinning the TNT/PETN handbook D_CJ and C, plus a convergence-suite check that1D_jwl_detonationruns at the CJ speed within ~1% (1500 cells, ~9000 steps, seconds). The 50-step golden only sees the initiator region; the three probes record nothing yet. - Allow stiffened-gas and ideal-gas material files, so the detonation example's reactant can use one too.
material_fileis not in the parameter registry, so./mfc.sh paramsand the generated tables do not know it exists. Registering it needs a toolchain-only (not namelist) flag.
Minor: the golden metadata cites commit 075932c, which is not in the PR history (harmless); the contribution-policy checkbox is unticked.
| if q is not None and params.get("reactive_burn", "F") == "T": | ||
| if phase != "2": | ||
| raise MFCException(f"{directive} contains Q, which belongs to the products (fluid 2) of a reactive burn") | ||
| products_q = (q, coefficients["jwl_rho0"]) |
There was a problem hiding this comment.
If the products file gives Q but omits rho0, this raises a bare KeyError: 'jwl_rho0' traceback rather than an MFCException. Reproduced with parameters: {A: 6.0, B: 0.15, R1: 4.0, R2: 1.0, omega: 0.3, Q: 2.0} and reactive_burn = T. It also fails when the case itself supplies fluid_pp(2)%jwl_rho0, which is otherwise a valid combination. Suggest reading rho0 from resolved after the merge, and raising a clear error if it is absent.
|
|
||
|
|
||
| def _reactant_qv(params: dict, q: float, rho0: float) -> float: | ||
| """Reactant qv that puts the unreacted state (rho0, p = 0) at energy Q on the products' JWL scale.""" |
There was a problem hiding this comment.
This evaluates the reactant's Pi at the products' jwl_rho0, which assumes the unreacted explosive sits at the fit's reference density. That holds for handbook fits, but if a patch initializes the reactant at another density (alpha_rho(1)/alpha(1) != jwl_rho0, e.g. a different pressing density), the derived qv is inconsistent with the initial state and nothing warns. Please either check the reactant patch densities against jwl_rho0 (the validator already walks the patches in _check_initial_states_inside_eos) or document the assumption in case.md.
Related: the third commit says a nonpositive Q is left to the schema and case validator, but Q is popped on line 54 and never reaches either. A negative Q is only caught indirectly by the qv1 > qv2 rule, and with a JWL reactant (negative Pi) it can pass. A direct Q > 0 check here would be clearer.
| ) from exc | ||
| if len(params) != len(entry["parameters"]): | ||
| raise MFCException(f"Material file '{path}' names a parameter twice") | ||
| # Names the family lacks, and Q outside JWL, reach the schema as unknown parameters. |
There was a problem hiding this comment.
Leaving unknown and missing names to the schema produces errors that do not mention the file. A c0 typo in a JWL file reports data must not contain {'fluid_pp(1)%jwl_c0'} properties, and a missing R2 reports fluid_pp(1)%eos = 'jwl' requires fluid_pp(1)%jwl_{a, b, r1, r2, omega, rho0}. The user wrote neither name. EOS_FAMILIES already carries required and optional, so checking here is about one line and lets the message name the file and the key as written.
| if {key: set(value) for key, value in data.items()} != _LAYOUT or set(entry) != {"family", "parameters"} or not all(isinstance(t, str) and t.strip() for t in texts): | ||
| raise ValueError | ||
| params = {str(name).lower(): _number(value) for name, value in entry["parameters"].items()} | ||
| except (OSError, UnicodeError, yaml.YAMLError, TypeError, KeyError, ValueError, AttributeError) as exc: |
There was a problem hiding this comment.
One message covers seven failure modes, so it often describes a file that already has everything listed. For example, adding provenance.doi or material.units fails with "requires material.name, material.eos.family, ...". Suggest separate messages (unexpected key, missing key, non-numeric value with its name), or a JSON schema via fastjsonschema generated from EOS_FAMILIES, which gives precise paths for free.
|
|
||
| _MATERIAL_KEY = re.compile(r"^fluid_pp\(([1-9][0-9]*)\)%material_file$") | ||
| _EOS = {family.suffix: family for family in EOS_FAMILIES if family.state_dependent} | ||
| _LAYOUT = {"material": {"name", "eos"}, "provenance": {"citation", "release_status"}} |
There was a problem hiding this comment.
Requiring exactly these keys makes the format reject the metadata a provenance file naturally carries (doi, source, notes, and above all units). Suggest requiring the minimum and allowing extra keys. Also, release_status is checked only for nonemptiness; either give it an enum with a behaviour (e.g. warn on anything but public) or drop it.
| matching the Fortran `case default`; a missing required one raises TypeError.""" | ||
| optional = {suffix for suffix, _math in family.optional} | ||
| values = {suffix: get(f"fluid_pp({i})%{family.prefix}_{suffix}") for suffix in family.coefficients_args} | ||
| return globals()[family.coefficients_fn](rho, *((value or 0.0) if suffix in optional else value for suffix, value in values.items())) |
There was a problem hiding this comment.
Moving this out of the validator to share it is good. The one-line dict plus generator is harder to read than the validator's old explicit loop, and globals()[...] reads oddly; the previous loop body with getattr(sys.modules[__name__], family.coefficients_fn) (or just the loop) would be clearer.
| ({"fluid_pp(2)%qv": 2.0}, {**_JWL, "qv": 1.0}, "public", "conflicts with fluid_pp\\(2\\)%qv"), | ||
| ({"fluid_pp(2)%eos": "ideal_gas"}, _JWL, "public", "conflicts with fluid_pp\\(2\\)%eos"), | ||
| ({}, _JWL, "", "nonempty provenance citation"), | ||
| ({}, {**_JWL, "A": ".nan"}, "public", "finite numeric"), |
There was a problem hiding this comment.
safe_dump writes ".nan" as a quoted string, so this hits the float() parse failure, not the isfinite check. Using float("nan") (dumped as .nan) tests the intended path; an inf case would be worth adding too.
| args = parser.parse_args() | ||
|
|
||
| rho0, p0, c0, s, gruneisen = 1.0, 1.0, 1.0, 1.5, 0.4 | ||
| rho0, p0, c0, gruneisen = 1.0, 1.0, 1.0, 0.4 |
There was a problem hiding this comment.
rho0, c0 and gruneisen now live both here (for c, dt and the perturbation amplitude) and in mg_material.yaml. Editing the YAML silently desyncs the CFL and the reference speed. The same applies to 1D_mg_impact (rho0, c0, s) and 1D_isentropic_release (rho0). Since case files run under the toolchain venv, yaml is available: reading the values from the file here would keep one source of truth.
| "probe_wrt": "T", | ||
| "fd_order": 1, | ||
| "num_probes": 3, | ||
| **{f"probe({i + 1})%x": x for i, x in enumerate([0.02, 0.03, 0.04])}, |
There was a problem hiding this comment.
In the 50-step Example test the front is still in the initiator region, so these probes record only the ambient state. The 7422 vs 7404 m/s agreement is the valuable check and is currently only in the PR description; a convergence-suite runner asserting the front speed within ~1% of the CJ value would make it permanent and is cheap at this size.
|
|
||
| #### External analytic EOS material files | ||
|
|
||
| `fluid_pp(i)%%material_file` loads analytic EOS parameters from YAML, searching the given path, the case directory, then `MFC_PUBLIC_MATERIAL_DIR`. Keep calibrated files outside the repository. |
There was a problem hiding this comment.
Please state units explicitly: the loader does no conversion, so coefficients must already be in the case's unit system. Handbook JWL fits are in Mbar, g/cc and cm/us, so this is an easy trap. Longer term, a required units field would let the loader check or convert.
|
please... Contribution Policy
[ ] I confirm this PR meets the contribution expectations and reflects my own understanding and real-world context. |
|
closed as AI contribution w/o human understanding. at least one human in the world needs to understand it |
Summary
MFC already accepts JWL, Mie–Grüneisen, and Vinet coefficients as runtime fluid parameters. This adds
fluid_pp(i)%material_file, a YAML file that keeps those coefficients together with their provenance. The toolchain expands it into the existing parameters before validation and namelist generation, so the solver is unchanged.The loader searches the given path, the case directory, then
MFC_PUBLIC_MATERIAL_DIR. It checks the file layout, a nonempty provenance, finite numeric values, and conflicts with parameters set in the case. Unknown and missing coefficients are left to the existing schema and case validator.cvandqvmap to the per-fluid parameters for all three families.A JWL products file may give
Q, the fit's detonation energy E0/rho0. In a reactive burn MFC then sets the productqvto zero and the reactantqvtoQ - Pi_r(rho0)/rho0, which places the unreacted explosive at energyQon the products' JWL scale. With the published TNT and PETN fits this reproduces the handbook CJ states (6929 vs 6930 m/s, 8301 vs 8300 m/s) and the handbook isentrope constantC, which therefore is not an input. A blast started from products does not useQ.The JWL/Vinet release and Mie–Grüneisen acoustic/impact examples now read synthetic material files, and the new
1D_jwl_detonationexample runs a stiffened-gas reactant into synthetic JWL products whoseqvvalues come fromQ. No calibrated data is added. Related to #1638.Testing
qvderived fromQfor stiffened-gas and JWL reactants.1D_jwl_detonationhas a new golden (DAC932D4). At full resolution its front runs at 7422 m/s against the fit's CJ speed of 7404 m/s../mfc.sh validateon every example pass locally. The SLURM monitor and submission unit tests were excluded because they depend on GNU tools.Contribution Policy