Skip to content

Commit

Permalink
Fix "all" region for FCI
Browse files Browse the repository at this point in the history
  • Loading branch information
dschwoerer committed Oct 25, 2021
1 parent 00e5388 commit e1b1e1c
Show file tree
Hide file tree
Showing 2 changed files with 12 additions and 6 deletions.
14 changes: 10 additions & 4 deletions xbout/geometries.py
Original file line number Diff line number Diff line change
Expand Up @@ -184,7 +184,7 @@ def apply_geometry(ds, geometry_name, *, coordinates=None, grid=None):
# In BOUT++ v5, dz is either a Field2D or Field3D.
# We can use it as a 1D coordinate if it's a Field3D, _or_ if nz == 1
bout_v5 = updated_ds.metadata["BOUT_VERSION"] > 5.0 or (
updated_ds.metadata["BOUT_VERSION"] == 5.0 and updated_ds["dz"].ndim == 2
updated_ds.metadata["BOUT_VERSION"] == 5.0 and updated_ds["dz"].ndim >= 2
)
use_metric_3d = updated_ds.metadata.get("use_metric_3d", False)
can_use_1d_z_coord = (nz == 1) or use_metric_3d
Expand All @@ -197,14 +197,20 @@ def apply_geometry(ds, geometry_name, *, coordinates=None, grid=None):
raise ValueError(
f"Spacing is not constant. Cannot create z coordinate"
)
dz = updated_ds["dz"][0, 0]

dz = updated_ds["dz"].min()
else:
dz = updated_ds["dz"]

z0 = 2 * np.pi * updated_ds.metadata["ZMIN"]
z1 = z0 + nz * dz
if not np.isclose(
z1, 2.0 * np.pi * updated_ds.metadata["ZMAX"], rtol=1.0e-15, atol=0.0
if not np.all(
np.isclose(
z1,
2.0 * np.pi * updated_ds.metadata["ZMAX"],
rtol=1.0e-15,
atol=0.0,
)
):
warn(
f"Size of toroidal domain as calculated from nz*dz ({str(z1 - z0)}"
Expand Down
4 changes: 2 additions & 2 deletions xbout/region.py
Original file line number Diff line number Diff line change
Expand Up @@ -1239,8 +1239,8 @@ def _create_single_region(ds, periodic_y=True):
"all": Region(
name="all",
ds=ds,
xouter_ind=0,
xinner_ind=nx,
xouter_ind=nx,
xinner_ind=0,
ylower_ind=0,
yupper_ind=ny,
connection_lower_y=connection,
Expand Down

0 comments on commit e1b1e1c

Please sign in to comment.