diff --git a/PySDM/exporters/parcel_vtk_exporter.py b/PySDM/exporters/parcel_vtk_exporter.py new file mode 100644 index 000000000..1b58732da --- /dev/null +++ b/PySDM/exporters/parcel_vtk_exporter.py @@ -0,0 +1,127 @@ +""" +VTK Exporter for Pyrcel PySDM simulations. + +This module defines `VTKExporterPyrcel`, a subclass of `PySDM.exporters.VTKExporter`, +that writes simulation outputs to VTK format using `pyevtk`. It exports product +profiles (e.g., relative humidity) as unstructured grids and particle attributes +as point clouds, along with `.pvd` collection files for time-series visualization +in ParaView. +""" + +from pyevtk.hl import unstructuredGridToVTK, pointsToVTK +from pyevtk.vtk import VtkHexahedron, VtkGroup +import numpy as np + +from PySDM.exporters import VTKExporter + + +class VTKExporterPyrcel(VTKExporter): + """ + Custom VTK exporter for Pyrcel PySDM, exporting products as grids + and attributes as point clouds for ParaView visualization. + """ + + def __init__(self, n_sd, output, mass_of_dry_air): + super().__init__() + self.output = output + self.x = np.random.random(n_sd) # pylint: disable=invalid-name + self.y = np.random.random(n_sd) # pylint: disable=invalid-name + self.z = np.random.random(n_sd) # pylint: disable=invalid-name + self.half_diagonal = [] + self.n_levels = len(self.output["products"]["z"]) + + _volume = mass_of_dry_air / output["products"]["rhod"][0] + _dz = output["products"]["z"][1] - output["products"]["z"][0] + for level in range(self.n_levels): + if level > 0: + prev_to_curr_density_ratio = ( + output["products"]["rhod"][level - 1] + / output["products"]["rhod"][level] + ) + _volume *= prev_to_curr_density_ratio + _dz = ( + output["products"]["z"][level] - output["products"]["z"][level - 1] + ) + _area = _volume / _dz + self.half_diagonal.append((2 * _area) ** 0.5) + + def write_pvd(self): + pvd_attributes = VtkGroup(self.attributes_file_path) + for key, value in self.exported_times["attributes"].items(): + pvd_attributes.addFile(key + ".vtu", sim_time=value) + pvd_attributes.save() + + pvd_products = VtkGroup(self.products_file_path) + for key, value in self.exported_times["products"].items(): + pvd_products.addFile(key + ".vtu", sim_time=value) + pvd_products.save() + + def export_products( + self, step, simulation + ): # pylint: disable=arguments-differ, too-many-locals + path = self.products_file_path + "_num" + self.add_leading_zeros(step) + self.exported_times["products"][path] = step * simulation.particulator.dt + + n_vertices = 4 * self.n_levels + x = np.zeros(n_vertices) # pylint: disable=invalid-name + y = np.zeros(n_vertices) # pylint: disable=invalid-name + z = np.zeros(n_vertices) # pylint: disable=invalid-name + conn = [0, 1, 2, 3] + ctype = np.zeros(self.n_levels - 1) + ctype[:] = VtkHexahedron.tid + + for level in range(self.n_levels): + hd = self.half_diagonal[level] + _z = self.output["products"]["z"][level] + x[i := level * 4], y[i], z[i] = -hd, -hd, _z + x[i := i + 1], y[i], z[i] = -hd, hd, _z + x[i := i + 1], y[i], z[i] = hd, hd, _z + x[i := i + 1], y[i], z[i] = hd, -hd, _z + conn += [*range(4 * (level + 1), 4 * (level + 2))] * 2 + conn = np.asarray(conn[:-4]) + offset = np.asarray(range(8, 8 * self.n_levels, 8)) + + _ = {"test_pd": np.array([44] * n_vertices)} # pointData + + _RH = self.output["products"]["S_max_percent"] # pylint: disable=invalid-name + cell_data = {"RH": np.full(shape=(len(_RH) - 1,), fill_value=np.nan)} + cell_data["RH"][:step] = (np.array(_RH[:-1] + np.diff(_RH) / 2))[:step] + unstructuredGridToVTK( + path, + x, + y, + z, + connectivity=conn, + offsets=offset, + cell_types=ctype, + cellData=cell_data, + # pointData=point_data, + # fieldData=field_data, + ) + + def export_attributes(self, step, simulation): # pylint: disable=arguments-differ + path = self.attributes_file_path + "_num" + self.add_leading_zeros(step) + self.exported_times["attributes"][path] = step * simulation.particulator.dt + + payload = {} + + for k in self.output["attributes"].keys(): + payload[k] = np.asarray(self.output["attributes"][k])[:, step].copy() + # payload["size"] = np.full(simulation.particulator.n_sd, 100.0) + if step != 0: + dz = ( + self.output["products"]["z"][step] + - self.output["products"]["z"][step - 1] + ) # pylint: disable=invalid-name + if step == 1: + self.z *= dz + else: + self.z += dz + + pointsToVTK( + path, + 2 * (self.x - 0.5) * self.half_diagonal[step], + 2 * (self.y - 0.5) * self.half_diagonal[step], + self.z, + data=payload, + ) diff --git a/examples/PySDM_examples/Pyrcel/paraview.ipynb b/examples/PySDM_examples/Pyrcel/paraview.ipynb new file mode 100644 index 000000000..94e19d933 --- /dev/null +++ b/examples/PySDM_examples/Pyrcel/paraview.ipynb @@ -0,0 +1,140 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 10, + "id": "6dcbcc6d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "c:\\Users\\strza\\Desktop\\PySDM\\examples\\PySDM_examples\\Pyrcel\\../../../PySDM/exporters/\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "from PySDM import Formulae\n", + "from PySDM.physics import si\n", + "from PySDM.initialisation.spectra import Lognormal\n", + "from PySDM.products import (\n", + " ParcelDisplacement,\n", + " AmbientTemperature,\n", + " AmbientDryAirDensity,\n", + " AmbientRelativeHumidity,\n", + " ParticleSizeSpectrumPerVolume,\n", + " ParticleVolumeVersusRadiusLogarithmSpectrum,\n", + ")\n", + "\n", + "from PySDM_examples.Pyrcel import Settings, Simulation\n", + "import sys\n", + "import os\n", + "\n", + "notebook_dir = os.path.dirname(os.path.abspath(\"__file__\"))\n", + "exporters_path = os.path.join(notebook_dir, \"../../../PySDM/exporters/\")\n", + "print(exporters_path)\n", + "sys.path.append(exporters_path)\n", + "\n", + "from parcel_vtk_exporter import VTKExporterPyrcel" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "cd6288b5", + "metadata": {}, + "outputs": [], + "source": [ + "settings = Settings(\n", + " dz=10 * si.m,\n", + " n_sd_per_mode=(25, 25),\n", + " aerosol_modes_by_kappa={\n", + " 0.54: Lognormal(norm_factor=850 / si.cm**3, m_mode=15 * si.nm, s_geom=1.6),\n", + " 1.2: Lognormal(norm_factor=10 / si.cm**3, m_mode=850 * si.nm, s_geom=1.2),\n", + " },\n", + " vertical_velocity=1.0 * si.m / si.s,\n", + " initial_pressure=775 * si.mbar,\n", + " initial_temperature=274 * si.K,\n", + " initial_relative_humidity=0.90,\n", + " displacement=1000 * si.m,\n", + " formulae=Formulae(constants={\"MAC\": 0.3}),\n", + ")\n", + "\n", + "dry_radius_bin_edges = np.logspace(\n", + " np.log10(1e-3 * si.um), np.log10(5e0 * si.um), 33, endpoint=False\n", + ")\n", + "\n", + "simulation = Simulation(\n", + " settings,\n", + " products=(\n", + " ParcelDisplacement(name=\"z\"),\n", + " AmbientRelativeHumidity(name=\"S_max_percent\", unit=\"%\", var=\"RH\"),\n", + " AmbientTemperature(name=\"T\"),\n", + " ParticleSizeSpectrumPerVolume(\n", + " name=\"dry:dN/dR\", radius_bins_edges=dry_radius_bin_edges, dry=True\n", + " ),\n", + " ParticleVolumeVersusRadiusLogarithmSpectrum(\n", + " name=\"dry:dV/dlnR\", radius_bins_edges=dry_radius_bin_edges, dry=True\n", + " ),\n", + " AmbientDryAirDensity(),\n", + " ),\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "3ff50a25", + "metadata": {}, + "outputs": [], + "source": [ + "output = simulation.run()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "15d59657", + "metadata": {}, + "outputs": [], + "source": [ + "e = VTKExporterPyrcel(n_sd=simulation.particulator.n_sd, output=output, mass_of_dry_air=simulation.particulator.environment.mass_of_dry_air)\n", + "for step in settings.output_steps:\n", + " e.export_products(step, simulation)\n", + " e.export_attributes(step, simulation)\n", + "e.write_pvd()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4a6ea0e6", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "base", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/PySDM_examples/Pyrcel/paraview_simulation.py b/examples/PySDM_examples/Pyrcel/paraview_simulation.py new file mode 100644 index 000000000..4ad22beb5 --- /dev/null +++ b/examples/PySDM_examples/Pyrcel/paraview_simulation.py @@ -0,0 +1,130 @@ +# to run pvpython script use command line: pvpython filename.py +from paraview import simple as pvs # Updated import statement + +pvs._DisableFirstRenderCameraReset() + +# get active view +renderView1 = pvs.GetActiveViewOrCreate("RenderView") + +# show data in view prod +sd_productspvd = pvs.PVDReader( + registrationName="sd_products.pvd", + FileName="C:\\Users\\strza\\Desktop\\PySDM\\examples\\PySDM_examples\\Pyrcel\\output\\sd_products.pvd", +) +sd_productspvdDisplay = pvs.Show( + sd_productspvd, renderView1, "UnstructuredGridRepresentation" +) +sd_productspvdDisplay.Representation = "Surface" +sd_productspvdDisplay.SetScalarBarVisibility(renderView1, True) +rHLUT = pvs.GetColorTransferFunction("RH") +rHLUT.RescaleTransferFunction(90.0, 101.0) + +# show data in view attr +sd_attributespvd = pvs.PVDReader( + registrationName="sd_attributes.pvd", + FileName="C:\\Users\\strza\\Desktop\\PySDM\\examples\\PySDM_examples\\Pyrcel\\output\\sd_attributes.pvd", +) +sd_attributespvdDisplay = pvs.Show( + sd_attributespvd, renderView1, "UnstructuredGridRepresentation" +) +sd_attributespvdDisplay.Representation = "Surface" +sd_attributespvdDisplay.SetScalarBarVisibility(renderView1, True) +volumeLUT = pvs.GetColorTransferFunction("volume") +volumeLUT.RescaleTransferFunction(1e-18, 1e-13) + +# Properties modified on rHLUTColorBar +rHLUTColorBar = pvs.GetScalarBar(rHLUT, renderView1) +rHLUTColorBar.LabelColor = [0.0, 0.0, 0.0] +rHLUTColorBar.DrawScalarBarOutline = 1 +rHLUTColorBar.ScalarBarOutlineColor = [0.0, 0.0, 0.0] +rHLUTColorBar.TitleColor = [0.0, 0.0, 0.0] + +# lightning and opacity sd_productspvdDisplay +sd_productspvdDisplay.Opacity = 0.4 +sd_productspvdDisplay.DisableLighting = 1 +sd_productspvdDisplay.Diffuse = 0.76 + +# Properties modified on sd_productspvdDisplay.DataAxesGrid +sd_productspvdDisplay.DataAxesGrid.GridAxesVisibility = 1 +sd_productspvdDisplay.DataAxesGrid.XTitle = "" +sd_productspvdDisplay.DataAxesGrid.YTitle = "" +sd_productspvdDisplay.DataAxesGrid.ZTitle = "" +sd_productspvdDisplay.DataAxesGrid.GridColor = [0.0, 0.0, 0.0] +sd_productspvdDisplay.DataAxesGrid.ShowGrid = 1 +sd_productspvdDisplay.DataAxesGrid.LabelUniqueEdgesOnly = 0 +sd_productspvdDisplay.DataAxesGrid.XAxisUseCustomLabels = 1 +sd_productspvdDisplay.DataAxesGrid.YAxisUseCustomLabels = 1 +sd_productspvdDisplay.DataAxesGrid.ZAxisUseCustomLabels = 1 + +# orientation axes +renderView1.OrientationAxesLabelColor = [0.0, 0.0, 0.0] +renderView1.OrientationAxesOutlineColor = [0.0, 0.0, 0.0] +renderView1.OrientationAxesXVisibility = 0 +renderView1.OrientationAxesYVisibility = 0 +renderView1.OrientationAxesZColor = [0.0, 0.0, 0.0] + +# background +renderView1.UseColorPaletteForBackground = 0 +renderView1.Background = [1.0, 1.0, 1.0] + +# rHLUT color +rHLUT.ApplyPreset("Black, Blue and White", True) +rHLUT.NanColor = [0.6666666666666666, 1.0, 1.0] + +# Properties modified on volumeLUTColorBar +volumeLUTColorBar = pvs.GetScalarBar(volumeLUT, renderView1) +volumeLUTColorBar.LabelColor = [0.0, 0.0, 0.0] +volumeLUTColorBar.DrawScalarBarOutline = 1 +volumeLUTColorBar.ScalarBarOutlineColor = [0.0, 0.0, 0.0] +volumeLUTColorBar.TitleColor = [0.0, 0.0, 0.0] + +# attr size and shape +sd_attributespvdDisplay.PointSize = 13.0 +sd_attributespvdDisplay.RenderPointsAsSpheres = 1 +sd_attributespvdDisplay.Interpolation = "PBR" + +# Properties modified on sd_attributespvdDisplay.DataAxesGrid +sd_attributespvdDisplay.DataAxesGrid.GridAxesVisibility = 1 +sd_attributespvdDisplay.DataAxesGrid.XTitle = "" +sd_attributespvdDisplay.DataAxesGrid.YTitle = "" +sd_attributespvdDisplay.DataAxesGrid.ZTitle = "Z [m]" +sd_attributespvdDisplay.DataAxesGrid.ZLabelColor = [0.0, 0.0, 0.0] +sd_attributespvdDisplay.DataAxesGrid.ZTitleColor = [0.0, 0.0, 0.0] +sd_attributespvdDisplay.DataAxesGrid.XAxisUseCustomLabels = 1 +sd_attributespvdDisplay.DataAxesGrid.YAxisUseCustomLabels = 1 + +# volumeLUT preset and scale +volumeLUT.ApplyPreset("Cold and Hot", True) +volumeLUT.MapControlPointsToLogSpace() +volumeLUT.UseLogScale = 1 +volumeLUT.NumberOfTableValues = 16 +volumeLUT.InvertTransferFunction() + +# layout +layout1 = pvs.GetLayout() +layout1.SetSize(1593, 1128) + +# current camera placement for renderView1 +renderView1.CameraPosition = [1548.945972263949, -1349.493616194682, 699.2699178747185] +renderView1.CameraFocalPoint = [ + -1.3686701146742275e-13, + 3.1755031232028544e-13, + 505.00000000000017, +] +renderView1.CameraViewUp = [ + -0.07221292769632315, + 0.06043215330151439, + 0.9955567527373153, +] +renderView1.CameraParallelScale = 534.07781536883 + +# get animation scene +animationScene1 = pvs.GetAnimationScene() +animationScene1.UpdateAnimationUsingDataTimeSteps() + +output_animation_path = "C:\\Users\\strza\\Desktop\\PySDM\\examples\\PySDM_examples\\Pyrcel\\output_animation.avi" # Change the extension as needed +output_screenshot_path = ( + "C:\\Users\\strza\\Desktop\\PySDM\\examples\\PySDM_examples\\Pyrcel\\last_frame.png" +) +pvs.SaveAnimation(output_animation_path, renderView1, FrameRate=15) +pvs.SaveScreenshot(output_screenshot_path, renderView1) diff --git a/examples/PySDM_examples/Pyrcel/simulation.py b/examples/PySDM_examples/Pyrcel/simulation.py index aa0f92b70..e631c347b 100644 --- a/examples/PySDM_examples/Pyrcel/simulation.py +++ b/examples/PySDM_examples/Pyrcel/simulation.py @@ -27,7 +27,7 @@ def __init__( initial_water_vapour_mixing_ratio=settings.initial_vapour_mixing_ratio, T0=settings.initial_temperature, w=settings.vertical_velocity, - mass_of_dry_air=44 * si.kg, + mass_of_dry_air=66666 * si.kg, ), ) builder.add_dynamic(AmbientThermodynamics()) diff --git a/examples/PySDM_examples/Strzabala_engineering_thesis_2025/calc2.py b/examples/PySDM_examples/Strzabala_engineering_thesis_2025/calc2.py new file mode 100644 index 000000000..ade9b36f8 --- /dev/null +++ b/examples/PySDM_examples/Strzabala_engineering_thesis_2025/calc2.py @@ -0,0 +1,343 @@ +#!/usr/bin/env pvpython +import argparse +from collections import namedtuple +import pathlib + +from paraview import simple as pvs # pylint: disable=import-error + +pvs._DisableFirstRenderCameraReset() + + +def cli_using_argparse(argp): + argp.add_argument("product_path", help="path to pvd products file") + argp.add_argument("attributes_path", help=" path to pvd attributes file") + argp.add_argument("output_path", help="path where to write output files") + argp.add_argument( + "--mode", + choices=["light", "dark"], + default="light", + help="Choose 'light' or 'dark' mode.", + ) + argp.add_argument( + "--multiplicity_preset", + default="Inferno (matplotlib)", + help="Preset for multiplicity", + ) + argp.add_argument( + "--multiplicity_logscale", + action="store_false", + help="Use log scale for multiplicity", + ) + argp.add_argument( + "--effectiveradius_preset", + default="Black, Blue and White", + help="Preset for effectiveradius", + ) + argp.add_argument( + "--effectiveradius_logscale", + action="store_false", + help="Use log scale for effectiveradius", + ) + argp.add_argument( + "--effectiveradius_nan_color", + nargs=3, + type=float, + default=[0.666, 0.333, 1.0], + help="Nan color in RGB format for effectiveradius", + ) + argp.add_argument( + "--sd_products_opacity", type=float, default=0.0, help="Opacity for sd_products" + ) + argp.add_argument( + "--calculator1_opacity", + type=float, + default=0.0, + help="Opacity for calculator1", + ) + argp.add_argument( + "--sd_attributes_opacity", + type=float, + default=0.0, + help="Opacity for sd_attributes", + ) + argp.add_argument( + "--animation_size", + nargs=2, + type=int, + default=[800, 800], + help="Animation size [x,y]", + ) + argp.add_argument( + "--animationframename", + type=str, + help="Name of the file with animation last frame", + ) + argp.add_argument( + "--animationname", + type=str, + help="Name of the file with animation", + ) + argp.add_argument( + "--framerate", type=int, help="Number of frame rates.", default=15 + ) + + +ap = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) +cli_using_argparse(ap) + +args = ap.parse_args() +sd_productspvd = pvs.OpenDataFile(args.product_path) +sd_attributespvd = pvs.OpenDataFile(args.attributes_path) + + +setup = { + "renderView1": pvs.GetActiveViewOrCreate("RenderView"), + "sd_attributespvdDisplay": pvs.GetDisplayProperties( + sd_attributespvd, view=pvs.GetActiveViewOrCreate("RenderView") + ), + "effectiveradiusLUT": pvs.GetColorTransferFunction("effectiveradius"), + "sd_productspvdDisplay": pvs.GetDisplayProperties( + sd_productspvd, view=pvs.GetActiveViewOrCreate("RenderView") + ), + "color": [1, 1, 1] if args.mode == "dark" else [0, 0, 0], + "inverted_color": [0.129, 0.145, 0.161] if args.mode == "dark" else [1, 1, 1], +} + +setup = namedtuple("Setup", setup.keys())(**setup) + +setup.effectiveradiusLUT.RescaleTransferFunction(0.1380997175392798, 207.7063518856934) +materialLibrary1 = pvs.GetMaterialLibrary() +setup.renderView1.Update() + + +def create_new_calculator( + calcinput, + representation, + function, + color_by1, + color_by2, + color_by3, + scalar_coloring=False, + hide=False, + *, + y, + registrationame, +): + calculator = pvs.Calculator(registrationName=registrationame, Input=calcinput) + display = pvs.Show(calculator, y.renderView1, representation) + calculator.Function = function + y.renderView1.Update() + if scalar_coloring is True: + pvs.ColorBy(display, (color_by1, color_by2, color_by3)) + if hide is True: + pvs.Hide(calculator, y.renderView1) + return y.renderView1.Update() + + +def scalar_bar(name, *, y, erLUT): + calculator1Display = pvs.Show( + calculator1, y.renderView1, "UnstructuredGridRepresentation" + ) + calculator1Display.SetScalarBarVisibility(y.renderView1, True) + scalarBar = pvs.GetScalarBar(erLUT.effectiveradiusLUT, y.renderView1) + scalarBar.ComponentTitle = "" + scalarBar.Title = name + scalarBar.TitleFontSize = 25 + scalarBar.LabelFontSize = 25 + scalarBar.LabelColor = setup.color + scalarBar.TitleColor = setup.color + pvs.Hide(scalarBar) + y.renderView1.Update() + + +def create_glyph( + registration_name, put, scale_array1, scale_array2, color_by=False, *, y +): + glyph = pvs.Glyph(registrationName=registration_name, Input=put, GlyphType="Arrow") + glyphDisplay = pvs.Show(glyph, y.renderView1, "GeometryRepresentation") + glyphDisplay.Representation = "Surface" + glyph.ScaleArray = [scale_array1, scale_array2] + glyph.ScaleFactor = 100 + glyphDisplay.SetScalarBarVisibility(y.renderView1, True) + multiplicityLUT = pvs.GetColorTransferFunction("multiplicity") + multiplicityLUTColorBar = pvs.GetScalarBar(multiplicityLUT, y.renderView1) + multiplicityLUTColorBar.Position = [0.8, 0.5] + multiplicityLUTColorBar.TitleFontSize = 25 + multiplicityLUTColorBar.LabelFontSize = 25 + multiplicityLUTColorBar.LabelColor = setup.color + multiplicityLUTColorBar.TitleColor = setup.color + if color_by is True: + glyphDisplay.ColorArrayName = ["POINTS", ""] + pvs.ColorBy(glyphDisplay, None) + y.renderView1.Update() + + +def apply_presets_logscale_opacity_and_update(*, y, attdisplay, erLUT, proddisplay): + multiplicityLUT = pvs.GetColorTransferFunction("multiplicity") + multiplicityLUT.RescaleTransferFunction(19951.0, 50461190157.0) + calculator1Display = pvs.Show( + calculator1, y.renderView1, "UnstructuredGridRepresentation" + ) + multiplicityLUT.ApplyPreset(args.multiplicity_preset, True) + if args.multiplicity_logscale: + multiplicityLUT.MapControlPointsToLogSpace() + multiplicityLUT.UseLogScale = 1 + else: + multiplicityLUT.MapControlPointsToLinearSpace() + multiplicityLUT.UseLogScale = 0 + + erLUT.effectiveradiusLUT.ApplyPreset(args.effectiveradius_preset, True) + if args.effectiveradius_logscale: + erLUT.effectiveradiusLUT.MapControlPointsToLogSpace() + erLUT.effectiveradiusLUT.UseLogScale = 1 + else: + erLUT.effectiveradiusLUT.MapControlPointsToLinearSpace() + erLUT.effectiveradiusLUT.UseLogScale = 0 + + erLUT.effectiveradiusLUT.NanColor = args.effectiveradius_nan_color + # proddisplay.sd_productspvdDisplay.SetRepresentationType("Surface With Edges") + proddisplay.sd_productspvdDisplay.Opacity = args.sd_products_opacity + proddisplay.sd_productspvdDisplay.SetScalarBarVisibility(y.renderView1, False) + calculator1Display.Opacity = args.calculator1_opacity + attdisplay.sd_attributespvdDisplay.Opacity = args.sd_attributes_opacity + y.renderView1.Update() + + +def get_layout(*, y): + pvs.SetViewProperties( + Background=setup.inverted_color, UseColorPaletteForBackground=0 + ) + pvs.Render(setup.renderView1) + layout1 = pvs.GetLayout() + layout1.SetSize(args.animation_size) + layout1.PreviewMode = args.animation_size + y.renderView1.Update() + + +def set_current_camera_placement(*, y): + y.renderView1.InteractionMode = "2D" + y.renderView1.CameraPosition = [ + 836, + 677, + -4098, + ] + y.renderView1.CameraFocalPoint = [636, 1030, 0.0] + y.renderView1.CameraViewUp = [1.0, 0.0, 0.0] + y.renderView1.CameraParallelScale = 1560 + y.renderView1.Update() + + +def axes_settings(*, view): + # setup.renderView1.Background = [1,0.5,0.2] + view.CenterAxesVisibility = True + view.OrientationAxesVisibility = False + axesGrid = view.AxesGrid + axesGrid.Visibility = True + axesGrid.XTitle = "Z [m]" + axesGrid.YTitle = "X [m]" + + axesGrid.XAxisUseCustomLabels = True + axesGrid.XAxisLabels = [300, 600, 900, 1200] + axesGrid.YAxisUseCustomLabels = True + axesGrid.YAxisLabels = [300, 600, 900, 1200] + + axesGrid.XTitleFontSize = 30 + axesGrid.XLabelFontSize = 30 + axesGrid.YTitleFontSize = 30 + axesGrid.YLabelFontSize = 30 + + axesGrid.XTitleColor = setup.color + axesGrid.XLabelColor = setup.color + axesGrid.YTitleColor = setup.color + axesGrid.YLabelColor = setup.color + axesGrid.GridColor = [0.1, 0.1, 0.1] + view.CenterAxesVisibility = False + view.Update() + + +def time_annotation(*, y): + time = pvs.AnnotateTimeFilter( + guiName="AnnotateTimeFilter1", Scale=1 / 60, Format="Time:{time:g}min" + ) + timedisplay = pvs.Show(time, y.renderView1) + timedisplay.FontSize = 25 + timedisplay.WindowLocation = "Any Location" + timedisplay.FontSize = 30 + timedisplay.Position = [0.4, 0.9] + timedisplay.Color = setup.color + y.renderView1.Update() + + +def text(text_in, position_y, *, view): + sentence = pvs.Text() + sentence.Text = text_in + textDisplay = pvs.Show(sentence, view) + textDisplay.Color = setup.color + textDisplay.WindowLocation = "Any Location" + textDisplay.FontSize = 28 + textDisplay.Position = [0.17, position_y] + + +def last_anim_frame(animation_frame_name): + time_steps = sd_productspvd.TimestepValues + last_time = time_steps[90] + setup.renderView1.ViewTime = last_time + for reader in (sd_productspvd,): + reader.UpdatePipeline(last_time) + pvs.ExportView( + filename=str(pathlib.Path(args.output_path) / animation_frame_name), + view=setup.renderView1, + Rasterize3Dgeometry=False, + GL2PSdepthsortmethod="BSP sorting (slow, best)", + ) + pvs.RenderAllViews() + + +calculator1 = create_new_calculator( + sd_attributespvd, + "UnstructuredGridRepresentation", + '"relative fall velocity"*(-iHat)', + "None", + "None", + "None", + y=setup, + registrationame="Calculator1", +) +# scalar_bar("effective radius [um]", y=setup, erLUT=setup) +glyph1 = create_glyph( + "Glyph1", calculator1, "POINTS", "relative fall velocity", y=setup +) +pvs.Hide(glyph1) +"""calculator2 = create_new_calculator( + sd_productspvd, + "StructuredGridRepresentation", + "cx*jHat+cy*iHat", + "CELLS", + "Result", + "Magnitude", + True, + True, + y=setup, + registrationame="Calculator2", +)""" +# pvs.Hide(calculator2) +apply_presets_logscale_opacity_and_update( + y=setup, attdisplay=setup, erLUT=setup, proddisplay=setup +) +# glyph2 = create_glyph("Glyph2", calculator2, "CELLS", "Result", True, y=setup) +# pvs.Hide(glyph2) +get_layout(y=setup) +set_current_camera_placement(y=setup) +axes_settings(view=setup.renderView1) +time_annotation(y=setup) +if args.animationframename is not None: + last_anim_frame(animation_frame_name=args.animationframename) +scene = pvs.GetAnimationScene() +scene.UpdateAnimationUsingDataTimeSteps() +pvs.Render(setup.renderView1) +pvs.SaveAnimation( + str(pathlib.Path(args.output_path) / args.animationname), + setup.renderView1, + FrameRate=args.framerate, +) +pvs.RenderAllViews() diff --git a/examples/PySDM_examples/Strzabala_engineering_thesis_2025/no_glyphs.py b/examples/PySDM_examples/Strzabala_engineering_thesis_2025/no_glyphs.py new file mode 100644 index 000000000..803819ecb --- /dev/null +++ b/examples/PySDM_examples/Strzabala_engineering_thesis_2025/no_glyphs.py @@ -0,0 +1,341 @@ +#!/usr/bin/env pvpython +import argparse +from collections import namedtuple +import pathlib + +from paraview import simple as pvs # pylint: disable=import-error + +pvs._DisableFirstRenderCameraReset() + + +def cli_using_argparse(argp): + argp.add_argument("product_path", help="path to pvd products file") + argp.add_argument("attributes_path", help=" path to pvd attributes file") + argp.add_argument("output_path", help="path where to write output files") + argp.add_argument( + "--mode", + choices=["light", "dark"], + default="dark", + help="Choose 'light' or 'dark' mode.", + ) + argp.add_argument( + "--multiplicity_preset", + default="Inferno (matplotlib)", + help="Preset for multiplicity", + ) + argp.add_argument( + "--multiplicity_logscale", + action="store_false", + help="Use log scale for multiplicity", + ) + argp.add_argument( + "--effectiveradius_preset", + default="Black, Blue and White", + help="Preset for effectiveradius", + ) + argp.add_argument( + "--effectiveradius_logscale", + action="store_false", + help="Use log scale for effectiveradius", + ) + argp.add_argument( + "--effectiveradius_nan_color", + nargs=3, + type=float, + default=[0.666, 0.333, 1.0], + help="Nan color in RGB format for effectiveradius", + ) + argp.add_argument( + "--sd_products_opacity", type=float, default=0.9, help="Opacity for sd_products" + ) + argp.add_argument( + "--calculator1_opacity", + type=float, + default=0.19, + help="Opacity for calculator1", + ) + argp.add_argument( + "--sd_attributes_opacity", + type=float, + default=0.77, + help="Opacity for sd_attributes", + ) + argp.add_argument( + "--animation_size", + nargs=2, + type=int, + default=[800, 800], + help="Animation size [x,y]", + ) + argp.add_argument( + "--animationframename", + type=str, + help="Name of the file with animation last frame", + ) + argp.add_argument( + "--animationname", + type=str, + help="Name of the file with animation", + ) + argp.add_argument( + "--framerate", type=int, help="Number of frame rates.", default=15 + ) + + +ap = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) +cli_using_argparse(ap) + +args = ap.parse_args() +sd_productspvd = pvs.OpenDataFile(args.product_path) +sd_attributespvd = pvs.OpenDataFile(args.attributes_path) + + +setup = { + "renderView1": pvs.GetActiveViewOrCreate("RenderView"), + "sd_attributespvdDisplay": pvs.GetDisplayProperties( + sd_attributespvd, view=pvs.GetActiveViewOrCreate("RenderView") + ), + "effectiveradiusLUT": pvs.GetColorTransferFunction("effectiveradius"), + "sd_productspvdDisplay": pvs.GetDisplayProperties( + sd_productspvd, view=pvs.GetActiveViewOrCreate("RenderView") + ), + "color": [1, 1, 1] if args.mode == "dark" else [0, 0, 0], + "inverted_color": [0.129, 0.145, 0.161] if args.mode == "dark" else [1, 1, 1], +} + +setup = namedtuple("Setup", setup.keys())(**setup) + +setup.effectiveradiusLUT.RescaleTransferFunction(0.1380997175392798, 207.7063518856934) +materialLibrary1 = pvs.GetMaterialLibrary() +setup.renderView1.Update() + + +def create_new_calculator( + calcinput, + representation, + function, + color_by1, + color_by2, + color_by3, + scalar_coloring=False, + hide=False, + *, + y, + registrationame, +): + calculator = pvs.Calculator(registrationName=registrationame, Input=calcinput) + display = pvs.Show(calculator, y.renderView1, representation) + calculator.Function = function + y.renderView1.Update() + if scalar_coloring is True: + pvs.ColorBy(display, (color_by1, color_by2, color_by3)) + if hide is True: + pvs.Hide(calculator, y.renderView1) + return y.renderView1.Update() + + +def scalar_bar(name, *, y, erLUT): + calculator1Display = pvs.Show( + calculator1, y.renderView1, "UnstructuredGridRepresentation" + ) + calculator1Display.SetScalarBarVisibility(y.renderView1, True) + scalarBar = pvs.GetScalarBar(erLUT.effectiveradiusLUT, y.renderView1) + scalarBar.ComponentTitle = "" + scalarBar.Title = name + scalarBar.TitleFontSize = 25 + scalarBar.LabelFontSize = 25 + scalarBar.LabelColor = setup.color + scalarBar.TitleColor = setup.color + y.renderView1.Update() + + +def create_glyph( + registration_name, put, scale_array1, scale_array2, color_by=False, *, y +): + glyph = pvs.Glyph(registrationName=registration_name, Input=put, GlyphType="Arrow") + glyphDisplay = pvs.Show(glyph, y.renderView1, "GeometryRepresentation") + glyphDisplay.Representation = "Surface" + glyph.ScaleArray = [scale_array1, scale_array2] + glyph.ScaleFactor = 100 + glyphDisplay.SetScalarBarVisibility(y.renderView1, True) + multiplicityLUT = pvs.GetColorTransferFunction("multiplicity") + multiplicityLUTColorBar = pvs.GetScalarBar(multiplicityLUT, y.renderView1) + multiplicityLUTColorBar.TitleFontSize = 25 + multiplicityLUTColorBar.LabelFontSize = 25 + multiplicityLUTColorBar.LabelColor = setup.color + multiplicityLUTColorBar.TitleColor = setup.color + if color_by is True: + glyphDisplay.ColorArrayName = ["POINTS", ""] + pvs.ColorBy(glyphDisplay, None) + y.renderView1.Update() + + +def apply_presets_logscale_opacity_and_update(*, y, attdisplay, erLUT, proddisplay): + multiplicityLUT = pvs.GetColorTransferFunction("multiplicity") + multiplicityLUT.RescaleTransferFunction(19951.0, 50461190157.0) + calculator1Display = pvs.Show( + calculator1, y.renderView1, "UnstructuredGridRepresentation" + ) + multiplicityLUT.ApplyPreset(args.multiplicity_preset, True) + if args.multiplicity_logscale: + multiplicityLUT.MapControlPointsToLogSpace() + multiplicityLUT.UseLogScale = 1 + else: + multiplicityLUT.MapControlPointsToLinearSpace() + multiplicityLUT.UseLogScale = 0 + + erLUT.effectiveradiusLUT.ApplyPreset(args.effectiveradius_preset, True) + if args.effectiveradius_logscale: + erLUT.effectiveradiusLUT.MapControlPointsToLogSpace() + erLUT.effectiveradiusLUT.UseLogScale = 1 + else: + erLUT.effectiveradiusLUT.MapControlPointsToLinearSpace() + erLUT.effectiveradiusLUT.UseLogScale = 0 + + erLUT.effectiveradiusLUT.NanColor = args.effectiveradius_nan_color + + proddisplay.sd_productspvdDisplay.SetRepresentationType("Surface With Edges") + proddisplay.sd_productspvdDisplay.Opacity = args.sd_products_opacity + calculator1Display.Opacity = args.calculator1_opacity + attdisplay.sd_attributespvdDisplay.Opacity = args.sd_attributes_opacity + + y.renderView1.Update() + + +def get_layout(*, y): + pvs.SetViewProperties( + Background=setup.inverted_color, UseColorPaletteForBackground=0 + ) + pvs.Render(setup.renderView1) + layout1 = pvs.GetLayout() + layout1.SetSize(args.animation_size) + layout1.PreviewMode = args.animation_size + y.renderView1.Update() + + +def set_current_camera_placement(*, y): + y.renderView1.InteractionMode = "2D" + y.renderView1.CameraPosition = [ + 836, + 677, + -4098, + ] + y.renderView1.CameraFocalPoint = [636, 1030, 0.0] + y.renderView1.CameraViewUp = [1.0, 0.0, 0.0] + y.renderView1.CameraParallelScale = 1560 + y.renderView1.Update() + + +def axes_settings(*, view): + # setup.renderView1.Background = [1,0.5,0.2] + view.CenterAxesVisibility = True + view.OrientationAxesVisibility = False + axesGrid = view.AxesGrid + axesGrid.Visibility = True + axesGrid.XTitle = "Z [m]" + axesGrid.YTitle = "X [m]" + + axesGrid.XAxisUseCustomLabels = True + axesGrid.XAxisLabels = [300, 600, 900, 1200] + axesGrid.YAxisUseCustomLabels = True + axesGrid.YAxisLabels = [300, 600, 900, 1200] + + axesGrid.XTitleFontSize = 30 + axesGrid.XLabelFontSize = 30 + axesGrid.YTitleFontSize = 30 + axesGrid.YLabelFontSize = 30 + + axesGrid.XTitleColor = setup.color + axesGrid.XLabelColor = setup.color + axesGrid.YTitleColor = setup.color + axesGrid.YLabelColor = setup.color + axesGrid.GridColor = [0.1, 0.1, 0.1] + view.CenterAxesVisibility = False + view.Update() + + +def time_annotation(*, y): + time = pvs.AnnotateTimeFilter( + guiName="AnnotateTimeFilter1", Scale=1 / 60, Format="Time:{time:g}min" + ) + timedisplay = pvs.Show(time, y.renderView1) + timedisplay.FontSize = 25 + timedisplay.WindowLocation = "Any Location" + timedisplay.FontSize = 30 + timedisplay.Position = [0.4, 0.9] + timedisplay.Color = setup.color + y.renderView1.Update() + + +def text(text_in, position_y, *, view): + sentence = pvs.Text() + sentence.Text = text_in + textDisplay = pvs.Show(sentence, view) + textDisplay.Color = setup.color + textDisplay.WindowLocation = "Any Location" + textDisplay.FontSize = 28 + textDisplay.Position = [0.17, position_y] + + +def last_anim_frame(animation_frame_name): + time_steps = sd_productspvd.TimestepValues + last_time = time_steps[90] + setup.renderView1.ViewTime = last_time + for reader in (sd_productspvd,): + reader.UpdatePipeline(last_time) + pvs.ExportView( + filename=str(pathlib.Path(args.output_path) / animation_frame_name), + view=setup.renderView1, + Rasterize3Dgeometry=False, + GL2PSdepthsortmethod="BSP sorting (slow, best)", + ) + pvs.RenderAllViews() + + +calculator1 = create_new_calculator( + sd_attributespvd, + "UnstructuredGridRepresentation", + '"relative fall velocity"*(-iHat)', + "None", + "None", + "None", + y=setup, + registrationame="Calculator1", +) +scalar_bar("effective radius [um]", y=setup, erLUT=setup) +glyph1 = create_glyph( + "Glyph1", calculator1, "POINTS", "relative fall velocity", y=setup +) +pvs.Hide(glyph1) +calculator2 = create_new_calculator( + sd_productspvd, + "StructuredGridRepresentation", + "cx*jHat+cy*iHat", + "CELLS", + "Result", + "Magnitude", + True, + True, + y=setup, + registrationame="Calculator2", +) +apply_presets_logscale_opacity_and_update( + y=setup, attdisplay=setup, erLUT=setup, proddisplay=setup +) +glyph2 = create_glyph("Glyph2", calculator2, "CELLS", "Result", True, y=setup) +pvs.Hide(glyph2) +get_layout(y=setup) +set_current_camera_placement(y=setup) +axes_settings(view=setup.renderView1) +time_annotation(y=setup) +if args.animationframename is not None: + last_anim_frame(animation_frame_name=args.animationframename) +scene = pvs.GetAnimationScene() +scene.UpdateAnimationUsingDataTimeSteps() +pvs.Render(setup.renderView1) +pvs.SaveAnimation( + str(pathlib.Path(args.output_path) / args.animationname), + setup.renderView1, + FrameRate=args.framerate, +) +pvs.RenderAllViews()