diff --git a/README.md b/README.md index 611164b..731b5e7 100644 --- a/README.md +++ b/README.md @@ -161,7 +161,9 @@ options: --translate TRANSLATE TRANSLATE Translate the (x, y) coordinates by this amount --mapaxes Use the MAPAXES coordinates for an Eclipse file - --pinch Remove pinched elements + --pinch Remove pinched elements (coincident corners within + --pinch-tol). Off by default; faithful conversion is + produced without this flag. --pinch-tol PINCH_TOL Tolerance for coincident corners when removing pinched elements (default: 1e-3) @@ -202,6 +204,10 @@ options: --strict-jacobians Treat any non-positive element Jacobian as a fatal error and exit non-zero. By default such elements only produce a warning. Useful for CI / scripted workflows. + --remove-distorted Remove elements with non-positive Jacobians (degenerate + or inverted) from the output mesh, reporting a count of + those removed. By default such elements are kept and + only a warning is printed. ``` ### Lateral refinement (Eclipse only) @@ -311,32 +317,43 @@ A few practical notes: ### Element Jacobian check -After the mesh is built, em2ex automatically computes the Jacobian at all 8 corners of every HEX8 element and prints a one-line summary: +After conversion, em2ex evaluates the Jacobian at all 8 corners of every HEX8 element and prints a one-line summary: ``` Element Jacobian check: 1000000 / 1000000 elements OK ``` -A non-positive Jacobian (negative = inverted element, zero = degenerate) almost always means one of two things: +A non-positive Jacobian (zero = degenerate, negative = inverted) almost always means one of two things: - **The input file uses a `z increases upward` convention** rather than Eclipse's default `z increases downward`. The cells come out "upside down" and need `--flip` to invert the z values. -- **There's something genuinely wrong with the input geometry** — corner ordering inconsistent within a cell, self-intersecting cells that survived pinch removal, etc. +- **There are a few genuinely distorted cells** — for example near faults that are so skewed the element is inverted, or near-pinched cells that slipped through the tolerance filter. -The relevance is downstream: most finite-element solvers (MOOSE, libMesh, MFEM, deal.II) reject elements with non-positive Jacobian at the assembly stage. Catching the problem at conversion time is cheaper than chasing it through a failed simulation. +When all elements are non-positive the output is expanded with a hint: -When the check fails, em2ex prints an expanded report with element IDs and centroid locations: +``` +Element Jacobian check: 27 negative, 0 zero, 0 OK (out of 27) + All elements have non-positive Jacobians. Possible causes: + - z-up coordinate convention: try --flip to invert the z-axis + - Orientation-reversing coordinate system (e.g. MAPAXES handedness) + Use --remove-distorted to remove these elements and proceed anyway. +``` + +When only some elements are bad, the report includes element IDs and centroid locations: ``` -Element Jacobian check: 1 negative, 0 zero, 0 OK (out of 1) - Examples of negative-Jacobian elements (showing up to 5 of 1): - element 1: centroid (0.5, 0.5, 0.5), min Jacobian = -1.000e+00 +Element Jacobian check: 3 negative, 0 zero, 24 OK (out of 27) + Examples of negative-Jacobian elements (showing up to 5 of 3): + element 7: centroid (0.5, 0.5, 0.5), min Jacobian = -1.000e-03 ``` -By default the conversion continues despite warnings (the Exodus file is still written; the user is informed). Two flags adjust this: +By default the output file is still written regardless of warnings. Three flags adjust this: -- `--strict-jacobians` upgrades any non-positive Jacobian to a fatal error (exit code 1). Useful in CI / scripted workflows where a bad mesh should stop the pipeline rather than propagate downstream. +- `--remove-distorted` removes elements with non-positive Jacobians before writing, printing a count of those dropped. Useful when a small number of distorted cells near faults would otherwise cause solver failures. +- `--strict-jacobians` upgrades any non-positive Jacobian to a fatal error (exit code 1). Useful in CI / scripted workflows where a bad mesh should stop the pipeline. - `--no-check-jacobians` skips the check entirely (a small time saving on very large grids, but disables a useful safety net). +**Relationship to `--pinch`.** `--remove-distorted` and `--pinch` are complementary rather than interchangeable. `--pinch` detects cells where any two corners are within `--pinch-tol` of each other — it catches *near*-coincident corners (e.g. a 0.5 m thick cell in a grid measured in metres) even when the Jacobian is still technically positive. `--remove-distorted` catches cells whose Jacobian has already reached zero or gone negative, which only happens once corners are *exactly* coincident or the cell has become inverted. In practice, running both is the safest option for grids with thin reservoir layers near faults: `--pinch` handles the near-zero cells that `--remove-distorted` would miss, and `--remove-distorted` catches any remaining inverted cells. + `em2ex` attempts to guess the reservoir model format from the file extension (see supported formats below). If the reservoir model has a non-standard file extension, the user can force `em2ex` to read the correct format using the `--filetype` commandline option. diff --git a/em2ex.py b/em2ex.py index 984d8d3..a1df4a1 100755 --- a/em2ex.py +++ b/em2ex.py @@ -117,7 +117,7 @@ def get_parser(): parser.add_argument('--flip', dest = 'flip_z', action = 'store_true', help = 'Flip the sign of the Z coordinates') parser.add_argument('--translate', nargs = 2, dest = 'translate', type = float, help = 'Translate the (x, y) coordinates by this amount') parser.add_argument('--mapaxes', dest = 'use_mapaxes', action = 'store_true', help = 'Use the MAPAXES coordinates for an Eclipse file') - parser.add_argument('--pinch', default = True, dest = 'no_pinch', action = 'store_true', help = 'Remove pinched elements') + parser.add_argument('--pinch', default = False, dest = 'no_pinch', action = 'store_true', help = 'Remove pinched elements (coincident corners within --pinch-tol). Off by default; faithful conversion is produced without this flag.') parser.add_argument('--pinch-tol', default = 1e-3, dest = 'pinch_tol', type = float, help = 'Tolerance for coincident corners when removing pinched elements (default: 1e-3)') parser.add_argument('--refine-xy', nargs = 2, dest = 'refine_xy', type = _positive_int, metavar = ('RX', 'RY'), @@ -142,6 +142,8 @@ def get_parser(): help = 'Skip the per-element Jacobian sanity check. By default em2ex computes the Jacobian at all 8 corners of every HEX8 element and warns if any are non-positive (degenerate or inverted), since those elements would be rejected by most FEM solvers.') parser.add_argument('--strict-jacobians', dest = 'strict_jacobians', action = 'store_true', help = 'Treat any non-positive element Jacobian as a fatal error and exit non-zero. By default such elements only produce a warning. Useful for CI / scripted workflows.') + parser.add_argument('--remove-distorted', dest = 'remove_distorted', action = 'store_true', + help = 'Remove elements with non-positive Jacobians (degenerate or inverted) from the output mesh, reporting a count of those removed. By default such elements are kept and only a warning is printed.') return parser def main(): diff --git a/readers/eclipse.py b/readers/eclipse.py index b57aec6..5ab5c1a 100644 --- a/readers/eclipse.py +++ b/readers/eclipse.py @@ -367,6 +367,14 @@ def parseEclipse(f, args): xvec = xvec / np.sqrt(xvec[0]**2 + xvec[1]**2) yvec = yvec / np.sqrt(yvec[0]**2 + yvec[1]**2) + # If the MAPAXES transform is orientation-reversing (det < 0), negate xvec so + # the resulting mesh has positive Jacobians instead of all-negative ones. + det = xvec[0] * yvec[1] - xvec[1] * yvec[0] + if det < 0: + print("MAPAXES transform is orientation-reversing (det={:.6g}); " + "negating x-axis direction to restore positive Jacobians".format(det)) + xvec = -xvec + # Transform top and bottom pillar x,y: world = origin + local_x*xhat + local_y*yhat for xi, yi in [(0, 1), (3, 4)]: xdata = coord[:,:,xi].copy() @@ -404,12 +412,40 @@ def parseEclipse(f, args): # All elements are active active_elements = np.ones((nz, ny, nx), dtype = int) - # If pinched or distorted elements are removed, these elements will also be set to - # be inactive. - # Check for pinched out elements (where the vertical distance between nodes is less than tol) + # Check for pinched elements (coincident corners within pinch_tol). + # Always detect so a count can be reported; only remove when --pinch is passed. + distorted = distortedElem(elemcornx, elemcorny, elemcornz, args.pinch_tol).reshape(nz, ny, nx) + n_pinched = int((distorted & (active_elements > 0)).sum()) if args.no_pinch: - distorted = distortedElem(elemcornx, elemcorny, elemcornz, args.pinch_tol).reshape(nz, ny, nx) + if n_pinched > 0: + print('{} pinched element(s) removed'.format(n_pinched)) active_elements[distorted] = 0 + elif n_pinched > 0: + print('Note: {} element(s) have coincident corners (pinched) and will be ' + 'invalid in FEM solvers. Use --pinch to remove them.'.format(n_pinched)) + + # Optionally remove elements with non-positive Jacobian (degenerate or inverted). + # Off by default (faithful conversion); enabled with --remove-distorted. + if getattr(args, 'remove_distorted', False): + from readers.reader_utils import nonPositiveJacobianElems + # When z is flipped, the top/bottom corner swap applied later to elemNodes + # (corners [4,5,6,7,0,1,2,3]) restores positive Jacobians. Apply the same + # permutation here so the check reflects the final assembled element orientation. + cx, cy, cz = elemcornx, elemcorny, elemcornz + if args.flip_z: + perm = [4, 5, 6, 7, 0, 1, 2, 3] + cx, cy, cz = cx[:, perm], cy[:, perm], cz[:, perm] + bad_jac = nonPositiveJacobianElems(cx, cy, cz).reshape(nz, ny, nx) + newly_removed = bad_jac & (active_elements > 0) + n_removed = int(newly_removed.sum()) + if n_removed > 0: + n_active = int((active_elements > 0).sum()) + if n_removed == n_active: + print('Warning: --remove-distorted would remove all {} active element(s). ' + 'Skipping removal — check grid orientation (try --flip).'.format(n_active)) + else: + print('{} element(s) with zero or negative element Jacobian removed'.format(n_removed)) + active_elements[bad_jac] = 0 # The number of active elements is num_active_elements = np.count_nonzero(active_elements) diff --git a/readers/reader_utils.py b/readers/reader_utils.py index e1d99dd..91ddce8 100644 --- a/readers/reader_utils.py +++ b/readers/reader_utils.py @@ -154,6 +154,12 @@ def checkElementJacobians(model, strict=False): print('Element Jacobian check: {} negative, {} zero, {} OK (out of {})'.format( num_neg, num_zero, num_ok, model.numElems)) + if num_ok == 0: + print(' All elements have non-positive Jacobians. Possible causes:') + print(' - z-up coordinate convention: try --flip to invert the z-axis') + print(' - Orientation-reversing coordinate system (e.g. MAPAXES handedness)') + print(' Use --remove-distorted to remove these elements and proceed anyway.') + # Map elemNodes row -> Exodus element ID. elemIds[k, j, i] holds the # 1-based Exodus ID (0 for inactive); the non-zero values in (k, j, i) # flat order align with elemNodes' row ordering. @@ -182,6 +188,23 @@ def checkElementJacobians(model, strict=False): return False +def nonPositiveJacobianElems(elemcornx, elemcorny, elemcornz): + ''' Return a boolean array (numelems,) that is True for elements whose + minimum per-corner Jacobian is <= 0 (degenerate or inverted). + + Inputs are per-element corner coordinate arrays of shape (numelems, 8), + ordered in the HEX8 element corner layout (matching elemCornerCoords output). + ''' + P = np.stack([elemcornx, elemcorny, elemcornz], axis=-1) # (N, 8, 3) + jacobians = np.empty((P.shape[0], 8)) + for c, ((xi_a, xi_b), (eta_a, eta_b), (zeta_a, zeta_b)) in enumerate(_HEX8_JAC_EDGES): + e_xi = P[:, xi_a] - P[:, xi_b] + e_eta = P[:, eta_a] - P[:, eta_b] + e_zeta = P[:, zeta_a] - P[:, zeta_b] + jacobians[:, c] = np.einsum('ij,ij->i', e_xi, np.cross(e_eta, e_zeta)) + return jacobians.min(axis=1) <= 0 + + def nonZeroValues(arr): ''' Utility to determine first non-zero values in the top plane of an array ''' diff --git a/run_tests.py b/run_tests.py index dd5048d..0a520fd 100755 --- a/run_tests.py +++ b/run_tests.py @@ -66,6 +66,13 @@ def test_em2ex(key, use_official_api, exodiff): assert tests[key]['expected_error'] in str(excinfo.value) + # If the test type is output, check that expected text appears in stdout + elif tests[key]['type'] == 'output': + if 'expected_output' not in tests[key].keys(): + pytest.skip(key + ': Skipped as expected_output not specified') + else: + check_output(key) + else: # Skip unknown test type pytest.skip(key + ': Skipped as unknown test type') @@ -125,3 +132,24 @@ def expected_error(key): raise Em2exException(result.stdout + result.stderr) return + +def check_output(key): + ''' Check that expected text appears in em2ex stdout on a successful run ''' + + filepath = tests[key]['filepath'] + filename = tests[key]['filename'] + testfilename = os.path.join(filepath, filename) + + command = ['./em2ex.py', '-f'] + if 'cli_args' in tests[key].keys(): + command.extend(tests[key]['cli_args'].split()) + command.append(testfilename) + + result = subprocess.run(command, capture_output=True, text=True) + combined = result.stdout + result.stderr + + assert tests[key]['expected_output'] in combined, \ + '{}: expected output not found.\nExpected: {}\nGot: {}'.format( + key, tests[key]['expected_output'], combined) + + return diff --git a/test/eclipse/gold/simple_cube_mapaxes.e b/test/eclipse/gold/simple_cube_mapaxes.e index 88211e6..c5b1c0a 100644 Binary files a/test/eclipse/gold/simple_cube_mapaxes.e and b/test/eclipse/gold/simple_cube_mapaxes.e differ diff --git a/test/eclipse/tests b/test/eclipse/tests index 2565a9e..c006c01 100644 --- a/test/eclipse/tests +++ b/test/eclipse/tests @@ -38,6 +38,8 @@ simple_cube_translated: simple_cube_mapaxes: filename: simple_cube_mapaxes.grdecl type: exodiff + # The MAPAXES spec (1.0 0.0 0.0 0.0 0.0 1.0) is orientation-reversing (det=-1). + # The reader detects this and negates xvec so the output has positive Jacobians. cli_args: --mapaxes gold: simple_cube_mapaxes.e @@ -151,15 +153,22 @@ simple_cube_quoted_unit: type: exodiff gold: simple_cube.e -# Z-up convention file (pillar tops at z=1, bots at z=0). em2ex's default -# corner ordering treats this as inverted; --strict-jacobians should error -# out rather than write a mesh +# Z-up convention file (pillar tops at z=1, bots at z=0). em2ex produces a +# faithful conversion but the Jacobian check sees all-negative elements; +# --strict-jacobians upgrades that to a fatal error. strict_jacobians_inverted: filename: simple_cube_zup.grdecl type: exception cli_args: --strict-jacobians expected_error: exiting due to invalid Jacobians +# Same z-up file without --strict-jacobians: conversion succeeds but the +# all-negative hint (suggesting --flip) should be printed. +all_negative_jacobians_hint: + filename: simple_cube_zup.grdecl + type: output + expected_output: All elements have non-positive Jacobians + # Driven entirely by a YAML config file (extract + refine settings). The # output must match the same gold as the CLI-driven simple_cube_extract_refine # test, proving the config loader applies values correctly.