`write_data_to_vtk` volume normalization correction by pshriwise · Pull Request #2397 · openmc-dev/openmc

pshriwise · 2023-02-24T05:56:14Z

I recently noticed that the volume normalization for structured mesh types (other than regular mesh) was producing strange results while working on #2359. This turned out to be because the method expects data in kji ordering for the elements (with i changing fastest), but when normalizing these results by the element volumes the MeshBase.volumes property returns them in ijk ordering. This PR fixes that by transposing the volumes array in this method.

I also noticed that vertices were being created manually, so I've updated the vertices arrays to leverage the MeshBase.vertices property instead.

… normalization

shimwell · 2023-02-24T10:06:36Z

Patrick you are an absolute star thanks so much for finding this. I think I've been suffering from this bug without realizing. Keen to see this PR go in

paulromano · 2023-02-24T19:44:59Z

@shimwell if you don't mind, it'd be great if you can take lead on reviewing this fix!

shimwell · 2023-02-27T22:33:09Z

I made a quick test that writes a vtk mesh and then reads it in and checks the values. It passes on the develop branch, but no longer passes on this branch. I would be keen to adapt this test so that it works on this PR and then include it in the PR to protect the write to vtk method for the future. What do you think @pshriwise

import vtk
import numpy as np
import openmc


def test_write_data_to_vtk_round_trip():
    cmesh = openmc.CylindricalMesh()
    cmesh.r_grid = (0.0, 1.0, 2.0)
    cmesh.z_grid = (0.0, 2.0, 4.0, 5.0)
    cmesh.phi_grid = (0.0, 3.0, 6.0)

    smesh = openmc.SphericalMesh()
    smesh.r_grid = (0.0, 1.0, 2.0)
    smesh.theta_grid = (0.0, 2.0, 4.0, 5.0)
    smesh.phi_grid = (0.0, 3.0, 6.0)

    rmesh = openmc.RegularMesh()
    rmesh.lower_left = (0.0, 0.0, 0.0)
    rmesh.upper_right = (1.0, 3.0, 5.0)
    rmesh.dimension = (2, 1, 6)

    for mesh in [smesh, cmesh, rmesh]:

        filename = "mesh.vtk"
        data = np.array([1.0] * 12)  # there are 12 voxels in each mesh
        mesh.write_data_to_vtk(
            filename=filename, datasets={"normalized": data}, volume_normalization=True
        )

        reader = vtk.vtkStructuredGridReader()
        reader.SetFileName(filename)
        reader.ReadAllFieldsOn()
        reader.Update()

        cell_data = reader.GetOutput().GetCellData()
        uniform_array = cell_data.GetArray("normalized")
        num_tuples = uniform_array.GetNumberOfTuples()
        vtk_values = [uniform_array.GetValue(i) for i in range(num_tuples)]

        # checks that the vtk cell values are equal to the data / mesh volumes
        assert np.allclose(vtk_values, data / mesh.volumes.flatten())

        mesh.write_data_to_vtk(
            filename=filename,
            datasets={"not_normalized": data},
            volume_normalization=False,
        )

        reader = vtk.vtkStructuredGridReader()
        reader.SetFileName(filename)
        reader.ReadAllFieldsOn()
        reader.Update()

        cell_data = reader.GetOutput().GetCellData()
        uniform_array = cell_data.GetArray("not_normalized")
        num_tuples = uniform_array.GetNumberOfTuples()
        vtk_values = [uniform_array.GetValue(i) for i in range(num_tuples)]

        # checks that the vtk cell values are equal to the data
        assert np.array_equal(vtk_values, data)

shimwell · 2023-02-28T10:24:42Z

I've been simulating a simple CylindricalMesh with a source at the bottom.

This fix is producing much more sensible results

Before fix (current develop branch)

After fix (this PR branch)

code if anyone wants to reproduce this

import openmc

surf1 = openmc.ZPlane(z0=10, boundary_type='vacuum')
surf2 = openmc.ZPlane(z0=0, boundary_type='vacuum')
surf3 = openmc.ZCylinder(r=10, boundary_type='vacuum')

cell1 = openmc.Cell(region=-surf1 & + surf2 & -surf3)

my_geometry = openmc.Geometry([cell1])

my_settings = openmc.Settings()
batches = 2
my_settings.batches = batches
my_settings.inactive = 0
my_settings.particles = 50000
my_settings.run_mode = 'fixed source'

source = openmc.Source()
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([14e6], [1])
source.space = openmc.stats.Point((0, 0, 0.1))

my_settings.source = source

mesh = openmc.CylindricalMesh().from_domain(cell1)
mesh_filter = openmc.MeshFilter(mesh)

mesh_tally_1 = openmc.Tally(name='flux_on_mesh')
mesh_tally_1.filters = [mesh_filter]
mesh_tally_1.scores = ['flux']  # where X is a wildcard
my_tallies = openmc.Tallies([mesh_tally_1])

model = openmc.Model(my_geometry, openmc.Materials([]), my_settings, my_tallies)
sp_filename = model.run()

statepoint = openmc.StatePoint(sp_filename)

my_mesh_tally = statepoint.get_tally(name='flux_on_mesh')

mesh.write_data_to_vtk(
    filename="norm_fix.vtk",
    datasets={"mean": my_mesh_tally.mean},
    volume_normalization=True
)

shimwell · 2023-02-28T11:37:48Z

I've also run some SphericalMesh simulations with this new branch

This fix is producing much more sensible results

Before fix (current develop branch)

After fix (this PR branch)

shimwell · 2023-02-28T11:43:19Z

In the case where the volume normalization is set to False for a SphericalMesh the results are still looking a bit strange

mesh.write_data_to_vtk(
    filename="sphere_no_norm_fix.vtk",
    datasets={"mean": my_mesh_tally.mean},
    volume_normalization=False
)

pshriwise · 2023-03-01T03:59:57Z

Whoops, did not mean to close this.

Thanks for all of the great examples and verification @shimwell!

I like the testing approach you provided. I'll include that in this branch.

I agree that the results in your last image look al little funny, but when a tally mean isn't adjusted by the element volumes I'd expect to see something that's hard to interpret for a spherical mesh. The larger elements on the outside are going to skew the visualization quite a bit I think. One way to check that the data is being ordered correctly would be to make the tally mesh larger than the geometry. This way we can verify that mesh elements that are entirely outside of the mesh are zero and checking the results by inspection will be easier. What do you think?

Co-authored-by: Jonathan Shimwell <jonathan.shimwell@firstlightfusion.com>

shimwell · 2023-03-01T16:42:50Z

I've been running this script to test the spherical mesh with and without normalisation

import openmc
import math
import numpy as np
surf1 = openmc.Sphere(r=10, boundary_type='vacuum')

cell1 = openmc.Cell(region=-surf1)

my_geometry = openmc.Geometry([cell1])

my_settings = openmc.Settings()
batches = 2
my_settings.batches = batches
my_settings.inactive = 0
my_settings.particles = 500000
my_settings.run_mode = 'fixed source'

source = openmc.Source()
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([14e6], [1])
source.space = openmc.stats.Point((0, 0, 0.1))

my_settings.source = source

mesh = openmc.SphericalMesh()
mesh.r_grid = (0,1,2,3,4,5,6,7,8,9,10)
mesh.theta_grid = np.linspace(0, math.pi, 20)
mesh.phi_grid = np.linspace(0, 2*math.pi, 20)
mesh_filter = openmc.MeshFilter(mesh)

mesh_tally_1 = openmc.Tally(name='flux_on_mesh')
mesh_tally_1.filters = [mesh_filter]
mesh_tally_1.scores = ['flux']  # where X is a wildcard
my_tallies = openmc.Tallies([mesh_tally_1])

model = openmc.Model(my_geometry, openmc.Materials([]), my_settings, my_tallies)
sp_filename = model.run()

statepoint = openmc.StatePoint(sp_filename)

my_mesh_tally = statepoint.get_tally(name='flux_on_mesh')

mesh.write_data_to_vtk(
    filename="sphere_no_norm_fixPR.vtk",
    datasets={"mean": my_mesh_tally.mean},
    volume_normalization=False
)
mesh.write_data_to_vtk(
    filename="sphere_norm_fixPR.vtk",
    datasets={"mean": my_mesh_tally.mean},
    volume_normalization=True
)

The sphere_no_norm_fixPR.vtk one is the one that shows odd results.
It looks like I'm getting a stripe of values across one axis so I thought this might be a systematic bug

pshriwise · 2023-03-06T19:57:43Z

Thanks for digging into that @shimwell! I was pretty perplexed by this too. In the end, I think you've found a bug in how we're performing tallies on spherical meshes. When I apply the same simulation with a cylindrical or rectilinear mesh I see more reasonable results for both the normalized and un-normalized cases.

Cylindrical

Rectilinear

So I'll propose that we move forward with this PR and I fix the spherical mesh problems in another PR. I have another one from @eepeterson to tackle as well. Is that alright with you @shimwell?

shimwell

Looks good to me, improvement in a few mesh types with and without volume normalization. Discovered another potential bug in spherical meshes that will be fixed in another PR

shimwell · 2023-03-06T20:41:01Z

Happy to merge this tomorrow if there are no objections

shimwell · 2023-03-07T20:19:26Z

Thanks @pshriwise

`write_data_to_vtk` volume normalization correction

pshriwise added 2 commits February 23, 2023 23:39

Using vertices property where possible. Correcting volume ordering in…

a270d4c

… normalization

Add test of volume normalization ordering

70b743e

pshriwise mentioned this pull request Feb 24, 2023

Global Weight Window Generation via CAPI #2359

Merged

2 of 3 tasks

paulromano requested a review from shimwell February 24, 2023 19:44

shimwell self-assigned this Feb 24, 2023

pshriwise closed this Mar 1, 2023

pshriwise reopened this Mar 1, 2023

Adding another test that uses VTK to gather the data directly.

728f48c

Co-authored-by: Jonathan Shimwell <jonathan.shimwell@firstlightfusion.com>

shimwell approved these changes Mar 6, 2023

View reviewed changes

shimwell added the Merging Soon PR will be merged in < 24 hrs if no further comments are made. label Mar 6, 2023

shimwell merged commit 9768e34 into openmc-dev:develop Mar 7, 2023

shimwell mentioned this pull request Mar 7, 2023

spherical mesh tally bug #2414

Closed

pshriwise mentioned this pull request Mar 8, 2023

Spherical mesh VTK coordinates bug #2416

Closed

pshriwise removed the Merging Soon PR will be merged in < 24 hrs if no further comments are made. label Apr 21, 2025

apingegno pushed a commit to apingegno/openmc that referenced this pull request May 7, 2026

Merge pull request openmc-dev#2397 from pshriwise/mesh-vtk-fixes

913b453

`write_data_to_vtk` volume normalization correction

Navigation Menu

Repository navigation

`write_data_to_vtk` volume normalization correction - #2397

`write_data_to_vtk` volume normalization correction#2397
shimwell merged 3 commits into
openmc-dev:developfrom
pshriwise:mesh-vtk-fixes

pshriwise commented Feb 24, 2023

Uh oh!

shimwell commented Feb 24, 2023

Uh oh!

paulromano commented Feb 24, 2023

Uh oh!

shimwell commented Feb 27, 2023

Uh oh!

shimwell commented Feb 28, 2023 •

edited

Loading

Uh oh!

shimwell commented Feb 28, 2023

Uh oh!

shimwell commented Feb 28, 2023

Uh oh!

pshriwise commented Mar 1, 2023

Uh oh!

shimwell commented Mar 1, 2023

Uh oh!

pshriwise commented Mar 6, 2023

Uh oh!

shimwell left a comment

Uh oh!

shimwell commented Mar 6, 2023

Uh oh!

shimwell commented Mar 7, 2023

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Conversation

pshriwise commented Feb 24, 2023

Uh oh!

shimwell commented Feb 24, 2023

Uh oh!

paulromano commented Feb 24, 2023

Uh oh!

shimwell commented Feb 27, 2023

Uh oh!

shimwell commented Feb 28, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

shimwell commented Feb 28, 2023

Uh oh!

shimwell commented Feb 28, 2023

Uh oh!

pshriwise commented Mar 1, 2023

Uh oh!

shimwell commented Mar 1, 2023

Uh oh!

pshriwise commented Mar 6, 2023

Cylindrical

Rectilinear

Uh oh!

shimwell left a comment

Choose a reason for hiding this comment

Uh oh!

shimwell commented Mar 6, 2023

Uh oh!

shimwell commented Mar 7, 2023

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

shimwell commented Feb 28, 2023 •

edited

Loading