Negatively oriented elements in oneX mesh

[matthiasdiener]

  File "/shared/home/mdiener/Work/efuse/meshmode/meshmode/mesh/__init__.py", line 1037, in __init__
    raise ValueError("negatively oriented elements found")
ValueError: negatively oriented elements found

[nchristensen]

I'm seeing this for the other 3D meshes too.

[lukeolson]

which meshes specifically?

[nchristensen]

test_mesh/eager/eigthX/data/isolator.msh
test_mesh/eager/quarterX/data/isolator.msh
test_mesh/eager/oneX/data/isolator.msh

I haven't tested the others but they likely have the same problem.

[anderson2981]

How/when were the meshes generated. Which version of gmsh? To my knowledge the meshes are not checked into the repo due to size. There was some issue that has been fixed when I switched to the OpenCascade generator inside gmsh. But that fix was made some time ago and as far as I know only pertained to 2D meshes, not 3D.

[nchristensen]

I generated them right before running yesterday using gmsh 4.9.5.

[anderson2981]

which platform? and what version of meshmode?

I've rebuilt on Lassen and I do not see this issue, meshmode is pointing at 6bfe6a, pytato-array-context-transforms

[nchristensen]

I'm running on Crusher. I just tried with that branch and I still see the problem. I can give it a try on Lassen to see if I can replicate the problem there.

[nchristensen]

I also encounter this problem on Lassen. Output below.

(dev) [njchris@lassen559:eigthX]$ mpirun -np 1 python -m mpi4py isolator.py 
Default casename isolator
2022-03-23 13:57:21,625 - INFO - grudge.array_context - get_reasonable_array_context_class: MPIPyOpenCLArrayContext lazy=False distributed=True device-parallel=True
No user input file, using default values
Running isolator.py

Shock capturing parameters: alpha 0.3, s0 -4.0, kappa 0.5

#### Simluation control data: ####
	nviz = 500
	nrestart = 5000
	nhealth = 1
	nstatus = 1
	current_dt = 1e-08
	t_final = 1e-07
	order = 1
	dimen = 2
	Time integration rk4
#### Simluation control data: ####


#### Simluation material properties: ####
	mu = 4.216360056e-05
	kappa = 0.05621788139856423
	Prandtl Number  = 0.75
#### Simluation initialization data: ####
	inlet Mach number 0.038727226466056475
	inlet temperature 2075.8073416311686
	inlet pressure 274212.0081934861
	inlet rho 0.4623484671063763
	inlet velocity 35.28896626974282
	outlet Mach number 3.661995139955526
	outlet temperature 563.9344091915701
	outlet pressure 2865.7077436508835
	outlet rho 0.01778579568979098
	outlet velocity 1739.2462693097807
#### Simluation initialization data: ####

Reading mesh from data/isolator.msh
2022-03-23 13:57:26,385 - INFO - meshmode.distributed - rank 0: sent all mesh partitions
Making discretization
Done making discretization
Before restart/init
Initializing soln.
Initialization for Case(isolator)
===
Num 2d order-1 elements: 51730
Num global elements: 51730
Timestep:        1e-08
Final time:      1e-07
CFL:             1.0
Constant CFL:    False
Initialization:  isolator
EOS:             IdealSingleGas

isolator.py:1098: DeprecationWarning: EagerDGDiscretization is deprecated and will go away in 2022. Use the base DiscretizationCollection with grudge.op instead.
  discr = EagerDGDiscretization(
/g/g15/njchris/miniforge3/envs/dev/lib/python3.8/site-packages/pytools/__init__.py:920: UserWarning: The geometric factor for the characteristic length scale in time step estimation is not necessarily valid for non-volume-filling discretizations. Continuing anyway.
  result = inner(*args)
2022-03-23 13:57:33,623 - INFO - meshmode.discretization.connection.face - building face restriction: start
Traceback (most recent call last):
  File "/g/g15/njchris/miniforge3/envs/dev/lib/python3.8/site-packages/pytools/__init__.py", line 918, in new_inner
    return self.cache_dict[args]
KeyError: ()

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/g/g15/njchris/miniforge3/envs/dev/lib/python3.8/runpy.py", line 194, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "/g/g15/njchris/miniforge3/envs/dev/lib/python3.8/runpy.py", line 87, in _run_code
    exec(code, run_globals)
  File "/g/g15/njchris/miniforge3/envs/dev/lib/python3.8/site-packages/mpi4py/__main__.py", line 7, in <module>
    main()
  File "/g/g15/njchris/miniforge3/envs/dev/lib/python3.8/site-packages/mpi4py/run.py", line 198, in main
    run_command_line(args)
  File "/g/g15/njchris/miniforge3/envs/dev/lib/python3.8/site-packages/mpi4py/run.py", line 47, in run_command_line
    run_path(sys.argv[0], run_name='__main__')
  File "/g/g15/njchris/miniforge3/envs/dev/lib/python3.8/runpy.py", line 265, in run_path
    return _run_module_code(code, init_globals, run_name,
  File "/g/g15/njchris/miniforge3/envs/dev/lib/python3.8/runpy.py", line 97, in _run_module_code
    _run_code(code, mod_globals, init_globals,
  File "/g/g15/njchris/miniforge3/envs/dev/lib/python3.8/runpy.py", line 87, in _run_code
    exec(code, run_globals)
  File "isolator.py", line 1527, in <module>
    main(restart_filename=restart_filename, user_input_file=input_file,
  File "/g/g15/njchris/Workspace/mirgecom-upstream/mirgecom/mpi.py", line 157, in wrapped_func
    func(*args, **kwargs)
  File "isolator.py", line 1435, in main
    advance_state(rhs=my_rhs, timestepper=timestepper,
  File "/g/g15/njchris/Workspace/mirgecom-upstream/mirgecom/steppers.py", line 330, in advance_state
    _advance_state_stepper_func(
  File "/g/g15/njchris/Workspace/mirgecom-upstream/mirgecom/steppers.py", line 127, in _advance_state_stepper_func
    state, dt = pre_step_callback(state=state, step=istep, t=t, dt=dt)
  File "isolator.py", line 1369, in my_pre_step
    alpha_field = my_get_alpha(discr, fluid_state, alpha_sc)
  File "isolator.py", line 1346, in my_get_alpha
    length_scales = characteristic_lengthscales(array_context, discr)
  File "/g/g15/njchris/Workspace/grudge-upstream/grudge/dt_utils.py", line 107, in characteristic_lengthscales
    return thaw(_compute_characteristic_lengthscales(), actx)
  File "/g/g15/njchris/miniforge3/envs/dev/lib/python3.8/site-packages/pytools/__init__.py", line 920, in new_inner
    result = inner(*args)
  File "/g/g15/njchris/Workspace/grudge-upstream/grudge/dt_utils.py", line 102, in _compute_characteristic_lengthscales
    thaw(dt_geometric_factors(dcoll), actx)
  File "/g/g15/njchris/Workspace/grudge-upstream/grudge/dt_utils.py", line 273, in dt_geometric_factors
    face_discr = dcoll.discr_from_dd(dd_face)
  File "/g/g15/njchris/miniforge3/envs/dev/lib/python3.8/site-packages/pytools/__init__.py", line 706, in wrapper
    result = function(obj, *args, **kwargs)
  File "/g/g15/njchris/Workspace/grudge-upstream/grudge/discretization.py", line 336, in discr_from_dd
    return self._all_faces_volume_connection().to_discr
  File "/g/g15/njchris/miniforge3/envs/dev/lib/python3.8/site-packages/pytools/__init__.py", line 706, in wrapper
    result = function(obj, *args, **kwargs)
  File "/g/g15/njchris/Workspace/grudge-upstream/grudge/discretization.py", line 653, in _all_faces_volume_connection
    return make_face_restriction(
  File "/g/g15/njchris/Workspace/meshmode-upstream/meshmode/discretization/connection/face.py", line 345, in make_face_restriction
    bdry_mesh = Mesh(bdry_vertices, bdry_mesh_groups)
  File "/g/g15/njchris/Workspace/meshmode-upstream/meshmode/mesh/__init__.py", line 1037, in __init__
    raise ValueError("negatively oriented elements found")
ValueError: negatively oriented elements found

[MTCam]

I also encounter this problem on Lassen. Output below.
[omitted error output]

I concur. The errors happen with fresh checkouts and builds of mirgecom on Lassen. I've just tested it this morning.

[lukeolson]

What is the check in meshmode doing? https://github.com/inducer/meshmode/blob/48abe4bedf5e54e57fca392756423b6f45abff9a/meshmode/mesh/processing.py#L511

I think we can check the signed area/volume using the determinant (equivalent to curls):
https://github.com/lukeolson/mesh_orientation/blob/main/check.py#L24

Look it over... I dumped a few tests in there and it appears to work, but it's not exhaustive. When I generate isolator.msh it is all negative. When I check annulus.msh in meshmode, it's all positive.

[MTCam]

Look it over... I dumped a few tests in there and it appears to work, but it's not exhaustive. When I generate isolator.msh it is all negative.
✅

[mtcampbe@lassen709:data]$ python ~/check_mesh.py isolator.msh

Checking 3D mesh
found 0 positive elements
found 51730 negative elements
found 0 co-planar elements

Thanks! Now we can check meshes regardless of size before we try to run them in our allocations.

[anderson2981]

I don't get the error message from mesh mode for this mesh, but luke's code snippet indicates that it is inside out.

2d mesh

@inducer

[inducer]

I think there's an emerging consensus over in inducer/meshmode#314. I'd welcome thoughts in inducer/meshmode#315.

[nchristensen]

When I run Luke's script on eigthX/data/isolator.msh, it finds all elements are positive:

Checking 3D mesh
found 51840 positive  elements
found 0 negative  elements
found 0 co-planar elements

However, when I run eigthX/isolator.py it still finds negative elements. in this check.

Checking the values in area_elements[valid] gives the following:

positive valid area_elements: 100039
negative valid area_elements: 100064
0-value valid area_elements: 6429

[lukeolson]

@nchristensen does the new meshmode check throw an error?

[nchristensen]

Yes, this is with the current meshmode HEAD, d63d12a840101d20737b18dd8dd78b30143297de.

[nchristensen]

Ah, misinterpreted your question. I'll try the meshmode check.

[nchristensen]

Using that check on isolator.msh succeeds. https://gist.github.com/nchristensen/a3b47be19ac07f35fafe10bc1cd58d94

positive valid area_elements: 51633
negative valid area_elements: 0
0-value valid area_elements: 0

[nchristensen]

  File "/g/g15/njchris/Workspace/meshmode-upstream/meshmode/discretization/connection/face.py", line 345, in make_face_restriction
    bdry_mesh = Mesh(bdry_vertices, bdry_mesh_groups)

The problem appears to be in the boundary mesh.

[inducer]

But that's very strange. Mathematically, if the volume is positively oriented, so are the faces...?

[nchristensen]

    bdry_mesh = Mesh(bdry_vertices, bdry_mesh_groups)
  File "/home/njchris2/Workspace/emirge/meshmode/meshmode/mesh/__init__.py", line 1066, in __init__
    my_test_element_orientation_via_flipping_v2(self)
  File "/home/njchris2/Workspace/emirge/meshmode/meshmode/mesh/__init__.py", line 1061, in my_test_element_orientation_via_flipping_v2
    assert ((mesh_orient < 0) == (flippy > 0)).all()
AssertionError
flippy > 0: 53607
flippy < 0: 0
flippy == 0: 152925
mesh_orient < 0 99928
mesh_orient > 0 100174
mesh_orient == 0 6430

flippy and mesh_orient also seem to disagree for the boundary mesh.

[inducer]

Could upload one of the .msh files with this symptom and a script to reproduce the issue to a gist and link it here?

[MTCam]

I just wanted to add another data point. I said I would re-try these yesterday.

For me, with a fresh install using emirge/install.sh --branch=fusion_actx, I can no longer reproduce the problem with the meshes in drivers_y2-isolator. It would be shocking if we get different results on Lassen after:

git clone git@github.com:/illinois-ceesd/emirge
cd emirge && ./install.sh --branch=fusion_actx
# get nodes
# lalloc 1 -q pdebug 
source config/activate_env.sh
cd mirgecom
git clone git@github.com:/illinois-ceesd/drivers_y2-isolator
cd drivers_y2-isolator/smoke_test_3d
cd data && rm *.msh && ./make_mesh.sh && cd ../
# perform your favorite mesh checks -and/or- run mirgecom
python -m mpi4py ./isolator.py -i run_params.yaml

[nchristensen]

python -m mpi4py ./isolator.py -i run_params.yaml

The absence of -i run_params.yml seems to cause the negative orientation warning. With that included it runs for me.