diff --git a/.github/workflows/wheels.yaml b/.github/workflows/wheels.yaml index 494eed45a2..c9ee27601f 100644 --- a/.github/workflows/wheels.yaml +++ b/.github/workflows/wheels.yaml @@ -15,10 +15,13 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: + # pinning OS versions to oldest versions available + # because forward compatibility is much easier to + # guarantee than backward compatibility os: [ - macos-latest, - windows-latest, - ubuntu-18.04, # has to be the oldest possible for manylinux + ubuntu-18.04, + windows-2019, + macos-10.15, ] fail-fast: false @@ -27,6 +30,8 @@ jobs: - name: Build wheels for CPython uses: pypa/cibuildwheel@v2.3.0 + with: + output-dir: dist env: CIBW_BUILD: "cp36-* cp37-* cp38-* cp39-* cp310-*" CIBW_SKIP: "*-musllinux_*" # numpy doesn't have wheels for musllinux so we can't build some quickly and without bloating @@ -40,4 +45,39 @@ jobs: - uses: actions/upload-artifact@v2 with: name: wheels - path: ./wheelhouse/*.whl + path: ./dist/*.whl + + build_sdist: + name: Build source distribution + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v2 + + - name: Build sdist + run: pipx run build --sdist + + - uses: actions/upload-artifact@v2 + with: + name: tarball + path: dist/*.tar.gz + + upload_pypi: + needs: [build_wheels, build_sdist] + runs-on: ubuntu-latest + # upload to PyPI on every tag starting with 'yt-' + if: github.event_name == 'push' && startsWith(github.event.ref, 'refs/tags/yt-') + steps: + - uses: actions/download-artifact@v2 + with: + name: tarball + path: dist + + - uses: actions/download-artifact@v2 + with: + name: wheels + path: dist + + - uses: pypa/gh-action-pypi-publish@v1.4.2 + with: + user: __token__ + password: ${{ secrets.pypi_token }} diff --git a/.gitignore b/.gitignore index bdd66181ea..b251fe0d79 100644 --- a/.gitignore +++ b/.gitignore @@ -7,6 +7,7 @@ png.cfg rockstar.cfg yt_updater.log yt/frontends/artio/_artio_caller.c +yt/frontends/gamer/cfields.c yt/analysis_modules/halo_finding/rockstar/rockstar_groupies.c yt/analysis_modules/halo_finding/rockstar/rockstar_interface.c yt/analysis_modules/ppv_cube/ppv_utils.c diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 841ce61280..f75c2ac06d 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -27,6 +27,7 @@ repos: - id: end-of-file-fixer - id: no-commit-to-branch args: [--branch, main] + - id: check-yaml - repo: https://github.com/asottile/pyupgrade rev: v2.19.4 hooks: diff --git a/doc/source/developing/creating_frontend.rst b/doc/source/developing/creating_frontend.rst index 7790a95d49..96637b5333 100644 --- a/doc/source/developing/creating_frontend.rst +++ b/doc/source/developing/creating_frontend.rst @@ -32,8 +32,8 @@ If you are interested in adding a new code, be sure to drop us a line on `yt-dev `_! -Boostraping a new frontend --------------------------- +Bootstrapping a new frontend +---------------------------- To get started diff --git a/doc/source/developing/deprecating_features.rst b/doc/source/developing/deprecating_features.rst index 025f93f892..34fbbd03e9 100644 --- a/doc/source/developing/deprecating_features.rst +++ b/doc/source/developing/deprecating_features.rst @@ -21,7 +21,7 @@ Here's an example call. issue_deprecation_warning( "`old_function` is deprecated, use `replacement_function` instead." since="4.0.0", - removal="4.1.0" + removal="4.1.0", ) ... diff --git a/doc/source/examining/Loading_Generic_Array_Data.ipynb b/doc/source/examining/Loading_Generic_Array_Data.ipynb index 72ffade3e3..400be1fb99 100644 --- a/doc/source/examining/Loading_Generic_Array_Data.ipynb +++ b/doc/source/examining/Loading_Generic_Array_Data.ipynb @@ -103,7 +103,7 @@ "* `mass_unit` : The unit that corresponds to `code_mass`, can be a string, tuple, or floating-point number\n", "* `time_unit` : The unit that corresponds to `code_time`, can be a string, tuple, or floating-point number\n", "* `velocity_unit` : The unit that corresponds to `code_velocity`\n", - "* `magnetic_unit` : The unit that corresponds to `code_magnetic`, i.e. the internal units used to represent magnetic field strengths.\n", + "* `magnetic_unit` : The unit that corresponds to `code_magnetic`, i.e. the internal units used to represent magnetic field strengths. NOTE: if you want magnetic field units to be in the SI unit system, you must specify it here, e.g. `magnetic_unit=(1.0, \"T\")`\n", "* `periodicity` : A tuple of booleans that determines whether the data will be treated as periodic along each axis\n", "\n", "This example creates a yt-native dataset `ds` that will treat your array as a\n", diff --git a/setup.cfg b/setup.cfg index 88b70cfa4a..f7644a2e42 100644 --- a/setup.cfg +++ b/setup.cfg @@ -65,6 +65,7 @@ nose.plugins.0.10 = doc = alabaster bottle + jupyter-client<7.0 nbconvert==5.6.1 pyregion pyx>=0.15 diff --git a/tests/tests.yaml b/tests/tests.yaml index 07aa928b6f..4e9d735ed6 100644 --- a/tests/tests.yaml +++ b/tests/tests.yaml @@ -28,8 +28,10 @@ answer_tests: - yt/frontends/athena/tests/test_outputs.py:test_blast - yt/frontends/athena/tests/test_outputs.py:test_stripping - local_athena_pp_004: - - yt/frontends/athena_pp/tests/test_outputs.py:test_disk + local_athena_pp_006: + # disabling disk test for now until we have better support + # for this dataset + #- yt/frontends/athena_pp/tests/test_outputs.py:test_disk - yt/frontends/athena_pp/tests/test_outputs.py:test_AM06 local_chombo_005: diff --git a/yt/data_objects/static_output.py b/yt/data_objects/static_output.py index 0a38e5ed90..8c8b862625 100644 --- a/yt/data_objects/static_output.py +++ b/yt/data_objects/static_output.py @@ -1119,27 +1119,47 @@ def relative_refinement(self, l0, l1): return self.refine_by ** (l1 - l0) def _assign_unit_system(self, unit_system): - if unit_system == "cgs": - current_mks_unit = None + # we need to determine if the requested unit system + # is mks-like: i.e., it has a current with the same + # dimensions as amperes. + mks_system = False + if getattr(self, "magnetic_unit", None): + mag_dims = self.magnetic_unit.units.dimensions.free_symbols else: - current_mks_unit = "A" - magnetic_unit = getattr(self, "magnetic_unit", None) - if magnetic_unit is not None: - if unit_system == "mks": - if current_mks not in self.magnetic_unit.units.dimensions.free_symbols: - self.magnetic_unit = self.magnetic_unit.to("gauss").to("T") - self.unit_registry.modify("code_magnetic", self.magnetic_unit.value) - else: - # if the magnetic unit is in T, we need to create the code unit - # system as an MKS-like system - if current_mks in self.magnetic_unit.units.dimensions.free_symbols: - self.magnetic_unit = self.magnetic_unit.to("T").to("gauss") + mag_dims = None + if unit_system != "code": + # if the unit system is known, we can check if it + # has a "current_mks" unit + us = unit_system_registry[str(unit_system).lower()] + mks_system = us.base_units[current_mks] is not None + elif mag_dims and current_mks in mag_dims: + # if we're using the not-yet defined code unit system, + # then we check if the magnetic field unit has a SI + # current dimension in it + mks_system = True + # Now we get to the tricky part. If we have an MKS-like system but + # we asked for a conversion to something CGS-like, or vice-versa, + # we have to convert the magnetic field + if mag_dims is not None: + if mks_system and current_mks not in mag_dims: + self.magnetic_unit = self.quan( + self.magnetic_unit.to_value("gauss") * 1.0e-4, "T" + ) + # The following modification ensures that we get the conversion to + # mks correct + self.unit_registry.modify( + "code_magnetic", self.magnetic_unit.value * 1.0e3 * 0.1 ** -0.5 + ) + elif not mks_system and current_mks in mag_dims: + self.magnetic_unit = self.quan( + self.magnetic_unit.to_value("T") * 1.0e4, "gauss" + ) # The following modification ensures that we get the conversion to # cgs correct self.unit_registry.modify( - "code_magnetic", self.magnetic_unit.value * 0.1 ** 0.5 + "code_magnetic", self.magnetic_unit.value * 1.0e-4 ) - + current_mks_unit = "A" if mks_system else None us = create_code_unit_system( self.unit_registry, current_mks_unit=current_mks_unit ) @@ -1165,25 +1185,34 @@ def _create_unit_registry(self, unit_system): # yt assumes a CGS unit system by default (for back compat reasons). # Since unyt is MKS by default we specify the MKS values of the base # units in the CGS system. So, for length, 1 cm = .01 m. And so on. + # Note that the values associated with the code units here will be + # modified once we actually determine what the code units are from + # the dataset + # NOTE that magnetic fields are not done here yet, see set_code_units self.unit_registry = UnitRegistry(unit_system=unit_system) + # 1 cm = 0.01 m self.unit_registry.add("code_length", 0.01, dimensions.length) + # 1 g = 0.001 kg self.unit_registry.add("code_mass", 0.001, dimensions.mass) + # 1 g/cm**3 = 1000 kg/m**3 self.unit_registry.add("code_density", 1000.0, dimensions.density) + # 1 erg/g = 1.0e-4 J/kg self.unit_registry.add( - "code_specific_energy", 1.0, dimensions.energy / dimensions.mass + "code_specific_energy", 1.0e-4, dimensions.energy / dimensions.mass ) + # 1 s = 1 s self.unit_registry.add("code_time", 1.0, dimensions.time) - if unit_system == "mks": - self.unit_registry.add("code_magnetic", 1.0, dimensions.magnetic_field) - else: - self.unit_registry.add( - "code_magnetic", 0.1 ** 0.5, dimensions.magnetic_field_cgs - ) + # 1 K = 1 K self.unit_registry.add("code_temperature", 1.0, dimensions.temperature) + # 1 dyn/cm**2 = 0.1 N/m**2 self.unit_registry.add("code_pressure", 0.1, dimensions.pressure) + # 1 cm/s = 0.01 m/s self.unit_registry.add("code_velocity", 0.01, dimensions.velocity) + # metallicity self.unit_registry.add("code_metallicity", 1.0, dimensions.dimensionless) + # dimensionless hubble parameter self.unit_registry.add("h", 1.0, dimensions.dimensionless, r"h") + # cosmological scale factor self.unit_registry.add("a", 1.0, dimensions.dimensionless) def set_units(self): @@ -1285,6 +1314,24 @@ def set_code_units(self): self.unit_registry.modify("code_pressure", pressure_unit) self.unit_registry.modify("code_density", density_unit) self.unit_registry.modify("code_specific_energy", specific_energy_unit) + # Defining code units for magnetic fields are tricky because + # they have different dimensions in different unit systems, so we have + # to handle them carefully + if hasattr(self, "magnetic_unit"): + if self.magnetic_unit.units.dimensions == dimensions.magnetic_field_cgs: + # We have to cast this explicitly to MKS base units, otherwise + # unyt will convert it automatically to Tesla + value = self.magnetic_unit.to_value("sqrt(kg)/(sqrt(m)*s)") + dims = dimensions.magnetic_field_cgs + else: + value = self.magnetic_unit.to_value("T") + dims = dimensions.magnetic_field_mks + else: + # Fallback to gauss if no magnetic unit is specified + # 1 gauss = 1 sqrt(g)/(sqrt(cm)*s) = 0.1**0.5 sqrt(kg)/(sqrt(m)*s) + value = 0.1 ** 0.5 + dims = dimensions.magnetic_field_cgs + self.unit_registry.add("code_magnetic", value, dims) # domain_width does not yet exist if self.domain_left_edge is not None and self.domain_right_edge is not None: DW = self.arr(self.domain_right_edge - self.domain_left_edge, "code_length") diff --git a/yt/fields/derived_field.py b/yt/fields/derived_field.py index 82ed27258d..5e0aad640c 100644 --- a/yt/fields/derived_field.py +++ b/yt/fields/derived_field.py @@ -324,7 +324,7 @@ def get_label(self, projected=False): units = Unit(self.units) # Add unit label if not units.is_dimensionless: - data_label += r"\ \ (%s)" % (units.latex_representation()) + data_label += r"\ \ \left(%s\right)" % (units.latex_representation()) data_label += r"$" return data_label diff --git a/yt/fields/field_info_container.py b/yt/fields/field_info_container.py index 8677544550..6caa86c59d 100644 --- a/yt/fields/field_info_container.py +++ b/yt/fields/field_info_container.py @@ -296,9 +296,8 @@ def _sanitize_sampling_type(sampling_type, particle_type=None): acceptable_samplings = ("cell", "particle", "local") if sampling_type not in acceptable_samplings: raise ValueError( - "Invalid sampling type %s. Valid sampling types are %s", - sampling_type, - ", ".join(acceptable_samplings), + f"Received invalid sampling type {sampling_type!r}. " + f"Expected any of {acceptable_samplings}" ) if particle_type: diff --git a/yt/fields/fluid_fields.py b/yt/fields/fluid_fields.py index acd13344dc..2e365d6ef2 100644 --- a/yt/fields/fluid_fields.py +++ b/yt/fields/fluid_fields.py @@ -1,6 +1,7 @@ import numpy as np -from yt.units.unit_object import Unit +from yt.units.dimensions import current_mks # type: ignore +from yt.units.unit_object import Unit # type: ignore from yt.utilities.chemical_formulas import compute_mu from yt.utilities.lib.misc_utilities import obtain_relative_velocity_vector @@ -31,7 +32,7 @@ def setup_fluid_fields(registry, ftype="gas", slice_info=None): unit_system = registry.ds.unit_system - if unit_system.name == "cgs": + if unit_system.base_units[current_mks] is None: mag_units = "magnetic_field_cgs" else: mag_units = "magnetic_field_mks" diff --git a/yt/fields/tests/test_ambiguous_fields.py b/yt/fields/tests/test_ambiguous_fields.py index a5dbf25dbf..2d4955505a 100644 --- a/yt/fields/tests/test_ambiguous_fields.py +++ b/yt/fields/tests/test_ambiguous_fields.py @@ -1,3 +1,4 @@ +import warnings from collections import defaultdict import pytest @@ -32,13 +33,12 @@ def test_ambiguous_fails(): # Test no warnings are issued for single fname access that aren't ambiguous for fname in unambiguous_fnames: - with pytest.warns(None) as record: + with warnings.catch_warnings(): + warnings.simplefilter("error") ds.r[fname] - assert len(record) == 0 # Test no warning are issued for tuple access for ftype, fname in ds.field_list: - with pytest.warns(None) as record: + with warnings.catch_warnings(): + warnings.simplefilter("error") ds.r[(ftype, fname)] - - assert len(record) == 0 diff --git a/yt/fields/tests/test_magnetic_fields.py b/yt/fields/tests/test_magnetic_fields.py index 4258eb4b73..709926e6a1 100644 --- a/yt/fields/tests/test_magnetic_fields.py +++ b/yt/fields/tests/test_magnetic_fields.py @@ -23,19 +23,27 @@ def test_magnetic_fields(): for field in data1: data2[field] = (data1[field][0] * 1.0e4, "gauss") - ds1 = load_uniform_grid(data1, ddims, unit_system="cgs") + ds0 = load_uniform_grid(data1, ddims, unit_system="cgs") + ds1 = load_uniform_grid(data1, ddims, magnetic_unit=(1.0, "T"), unit_system="cgs") ds2 = load_uniform_grid(data2, ddims, unit_system="mks") # For this test dataset, code units are cgs units ds3 = load_uniform_grid(data2, ddims, unit_system="code") + # For this test dataset, code units are SI units + ds4 = load_uniform_grid(data1, ddims, magnetic_unit=(1.0, "T"), unit_system="code") + ds0.index ds1.index ds2.index ds3.index + ds4.index + dd0 = ds0.all_data() dd1 = ds1.all_data() dd2 = ds2.all_data() dd3 = ds3.all_data() + dd4 = ds4.all_data() + assert ds0.fields.gas.magnetic_field_strength.units == "G" assert ds1.fields.gas.magnetic_field_strength.units == "G" assert ds1.fields.gas.magnetic_field_poloidal.units == "G" assert ds1.fields.gas.magnetic_field_toroidal.units == "G" @@ -45,6 +53,16 @@ def test_magnetic_fields(): assert ds3.fields.gas.magnetic_field_strength.units == "code_magnetic" assert ds3.fields.gas.magnetic_field_poloidal.units == "code_magnetic" assert ds3.fields.gas.magnetic_field_toroidal.units == "code_magnetic" + assert ds4.fields.gas.magnetic_field_strength.units == "code_magnetic" + assert ds4.fields.gas.magnetic_field_poloidal.units == "code_magnetic" + assert ds4.fields.gas.magnetic_field_toroidal.units == "code_magnetic" + + emag0 = ( + dd0[("gas", "magnetic_field_x")] ** 2 + + dd0[("gas", "magnetic_field_y")] ** 2 + + dd0[("gas", "magnetic_field_z")] ** 2 + ) / (8.0 * np.pi) + emag0.convert_to_units("dyne/cm**2") emag1 = ( dd1[("gas", "magnetic_field_x")] ** 2 @@ -64,16 +82,32 @@ def test_magnetic_fields(): dd3[("gas", "magnetic_field_x")] ** 2 + dd3[("gas", "magnetic_field_y")] ** 2 + dd3[("gas", "magnetic_field_z")] ** 2 - ) / (2.0 * mu_0) + ) / (8.0 * np.pi) emag3.convert_to_units("code_pressure") + emag4 = ( + dd4[("gas", "magnetic_field_x")] ** 2 + + dd4[("gas", "magnetic_field_y")] ** 2 + + dd4[("gas", "magnetic_field_z")] ** 2 + ) / (2.0 * mu_0) + emag4.convert_to_units("code_pressure") + + # note that "magnetic_energy_density" and "magnetic_pressure" are aliased + + assert_almost_equal(emag0, dd0[("gas", "magnetic_energy_density")]) assert_almost_equal(emag1, dd1[("gas", "magnetic_energy_density")]) assert_almost_equal(emag2, dd2[("gas", "magnetic_energy_density")]) assert_almost_equal(emag3, dd3[("gas", "magnetic_energy_density")]) + assert_almost_equal(emag4, dd4[("gas", "magnetic_energy_density")]) + assert str(emag0.units) == str(dd0[("gas", "magnetic_energy_density")].units) assert str(emag1.units) == str(dd1[("gas", "magnetic_energy_density")].units) assert str(emag2.units) == str(dd2[("gas", "magnetic_energy_density")].units) assert str(emag3.units) == str(dd3[("gas", "magnetic_energy_density")].units) + assert str(emag4.units) == str(dd4[("gas", "magnetic_energy_density")].units) + assert_almost_equal(emag1.in_cgs(), emag0.in_cgs()) + assert_almost_equal(emag2.in_cgs(), emag0.in_cgs()) assert_almost_equal(emag1.in_cgs(), emag2.in_cgs()) assert_almost_equal(emag1.in_cgs(), emag3.in_cgs()) + assert_almost_equal(emag1.in_cgs(), emag4.in_cgs()) diff --git a/yt/frontends/_skeleton/data_structures.py b/yt/frontends/_skeleton/data_structures.py index 3e7923910b..fc0e4d60e5 100644 --- a/yt/frontends/_skeleton/data_structures.py +++ b/yt/frontends/_skeleton/data_structures.py @@ -114,6 +114,10 @@ def _set_code_unit_attributes(self): # These can also be set: # self.velocity_unit = self.quan(1.0, "cm/s") # self.magnetic_unit = self.quan(1.0, "gauss") + # + # If your frontend uses SI EM units, set magnetic units like this + # instead: + # self.magnetic_unit = self.quan(1.0, "T") # this minimalistic implementation fills the requirements for # this frontend to run, change it to make it run _correctly_ ! diff --git a/yt/frontends/athena_pp/tests/test_outputs.py b/yt/frontends/athena_pp/tests/test_outputs.py index 16340299dc..1ef837bf81 100644 --- a/yt/frontends/athena_pp/tests/test_outputs.py +++ b/yt/frontends/athena_pp/tests/test_outputs.py @@ -1,20 +1,16 @@ -import numpy as np - from yt.frontends.athena_pp.api import AthenaPPDataset from yt.loaders import load -from yt.testing import ( - assert_allclose, - assert_equal, - requires_file, - units_override_check, -) +from yt.testing import assert_equal, requires_file, units_override_check +from yt.units import dimensions from yt.utilities.answer_testing.framework import ( - GenericArrayTest, data_dir_load, requires_ds, small_patch_amr, ) +# Deactivating this problematic test until the dataset type can be +# handled properly, see https://github.com/yt-project/yt/issues/3619 +""" _fields_disk = ("density", "velocity_r") disk = "KeplerianDisk/disk.out1.00000.athdf" @@ -35,7 +31,7 @@ def field_func(name): return dd[field] yield GenericArrayTest(ds, field_func, args=[field]) - +""" _fields_AM06 = ("temperature", "density", "velocity_magnitude", "magnetic_field_x") @@ -74,3 +70,10 @@ def test_units_override(): @requires_file(AM06) def test_AthenaPPDataset(): assert isinstance(data_dir_load(AM06), AthenaPPDataset) + + +@requires_file(AM06) +def test_magnetic_units(): + ds = load(AM06, unit_system="code") + assert ds.magnetic_unit.units.dimensions == dimensions.magnetic_field_cgs + assert (ds.magnetic_unit ** 2).units.dimensions == dimensions.pressure diff --git a/yt/frontends/boxlib/tests/test_outputs.py b/yt/frontends/boxlib/tests/test_outputs.py index 9ae07a7148..e426793fd2 100644 --- a/yt/frontends/boxlib/tests/test_outputs.py +++ b/yt/frontends/boxlib/tests/test_outputs.py @@ -8,7 +8,14 @@ OrionDataset, WarpXDataset, ) -from yt.testing import assert_equal, requires_file, units_override_check +from yt.loaders import load +from yt.testing import ( + assert_allclose, + assert_equal, + disable_dataset_cache, + requires_file, + units_override_check, +) from yt.utilities.answer_testing.framework import ( GridValuesTest, data_dir_load, @@ -315,11 +322,28 @@ def test_CastroDataset_2(): assert isinstance(data_dir_load("castro_sod_x_plt00036"), CastroDataset) -@requires_file(LyA) +@requires_file(plasma) def test_WarpXDataset(): assert isinstance(data_dir_load(plasma), WarpXDataset) +@disable_dataset_cache +@requires_file(plasma) +def test_magnetic_units(): + ds1 = load(plasma) + assert_allclose(ds1.magnetic_unit.value, 1.0) + assert str(ds1.magnetic_unit.units) == "T" + mag_unit1 = ds1.magnetic_unit.to("code_magnetic") + assert_allclose(mag_unit1.value, 1.0) + assert str(mag_unit1.units) == "code_magnetic" + ds2 = load(plasma, unit_system="cgs") + assert_allclose(ds2.magnetic_unit.value, 1.0e4) + assert str(ds2.magnetic_unit.units) == "G" + mag_unit2 = ds2.magnetic_unit.to("code_magnetic") + assert_allclose(mag_unit2.value, 1.0) + assert str(mag_unit2.units) == "code_magnetic" + + @requires_ds(laser) def test_WarpXDataset_2(): assert isinstance(data_dir_load(laser), WarpXDataset) diff --git a/yt/frontends/flash/tests/test_outputs.py b/yt/frontends/flash/tests/test_outputs.py index 176b85a457..778f31b05c 100644 --- a/yt/frontends/flash/tests/test_outputs.py +++ b/yt/frontends/flash/tests/test_outputs.py @@ -3,9 +3,12 @@ import numpy as np from yt.frontends.flash.api import FLASHDataset, FLASHParticleDataset +from yt.loaders import load from yt.testing import ( ParticleSelectionComparison, + assert_allclose, assert_equal, + disable_dataset_cache, requires_file, units_override_check, ) @@ -54,6 +57,23 @@ def test_units_override(): units_override_check(sloshing) +@disable_dataset_cache +@requires_file(sloshing) +def test_magnetic_units(): + ds1 = load(sloshing) + assert_allclose(ds1.magnetic_unit.value, np.sqrt(4.0 * np.pi)) + assert str(ds1.magnetic_unit.units) == "G" + mag_unit1 = ds1.magnetic_unit.to("code_magnetic") + assert_allclose(mag_unit1.value, 1.0) + assert str(mag_unit1.units) == "code_magnetic" + ds2 = load(sloshing, unit_system="mks") + assert_allclose(ds2.magnetic_unit.value, np.sqrt(4.0 * np.pi) * 1.0e-4) + assert str(ds2.magnetic_unit.units) == "T" + mag_unit2 = ds2.magnetic_unit.to("code_magnetic") + assert_allclose(mag_unit2.value, 1.0) + assert str(mag_unit2.units) == "code_magnetic" + + @requires_file(sloshing) def test_mu(): ds = data_dir_load(sloshing) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index ebfcfefd95..b6ef908adf 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -1,6 +1,7 @@ import os import weakref from collections import defaultdict +from itertools import product from pathlib import Path import numpy as np @@ -41,16 +42,20 @@ class RAMSESFileSanitizer: group_name = None # str | None: name of the first group folder (if any) def __init__(self, filename): - # Resolve so that it works with symlinks - filename = Path(filename).resolve() + + # Make the resolve optional, so that it works with symlinks + paths_to_try = (Path(filename), Path(filename).resolve()) self.original_filename = filename self.output_dir = None self.info_fname = None - for check_fun in (self.test_with_standard_file, self.test_with_folder_name): - ok, output_dir, info_fname = check_fun(filename) + check_functions = (self.test_with_standard_file, self.test_with_folder_name) + + # Loop on both the functions and the tested paths + for path, check_fun in product(paths_to_try, check_functions): + ok, output_dir, info_fname = check_fun(path) if ok: break diff --git a/yt/frontends/stream/data_structures.py b/yt/frontends/stream/data_structures.py index da3fa1e6ae..fa4e772bee 100644 --- a/yt/frontends/stream/data_structures.py +++ b/yt/frontends/stream/data_structures.py @@ -334,7 +334,12 @@ def _set_code_unit_attributes(self): cgs_units = ("cm", "g", "s", "cm/s", "gauss") for unit, attr, cgs_unit in zip(base_units, attrs, cgs_units): if isinstance(unit, str): - uq = self.quan(1.0, unit) + if unit == "code_magnetic": + # If no magnetic unit was explicitly specified + # we skip it now and take care of it at the bottom + continue + else: + uq = self.quan(1.0, unit) elif isinstance(unit, numeric_type): uq = self.quan(unit, cgs_unit) elif isinstance(unit, YTQuantity): @@ -344,6 +349,10 @@ def _set_code_unit_attributes(self): else: raise RuntimeError(f"{attr} ({unit}) is invalid.") setattr(self, attr, uq) + if not hasattr(self, "magnetic_unit"): + self.magnetic_unit = np.sqrt( + 4 * np.pi * self.mass_unit / (self.time_unit ** 2 * self.length_unit) + ) @classmethod def _is_valid(cls, filename, *args, **kwargs): diff --git a/yt/frontends/ytdata/data_structures.py b/yt/frontends/ytdata/data_structures.py index cfdda49ccd..0a520d9e64 100644 --- a/yt/frontends/ytdata/data_structures.py +++ b/yt/frontends/ytdata/data_structures.py @@ -78,7 +78,9 @@ def _parse_parameter_file(self): if cu not in self.unit_registry: self.unit_registry.add(cu, 1.0, getattr(dimensions, dim)) if "code_magnetic" not in self.unit_registry: - self.unit_registry.add("code_magnetic", 1.0, dimensions.magnetic_field) + self.unit_registry.add( + "code_magnetic", 0.1 ** 0.5, dimensions.magnetic_field_cgs + ) # if saved, set unit system if "unit_system_name" in self.parameters: diff --git a/yt/units/tests/test_magnetic_code_units.py b/yt/units/tests/test_magnetic_code_units.py index 7915c5069d..f7bdb9fe09 100644 --- a/yt/units/tests/test_magnetic_code_units.py +++ b/yt/units/tests/test_magnetic_code_units.py @@ -22,7 +22,6 @@ def test_magnetic_code_units(): assert str(mucu.units) == "code_magnetic" ds2 = load_uniform_grid(data, ddims, magnetic_unit=(1.0, "T"), unit_system="cgs") - assert_allclose(ds2.magnetic_unit.value, 10000.0) assert str(ds2.magnetic_unit.units) == "G" @@ -49,3 +48,25 @@ def test_magnetic_code_units(): mucu = ds4.magnetic_unit.to("code_magnetic") assert_allclose(mucu.value, 1.0) assert str(mucu.units) == "code_magnetic" + + ds5 = load_uniform_grid( + data, ddims, magnetic_unit=(1.0, "uG"), unit_system="mks" + ) + + assert_allclose(ds5.magnetic_unit.value, 1.0e-10) + assert str(ds5.magnetic_unit.units) == "T" + + mucu = ds5.magnetic_unit.to("code_magnetic") + assert_allclose(mucu.value, 1.0) + assert str(mucu.units) == "code_magnetic" + + ds6 = load_uniform_grid( + data, ddims, magnetic_unit=(1.0, "nT"), unit_system="cgs" + ) + + assert_allclose(ds6.magnetic_unit.value, 1.0e-5) + assert str(ds6.magnetic_unit.units) == "G" + + mucu = ds6.magnetic_unit.to("code_magnetic") + assert_allclose(mucu.value, 1.0) + assert str(mucu.units) == "code_magnetic" diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index a5cabea7c4..a50f236484 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -525,7 +525,8 @@ def pixelize_cylinder(np.float64_t[:,:] buff, cdef np.float64_t x, y, dx, dy, r0, theta0 cdef np.float64_t rmax, x0, y0, x1, y1 - cdef np.float64_t r_i, theta_i, dr_i, dtheta_i, dthetamin + cdef np.float64_t r_i, theta_i, dr_i, dtheta_i + cdef np.float64_t r_inc, theta_inc cdef np.float64_t costheta, sintheta cdef int i, pi, pj @@ -558,9 +559,9 @@ def pixelize_cylinder(np.float64_t[:,:] buff, rbounds[0] = 0.0 if y0 < 0 and y1 > 0: rbounds[0] = 0.0 - dthetamin = dx / rmax - for i in range(radius.shape[0]): + r_inc = 0.5 * fmin(dx, dy) + for i in range(radius.shape[0]): r0 = radius[i] theta0 = theta[i] dr_i = dradius[i] @@ -569,15 +570,15 @@ def pixelize_cylinder(np.float64_t[:,:] buff, if r0 + dr_i < rbounds[0] or r0 - dr_i > rbounds[1]: continue theta_i = theta0 - dtheta_i - # Buffer of 0.5 here - dthetamin = 0.5*dx/(r0 + dr_i) + theta_inc = r_inc / (r0 + dr_i) + while theta_i < theta0 + dtheta_i: r_i = r0 - dr_i costheta = math.cos(theta_i) sintheta = math.sin(theta_i) while r_i < r0 + dr_i: if rmax <= r_i: - r_i += 0.5*dx + r_i += r_inc continue y = r_i * costheta x = r_i * sintheta @@ -586,8 +587,8 @@ def pixelize_cylinder(np.float64_t[:,:] buff, if pi >= 0 and pi < buff.shape[0] and \ pj >= 0 and pj < buff.shape[1]: buff[pi, pj] = field[i] - r_i += 0.5*dx - theta_i += dthetamin + r_i += r_inc + theta_i += theta_inc cdef int aitoff_Lambda_btheta_to_xy(np.float64_t Lambda, np.float64_t btheta, np.float64_t *x, np.float64_t *y) except -1: diff --git a/yt/visualization/fits_image.py b/yt/visualization/fits_image.py index bc65793586..a3ceadb952 100644 --- a/yt/visualization/fits_image.py +++ b/yt/visualization/fits_image.py @@ -399,21 +399,22 @@ def _set_units(self, ds, base_units): else: uq = None - if uq is not None and hasattr(self, "hubble_constant"): - # Don't store cosmology units + if uq is not None: atoms = {str(a) for a in uq.units.expr.atoms()} - if str(uq.units).startswith("cm") or "h" in atoms or "a" in atoms: + if hasattr(self, "hubble_constant"): + # Don't store cosmology units + if str(uq.units).startswith("cm") or "h" in atoms or "a" in atoms: + uq.convert_to_cgs() + if any(a.startswith("code") for a in atoms): + # Don't store code units + mylog.warning( + "Cannot use code units of '%s' " + "when creating a FITSImageData instance! " + "Converting to a cgs equivalent.", + uq.units, + ) uq.convert_to_cgs() - if uq is not None and uq.units.is_code_unit: - mylog.warning( - "Cannot use code units of '%s' " - "when creating a FITSImageData instance! " - "Converting to a cgs equivalent.", - uq.units, - ) - uq.convert_to_cgs() - if attr == "length_unit" and uq.value != 1.0: mylog.warning("Converting length units from %s to %s.", uq, uq.units) uq = YTQuantity(1.0, uq.units) diff --git a/yt/visualization/plot_modifications.py b/yt/visualization/plot_modifications.py index 23aaa5091a..4e0d6bd64c 100644 --- a/yt/visualization/plot_modifications.py +++ b/yt/visualization/plot_modifications.py @@ -3051,7 +3051,8 @@ def __call__(self, plot): vectors = np.concatenate((pixX[..., np.newaxis], pixY[..., np.newaxis]), axis=2) if self.texture is None: - self.texture = np.random.rand(nx, ny).astype(np.double) + prng = np.random.RandomState(0x4D3D3D3) + self.texture = prng.random_sample((nx, ny)) elif self.texture.shape != (nx, ny): raise ValueError( "'texture' must have the same shape " @@ -3065,6 +3066,10 @@ def __call__(self, plot): lic_data = lic_data / lic_data.max() lic_data_clip = np.clip(lic_data, self.lim[0], self.lim[1]) + mask = ~(np.isfinite(pixX) & np.isfinite(pixY)) + lic_data[mask] = np.nan + lic_data_clip[mask] = np.nan + if self.const_alpha: plot._axes.imshow( lic_data_clip, diff --git a/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py b/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py index 1799a9d550..df808993e0 100644 --- a/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py +++ b/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py @@ -270,7 +270,7 @@ def test_center_3(): @requires_module("scipy") def find_compare_maxima(expected_maxima, buf, resolution, width): buf_ndarray = buf.ndarray_view() - max_filter_buf = ndimage.filters.maximum_filter(buf_ndarray, size=5) + max_filter_buf = ndimage.maximum_filter(buf_ndarray, size=5) maxima = np.isclose(max_filter_buf, buf_ndarray, rtol=1e-09) # ignore contributions from zones of no smoothing