diff --git a/iga/use_cases/toyota_3d_sbm/README.md b/iga/use_cases/toyota_3d_sbm/README.md index 0e6a9384..a3a168aa 100644 --- a/iga/use_cases/toyota_3d_sbm/README.md +++ b/iga/use_cases/toyota_3d_sbm/README.md @@ -1,128 +1,13 @@ -# 3D Toyota SBM Example +# 3D Toyota SBM Examples -This use case demonstrates a 3D immersed-body fluid analysis with the Surrogate Boundary Method (SBM) in the Kratos IGA Application. The Toyota body is not discretized with a body-fitted fluid mesh. Instead, the car surface is provided as an immersed skin in `toyota_immersed_surface_sbm.mdpa`, while the fluid is solved on a structured Cartesian background NURBS volume generated in the JSON settings. +This directory contains the shared geometry for the Toyota SBM examples together with two split use cases: -## SBM Philosophy +- [steady_stokes/README.md](/home/nantonelli/Examples/iga/use_cases/toyota_3d_sbm/steady_stokes/README.md): steady-state-style 3D Stokes example +- [transient_navier_stokes/README.md](/home/nantonelli/Examples/iga/use_cases/toyota_3d_sbm/transient_navier_stokes/README.md): transient 3D Navier-Stokes example -The main idea of SBM is to decouple the geometry description from the analysis mesh: +The shared files stored in this parent folder are: -- the immersed body is described only by a surface mesh -- the analysis domain is a simple Cartesian background volume -- boundary conditions are imposed on a surrogate boundary built from the background mesh, rather than on a body-fitted mesh +- `toyota_immersed_surface_sbm.mdpa`: common immersed Toyota surface mesh used by both subcases +- `FluidMaterials.json`: common material definition -For this kind of problem, that is useful because the geometry can come from an external triangulated surface, for example an STL-like workflow prepared outside Kratos and then exported to `.mdpa`. The analysis mesh can then be changed independently by editing the background box and its resolution directly in the parameter files. - -In this example, the no-slip condition on the immersed Toyota surface is imposed weakly on `IgaModelPart.SBM_Support_inner` with `SbmFluidConditionDirichlet`. The formulation includes skew-symmetric Nitsche terms, and the example is configured in penalty-free mode on the immersed surrogate boundary through `PENALTY_FACTOR: 0.0` in `FluidMaterials.json`. - -## Geometry and Modelers - -Both `ProjectParameters_3D_fluid.json` and `ProjectParameters_3D_fluid_steady_state.json` follow the same three-step geometry pipeline: - -1. `import_mdpa_modeler` imports `toyota_immersed_surface_sbm.mdpa` into `initial_skin_model_part_in`. -2. `NurbsGeometryModelerSbm` creates the background 3D Cartesian NURBS volume. -3. `IgaModelerSbm` builds the fluid domain and the fitted/surrogate support conditions needed by the SBM formulation. - -The main model parts created by the setup are: - -- `IgaModelPart.FluidDomain`: background fluid volume -- `IgaModelPart.SBM_Support_inner`: surrogate boundary associated with the immersed Toyota body -- `IgaModelPart.SBM_Support_outer`: outer wall support surfaces of the background box -- `IgaModelPart.Support_inlet`: inlet support surface -- `IgaModelPart.SupportOuterPressure`: outlet pressure support surface - -## What You Can Change in the JSON - -The most important SBM geometry controls are in the `NurbsGeometryModelerSbm` block: - -- `input_filename`: selects the immersed surface mesh to import -- `lower_point_xyz` and `upper_point_xyz`: define the background Cartesian box -- `polynomial_order`: defines the polynomial degree of the background NURBS volume -- `number_of_knot_spans`: defines the Cartesian resolution in each direction -- `lambda_outer`, `lambda_inner`, `number_of_inner_loops`: control the SBM surrogate-boundary construction - -In practice, if you want to analyze another immersed geometry, the usual workflow is: - -1. generate or export a triangulated surface externally -2. convert it to an `.mdpa` skin mesh -3. replace `toyota_immersed_surface_sbm.mdpa` -4. adapt the background box and the number of knot spans in the JSON files - -## Physics in This Folder - -This folder contains two related fluid setups: - -- `ProjectParameters_3D_fluid.json`: transient monolithic 3D Navier-Stokes case with `NavierStokesElement` and a ramped inlet velocity -- `ProjectParameters_3D_fluid_steady_state.json`: steady-state-style 3D case using `StokesElement` and a constant inlet velocity - -The boundary treatment is: - -- inlet velocity prescribed on `IgaModelPart.Support_inlet` -- zero outlet pressure prescribed on `IgaModelPart.SupportOuterPressure` -- outer-box support conditions imposed weakly on `IgaModelPart.SBM_Support_outer` -- immersed no-slip condition imposed weakly on `IgaModelPart.SBM_Support_inner` - -The material file uses a 3D Newtonian constitutive law. In the current setup: - -- `IgaModelPart` uses `DENSITY = 1.0`, `DYNAMIC_VISCOSITY = 1e-3`, and `PENALTY_FACTOR = 1e3` for the fitted support terms on the background box -- `IgaModelPart.SBM_Support_inner` uses `PENALTY_FACTOR = 0.0` for the penalty-free immersed-boundary configuration - -## Files - -- `ProjectParameters_3D_fluid.json`: transient fluid setup -- `ProjectParameters_3D_fluid_steady_state.json`: steady-state fluid setup -- `FluidMaterials.json`: 3D Newtonian material data and SBM penalty settings -- `toyota_immersed_surface_sbm.mdpa`: immersed Toyota skin used by the SBM modeler -- `run_and_post_nurbs_time.py`: runs the transient case and generates postprocessed GIFs -- `run_and_post_nurbs_steady_state.py`: runs the steady case and generates final plots -- `plot_conditions.py`: visualizes the surrogate faces created by the geometry/modeler setup - -## Python Dependencies - -The post-processing scripts use: - -- `numpy` -- `matplotlib` -- `imageio` for GIF writing in the transient script - -If a `latex` executable is available, the plotting scripts can use Computer Modern style math text. Otherwise they fall back to standard matplotlib text rendering. - -## Run - -Transient case: - -```bash -cd /home/nantonelli/Examples/iga/use_cases/toyota_3d_sbm -python3 run_and_post_nurbs_time.py -``` - -Steady-state case: - -```bash -cd /home/nantonelli/Examples/iga/use_cases/toyota_3d_sbm -python3 run_and_post_nurbs_steady_state.py -``` - -Geometry / surrogate-condition plot: - -```bash -cd /home/nantonelli/Examples/iga/use_cases/toyota_3d_sbm -python3 plot_conditions.py -``` - -## Outputs - -Transient script: - -- `particles_trails.gif` -- `cut_contour_x_lt_0.gif` -- `plane_x_eq_0_contour.gif` - -Steady-state script: - -- `steady_3d_velocity_pressure.png` -- `steady_plane_x0_velocity.png` -- `steady_plane_x0_pressure.png` - -Geometry helper: - -- `surrogate_faces_plot.png` +The split layout avoids duplicating the immersed `.mdpa` file while keeping the steady and transient examples documented separately. diff --git a/iga/use_cases/toyota_3d_sbm/ProjectParameters_3D_fluid_steady_state.json b/iga/use_cases/toyota_3d_sbm/steady_stokes/ProjectParameters_3D_fluid_steady_state.json similarity index 83% rename from iga/use_cases/toyota_3d_sbm/ProjectParameters_3D_fluid_steady_state.json rename to iga/use_cases/toyota_3d_sbm/steady_stokes/ProjectParameters_3D_fluid_steady_state.json index 0c79dbc9..887ea503 100644 --- a/iga/use_cases/toyota_3d_sbm/ProjectParameters_3D_fluid_steady_state.json +++ b/iga/use_cases/toyota_3d_sbm/steady_stokes/ProjectParameters_3D_fluid_steady_state.json @@ -6,7 +6,6 @@ "echo_level": 1, "parallel_type": "OpenMP", "start_time": 0.0, - // "end_time": 9.0 "end_time": 5.0 }, @@ -17,7 +16,7 @@ "model_part_name": "IgaModelPart", "domain_size": 3, "model_import_settings": { "input_type": "use_input_model_part" }, - "material_import_settings": { "materials_filename": "FluidMaterials.json" }, + "material_import_settings": { "materials_filename": "../FluidMaterials.json" }, "linear_solver_settings": { "solver_type": "bicgstab", "tolerance": 1.0e-14, "max_iteration": 1000, "scaling": true, "preconditioner_type": "ilu0"}, // "linear_solver_settings":{ // // "solver_type": "LinearSolversApplication.sparse_lu" @@ -44,7 +43,7 @@ "modeler_name" : "import_mdpa_modeler", "kratos_module" : "KratosMultiphysics", "Parameters" : { - "input_filename" : "toyota_immersed_surface_sbm", + "input_filename" : "../toyota_immersed_surface_sbm", "model_part_name" : "initial_skin_model_part_in" } }, @@ -139,8 +138,6 @@ { "variable_name": "BODY_FORCE", "value": ["0", "0", "0"] - // NAVIER STOKES - // "value": ["(cosh(x)*sinh(x)*cosh(y)*cosh(y)*cosh(z)*cosh(z) + cosh(x)*sinh(x)*sinh(y)*cosh(z)*sinh(z-y) - cosh(x)*sinh(x)*sinh(y)*cosh(y)*cosh(z)*sinh(z)) + cosh(x) - 3*cosh(x)*cosh(y)*cosh(z)", "(cosh(x)*cosh(x)*cosh(y)*cosh(z)*sinh(z-y) - sinh(x)*sinh(x)*sinh(z-y)*cosh(z-y) - sinh(x)*sinh(x)*sinh(y)*cosh(z)*cosh(z-y)) + sinh(y) - 3*sinh(x)*sinh(z-y)", "(-cosh(x)*cosh(x)*cosh(y)*sinh(y)*cosh(z)*cosh(z) - sinh(x)*sinh(x)*sinh(z-y)*cosh(y)*cosh(z) + sinh(x)*sinh(x)*sinh(y)*sinh(y)*cosh(z)*sinh(z)) + cosh(z) + 3*sinh(x)*sinh(y)*cosh(z)"] } ] }, @@ -154,7 +151,6 @@ "0.0", "0.0" ] - // "value": ["cosh(x)*cosh(y)*cosh(z)","sinh(x)*sinh(z-y)","-sinh(x)*sinh(y)*cosh(z)"] }] }, { @@ -167,7 +163,6 @@ "0.0", "2.0" ] - // "value": ["cosh(x)*cosh(y)*cosh(z)","sinh(x)*sinh(z-y)","-sinh(x)*sinh(y)*cosh(z)"] }] }, { @@ -176,7 +171,6 @@ { "variable_name": "PRESSURE", "value": "0.0" - // "value": "sinh(x) + cosh(y) + sinh(z)" }, { "variable_name": "NORMAL_STRESS", @@ -203,7 +197,6 @@ "0.0", "0.0" ] - // "value": ["cosh(x)*cosh(y)*cosh(z)","sinh(x)*sinh(z-y)","-sinh(x)*sinh(y)*cosh(z)"] } } ], diff --git a/iga/use_cases/toyota_3d_sbm/steady_stokes/README.md b/iga/use_cases/toyota_3d_sbm/steady_stokes/README.md new file mode 100644 index 00000000..b0542f9a --- /dev/null +++ b/iga/use_cases/toyota_3d_sbm/steady_stokes/README.md @@ -0,0 +1,335 @@ +# 3D Toyota SBM Steady Stokes Example + +Single-patch 3D SBM/IGA [1,2] example for a steady-state-style incompressible Stokes flow around an immersed Toyota-like car body. The body is not surrounded by a body-fitted fluid mesh. Instead, the geometry is provided through the shared immersed surface `../toyota_immersed_surface_sbm.mdpa`, while the fluid is solved on a structured Cartesian background NURBS volume generated directly from the JSON settings. + +This folder contains the steady Stokes variant of the Toyota SBM example: + +- steady-state-style 3D Stokes case in `ProjectParameters_3D_fluid_steady_state.json` +- post-processing driver `run_and_post_nurbs_steady_state.py` +- geometry helper `plot_conditions.py` + +The example uses: + +- the immersed Toyota surface mesh in `.mdpa` format +- a single Cartesian background NURBS volume created with `NurbsGeometryModelerSbm` +- `IgaModelerSbm` to generate the fluid domain and the surrogate support boundaries +- SBM/IGA to impose the immersed no-slip condition without body-fitted remeshing +- weak boundary conditions on the outer box walls, inlet, outlet, and immersed body + +The aim of the example is methodological: it demonstrates how a 3D external flow problem can be solved with the SBM inside the Kratos IGA Application. + +## Physical Parameters + +The material model is defined in the shared file `../FluidMaterials.json`: + +```json +{ + "name": "Newtonian3DLaw", + "Variables": { + "DENSITY": 1.0, + "DYNAMIC_VISCOSITY": 1e-3, + "PENALTY_FACTOR": 1e3 + } +} +``` + +For the immersed SBM boundary, a second property block is assigned: + +```json +{ + "name": "Newtonian3DLaw", + "Variables": { + "PENALTY_FACTOR": 0.0 + } +} +``` + +Therefore: + +- density: `rho = 1.0` +- dynamic viscosity: `mu = 1.0e-3` +- kinematic viscosity: `nu = mu / rho = 1.0e-3` + +The role of the penalty parameters is: + +- `IgaModelPart` uses `PENALTY_FACTOR = 1e3` for the fitted support conditions associated with the background box +- `IgaModelPart.SBM_Support_inner` uses `PENALTY_FACTOR = 0.0`, which activates the penalty-free skew-symmetric Nitsche-style SBM treatment on the immersed surrogate boundary [1] + +## Geometry + +### Immersed car surface + +The immersed geometry is stored in the parent Toyota folder: + +```text +../toyota_immersed_surface_sbm.mdpa +``` + +The imported surface contains: + +- `26939` nodes +- `29860` surface conditions + +From the coordinates in the `.mdpa` file, the axis-aligned extent of the immersed body is: + +```text +x in [-0.111937578, 0.074097578] -> width = 0.186035156 +y in [ 0.055194000, 0.233533669] -> height = 0.178339669 +z in [ 0.244411751, 0.699946535] -> length = 0.455534784 +``` + +In this setup, the imposed inlet velocity is aligned with the `z` direction, so `z` is the streamwise direction. + +### Background box + +The Cartesian background volume is: + +```json +"lower_point_xyz": [-0.2142325, -0.04246225, 0.01640625], +"upper_point_xyz": [ 0.1763925, 0.3481627, 1.5] +``` + +Hence the fluid box dimensions are: + +```text +Lx = 0.390625 +Ly = 0.39062495 +Lz = 1.48359375 +``` + +## Background Mesh and SBM Construction + +The background NURBS volume is created with `NurbsGeometryModelerSbm`. In this case the polynomial order is: + +```json +"polynomial_order": [1, 1, 1] +``` + +Higher-order B-splines can also be used by increasing `polynomial_order`. This usually improves approximation quality, but it also increases the number of integration points required in the elements and boundary conditions, and therefore the computational cost of the analysis. For low-Reynolds-number flows, it can make good sense to use a coarser background mesh together with a higher polynomial order, since the solution is typically smoother and the additional approximation power can be more effective than refining the mesh uniformly. + +The knot-span counts are: + +```json +"number_of_knot_spans": [25, 25, 70] +``` + +This corresponds to a structured Cartesian resolution of: + +```text +dx = 0.390625 / 25 = 0.015625 +dy = 0.39062495 / 25 ≈ 0.015625 +dz = 1.48359375 / 70 ≈ 0.0211942 +``` + +The surrogate-boundary controls are: + +- `lambda_outer = 0.5` +- `lambda_inner = 0.0` +- `number_of_inner_loops = 1` + +## Geometry Pipeline and Model Parts + +`ProjectParameters_3D_fluid_steady_state.json` follows this three-stage setup: + +1. `import_mdpa_modeler` imports `../toyota_immersed_surface_sbm.mdpa` into `initial_skin_model_part_in`. +2. `NurbsGeometryModelerSbm` creates the background Cartesian NURBS volume and classifies the immersed geometry relative to the background mesh. +3. `IgaModelerSbm` creates the analysis elements and support conditions needed by the SBM fluid formulation. + +The main model parts are: + +- `IgaModelPart.FluidDomain`: background volume where the fluid equations are solved +- `IgaModelPart.SBM_Support_inner`: immersed surrogate boundary associated with the Toyota body +- `IgaModelPart.SBM_Support_outer`: outer wall support surfaces of the box +- `IgaModelPart.Support_inlet`: inlet support surface +- `IgaModelPart.SupportOuterPressure`: outlet pressure support surface + +The immersed skin itself is also kept as a separate model part and its nodal velocity is fixed to zero. The actual fluid boundary condition, however, is imposed weakly on the surrogate boundary `IgaModelPart.SBM_Support_inner`, not by meshing the fluid to the true car surface. + +## Governing Equations and Solver + +This folder solves the Stokes variant: + +- `StokesElement` +- `solver_type = monolithic_iga` +- `time_scheme = bdf2_higher_order_vms` +- `analysis_type = non_linear` + +The nominal time data are: + +```json +"start_time": 0.0, +"end_time": 5.0, +"time_step": 1e20 +``` + +This is not a dedicated stationary solver. Instead, it reuses the monolithic IGA fluid infrastructure with a practically infinite time step, so the run behaves as a single-step Stokes solve for documentation and post-processing purposes. + +The linear solver is: + +```json +"solver_type": "bicgstab", +"preconditioner_type": "ilu0", +"tolerance": 1.0e-14, +"max_iteration": 1000, +"scaling": true +``` + +The nonlinear tolerances are: + +- relative velocity tolerance: `1e-11` +- absolute velocity tolerance: `1e-11` +- relative pressure tolerance: `1e-9` +- absolute pressure tolerance: `1e-9` + +## Boundary Conditions + +The box uses six faces. In the `IgaModelerSbm` setup: + +- `brep_ids [2]` define the inlet support +- `brep_ids [3]` define the outlet pressure support +- `brep_ids [4,5,6,7]` define the outer wall supports + +The immersed Toyota boundary is handled by: + +- `SbmFluidConditionDirichlet` on `IgaModelPart.SBM_Support_inner` + +### Outer box walls + +```text +u = (0, 0, 0) +``` + +is imposed weakly on `IgaModelPart.SBM_Support_outer` through `SupportFluidCondition`. + +### Immersed Toyota body + +```text +u = (0, 0, 0) +``` + +The true surface nodes in `skin_model_part.inner` are assigned zero velocity, and the fluid boundary condition is weakly transferred to the surrogate boundary `IgaModelPart.SBM_Support_inner`. + +### Outlet + +```text +p = 0 +NORMAL_STRESS = (0, 0, 0) +``` + +is imposed on `IgaModelPart.SupportOuterPressure`. + +### Inlet + +The Stokes case uses a constant inlet velocity: + +```text +u = (0, 0, 2.0) +``` + +## Reynolds Number + +Using: + +```text +U_ref = 2.0 +nu = 1.0e-3 +L_ref = x_max - x_min = 0.186035156 +``` + +gives the reference scale: + +```text +Re_L = U_ref * L_ref / nu + = 2.0 * 0.186035156 / 1.0e-3 + ≈ 372.07 +``` + +For this Stokes example, that Reynolds number is only a reference scale for the chosen geometry and boundary data. Inertia is not part of the equations actually solved in this folder. + +## Initial Condition + +The JSON file does not prescribe a separate nonzero initial field. In practice, the analysis starts from the default rest state. + +## Post-Processing Content + +### `plot_conditions.py` + +This helper script rebuilds the geometry/modeler setup and plots the surrogate boundary faces in 3D. It is useful to inspect which background faces have been selected by the SBM construction around the immersed vehicle and around the outer box. + +Its output file is: + +- `surrogate_faces_plot.png` + +### `run_and_post_nurbs_steady_state.py` + +The steady Stokes script writes: + +- a 3D cut figure of velocity magnitude and pressure +- a plane plot of velocity magnitude at `x ≈ 0` +- a plane plot of pressure at `x ≈ 0` + +The output files are: + +- `steady_3d_velocity_pressure.png` +- `steady_plane_x0_velocity.png` +- `steady_plane_x0_pressure.png` + +## Files + +- `ProjectParameters_3D_fluid_steady_state.json`: steady-state-style 3D Stokes setup +- `../FluidMaterials.json`: shared 3D Newtonian material data and penalty settings +- `../toyota_immersed_surface_sbm.mdpa`: shared immersed Toyota surface mesh +- `run_and_post_nurbs_steady_state.py`: steady Stokes solver driver and static post-processing +- `plot_conditions.py`: surrogate-face visualization utility + +## Python Dependencies + +The post-processing scripts use: + +- `numpy` +- `matplotlib` + +## What You Can Change + +The standard workflow for another immersed geometry is: + +1. export or generate the external surface elsewhere +2. convert it to a Kratos `.mdpa` skin +3. replace `../toyota_immersed_surface_sbm.mdpa` +4. adapt the background box so the new body is properly enclosed +5. refine `number_of_knot_spans` to obtain adequate resolution + +## Run + +Steady-state-style Stokes case: + +```bash +cd /home/nantonelli/Examples/iga/use_cases/toyota_3d_sbm/steady_stokes +python3 run_and_post_nurbs_steady_state.py +``` + +Geometry / surrogate-boundary plot: + +```bash +cd /home/nantonelli/Examples/iga/use_cases/toyota_3d_sbm/steady_stokes +python3 plot_conditions.py +``` + +## Local Refinement Requirement + +Like the 2D Turek SBM example, this case currently uses a globally structured Cartesian background mesh. This is useful for testing, debugging, and demonstrating the workflow, but it is not the most efficient strategy for serious external-aerodynamics simulations around an immersed vehicle. + +More efficient options include: + +- adaptive local refinement with THB-splines [3] +- local multipatch refinement around the vehicle and wake, coupled back to the background patch with Gap-SBM [4] + +## References + +[1] Antonelli, N., Aristio, R., Gorgi, A., Zorrilla, R., Rossi, R., Scovazzi, G., and Wüchner, R., *The Shifted Boundary Method in Isogeometric Analysis*, Computer Methods in Applied Mechanics and Engineering, Volume 430, 2024, 117228. https://doi.org/10.1016/j.cma.2024.117228 + +[2] Antonelli, N., Gorgi, A., Zorrilla, R., and Rossi, R., *Isogeometric analysis for non-Newtonian viscoplastic fluids: challenges for non-smooth solutions*, Computer Methods in Applied Mechanics and Engineering, Volume 447, 2025, 118386. https://doi.org/10.1016/j.cma.2025.118386 + +[3] Giannelli, C., Jüttler, B., and Speleers, H., *THB-splines: The truncated basis for hierarchical splines*, Computer Aided Geometric Design, Volume 29, Issue 7, 2012, pp. 485-498. https://doi.org/10.1016/j.cagd.2012.03.025 + +[4] Antonelli, N., Gorgi, A., Zorrilla, R., and Rossi, R., *Isogeometric multipatch coupling with arbitrary refinement and parametrization using the Gap-Shifted Boundary Method*, Computer Methods in Applied Mechanics and Engineering, Volume 456, 2026, 118913. https://doi.org/10.1016/j.cma.2026.118913 diff --git a/iga/use_cases/toyota_3d_sbm/plot_conditions.py b/iga/use_cases/toyota_3d_sbm/steady_stokes/plot_conditions.py similarity index 96% rename from iga/use_cases/toyota_3d_sbm/plot_conditions.py rename to iga/use_cases/toyota_3d_sbm/steady_stokes/plot_conditions.py index 6456ce99..39452423 100644 --- a/iga/use_cases/toyota_3d_sbm/plot_conditions.py +++ b/iga/use_cases/toyota_3d_sbm/steady_stokes/plot_conditions.py @@ -22,9 +22,19 @@ def _read_parameters(): script_dir = os.path.dirname(os.path.abspath(__file__)) - parameter_path = os.path.join(script_dir, "ProjectParameters_3D_fluid.json") - with open(parameter_path, "r") as parameter_file: - return KratosMultiphysics.Parameters(parameter_file.read()) + candidate_filenames = ( + "ProjectParameters_3D_fluid.json", + "ProjectParameters_3D_fluid_steady_state.json", + ) + for filename in candidate_filenames: + parameter_path = os.path.join(script_dir, filename) + if not os.path.exists(parameter_path): + continue + with open(parameter_path, "r") as parameter_file: + return KratosMultiphysics.Parameters(parameter_file.read()) + raise FileNotFoundError( + "No local ProjectParameters_3D_fluid*.json file found next to plot_conditions.py" + ) def _map_points(points): diff --git a/iga/use_cases/toyota_3d_sbm/run_and_post_nurbs_steady_state.py b/iga/use_cases/toyota_3d_sbm/steady_stokes/run_and_post_steady_state.py similarity index 98% rename from iga/use_cases/toyota_3d_sbm/run_and_post_nurbs_steady_state.py rename to iga/use_cases/toyota_3d_sbm/steady_stokes/run_and_post_steady_state.py index 8e7a8a65..0b94e9a4 100644 --- a/iga/use_cases/toyota_3d_sbm/run_and_post_nurbs_steady_state.py +++ b/iga/use_cases/toyota_3d_sbm/steady_stokes/run_and_post_steady_state.py @@ -409,7 +409,7 @@ def main(): label=r"$|v|$", ) plt.tight_layout() - plt.savefig(output_plane_velocity_path, dpi=200) + plt.savefig(output_plane_velocity_path, dpi=200, bbox_inches="tight", pad_inches=0.02) plt.close(fig_plane_velocity) plane_points_pressure, plane_values_pressure, plane_x_pressure = _extract_x_plane_layer( @@ -433,7 +433,7 @@ def main(): label=r"$p$", ) plt.tight_layout() - plt.savefig(output_plane_pressure_path, dpi=200) + plt.savefig(output_plane_pressure_path, dpi=200, bbox_inches="tight", pad_inches=0.02) plt.close(fig_plane_pressure) print(f"Saved {output_3d_path}") diff --git a/iga/use_cases/toyota_3d_sbm/ProjectParameters_3D_fluid.json b/iga/use_cases/toyota_3d_sbm/transient_navier_stokes/ProjectParameters_3D_fluid.json similarity index 74% rename from iga/use_cases/toyota_3d_sbm/ProjectParameters_3D_fluid.json rename to iga/use_cases/toyota_3d_sbm/transient_navier_stokes/ProjectParameters_3D_fluid.json index 4ff32221..6f2eb52e 100644 --- a/iga/use_cases/toyota_3d_sbm/ProjectParameters_3D_fluid.json +++ b/iga/use_cases/toyota_3d_sbm/transient_navier_stokes/ProjectParameters_3D_fluid.json @@ -7,7 +7,6 @@ "parallel_type": "OpenMP", "start_time": 0.0, "end_time": 4.0 - // "end_time": 1.0 }, "solver_settings": { @@ -17,7 +16,7 @@ "model_part_name": "IgaModelPart", "domain_size": 3, "model_import_settings": { "input_type": "use_input_model_part" }, - "material_import_settings": { "materials_filename": "FluidMaterials.json" }, + "material_import_settings": { "materials_filename": "../FluidMaterials.json" }, "linear_solver_settings": { "solver_type": "bicgstab", "tolerance": 1.0e-14, "max_iteration": 1000, "scaling": true, "preconditioner_type": "ilu0"}, // "linear_solver_settings":{ // // "solver_type": "LinearSolversApplication.sparse_lu" @@ -31,22 +30,8 @@ "absolute_velocity_tolerance": 1e-9, "relative_pressure_tolerance": 1e-7, "absolute_pressure_tolerance": 1e-7, - // "relative_velocity_tolerance": 1e-14, - // "absolute_velocity_tolerance": 1e-14, - // "relative_pressure_tolerance": 1e-12, - // "absolute_pressure_tolerance": 1e-12, "time_stepping": { - // "time_step" : 1e20, - // "time_step" : 0.0025, - // "time_step" : 1.0, - // "time_step" : 0.5, - // "time_step" : 0.2, - // "time_step" : 0.1, - // "time_step" : 0.05, // - // "time_step" : 0.02, "time_step" : 0.01, - // "time_step" : 0.005, - // "time_step" : 0.0025, "automatic_time_step": false } }, @@ -56,7 +41,7 @@ "modeler_name" : "import_mdpa_modeler", "kratos_module" : "KratosMultiphysics", "Parameters" : { - "input_filename" : "toyota_immersed_surface_sbm", + "input_filename" : "../toyota_immersed_surface_sbm", "model_part_name" : "initial_skin_model_part_in" } }, @@ -72,16 +57,15 @@ "polynomial_order" : [1, 1, 1], - // "number_of_knot_spans" : [30, - // 30, - // 90], - "number_of_knot_spans" : [25, + "number_of_knot_spans" : [25, // m1 25, 70], + // "number_of_knot_spans" : [41, // m2 + // 41, + // 120], "lambda_outer": 0.5, - "lambda_inner": 0.0, + "lambda_inner": 0.5, "number_of_inner_loops": 1, - // "skin_model_part_outer_initial_name": "initial_skin_model_part_out", "skin_model_part_inner_initial_name": "initial_skin_model_part_in", "skin_model_part_name": "skin_model_part" } @@ -153,9 +137,7 @@ "variables": [ { "variable_name": "BODY_FORCE", - "value": ["0", "0", "0"] // linear - // NAVIER STOKES - // "value": ["(cosh(x)*sinh(x)*cosh(y)*cosh(y)*cosh(z)*cosh(z) + cosh(x)*sinh(x)*sinh(y)*cosh(z)*sinh(z-y) - cosh(x)*sinh(x)*sinh(y)*cosh(y)*cosh(z)*sinh(z)) + cosh(x) - 3*cosh(x)*cosh(y)*cosh(z)", "(cosh(x)*cosh(x)*cosh(y)*cosh(z)*sinh(z-y) - sinh(x)*sinh(x)*sinh(z-y)*cosh(z-y) - sinh(x)*sinh(x)*sinh(y)*cosh(z)*cosh(z-y)) + sinh(y) - 3*sinh(x)*sinh(z-y)", "(-cosh(x)*cosh(x)*cosh(y)*sinh(y)*cosh(z)*cosh(z) - sinh(x)*sinh(x)*sinh(z-y)*cosh(y)*cosh(z) + sinh(x)*sinh(x)*sinh(y)*sinh(y)*cosh(z)*sinh(z)) + cosh(z) + 3*sinh(x)*sinh(y)*cosh(z)"] + "value": ["0", "0", "0"] } ] }, @@ -169,7 +151,6 @@ "0.0", "0.0" ] - // "value": ["cosh(x)*cosh(y)*cosh(z)","sinh(x)*sinh(z-y)","-sinh(x)*sinh(y)*cosh(z)"] }] }, { @@ -182,7 +163,6 @@ "0.0", "2.0 * 0.5*(t/2.0 + 1.0 - abs(t/2.0 - 1.0))" ] - // "value": ["cosh(x)*cosh(y)*cosh(z)","sinh(x)*sinh(z-y)","-sinh(x)*sinh(y)*cosh(z)"] }] }, { @@ -191,7 +171,6 @@ { "variable_name": "PRESSURE", "value": "0.0" - // "value": "sinh(x) + cosh(y) + sinh(z)" }, { "variable_name": "NORMAL_STRESS", @@ -218,7 +197,6 @@ "0.0", "0.0" ] - // "value": ["cosh(x)*cosh(y)*cosh(z)","sinh(x)*sinh(z-y)","-sinh(x)*sinh(y)*cosh(z)"] } } ], diff --git a/iga/use_cases/toyota_3d_sbm/transient_navier_stokes/README.md b/iga/use_cases/toyota_3d_sbm/transient_navier_stokes/README.md new file mode 100644 index 00000000..57c5a485 --- /dev/null +++ b/iga/use_cases/toyota_3d_sbm/transient_navier_stokes/README.md @@ -0,0 +1,362 @@ +# 3D Toyota SBM Transient Navier-Stokes Example + +Single-patch 3D SBM/IGA [1,2] example for incompressible transient Navier-Stokes flow around an immersed Toyota-like car body. The body is not surrounded by a body-fitted fluid mesh. Instead, the geometry is provided through the shared immersed surface `../toyota_immersed_surface_sbm.mdpa`, while the fluid is solved on a structured Cartesian background NURBS volume generated directly from the JSON settings. + +This folder contains the transient Navier-Stokes variant of the Toyota SBM example: + +- transient 3D Navier-Stokes case in `ProjectParameters_3D_fluid.json` +- post-processing driver `run_and_post_nurbs_time.py` +- geometry helper `plot_conditions.py` + +The example uses: + +- the immersed Toyota surface mesh in `.mdpa` format +- a single Cartesian background NURBS volume created with `NurbsGeometryModelerSbm` +- `IgaModelerSbm` to generate the fluid domain and the surrogate support boundaries +- SBM/IGA to impose the immersed no-slip condition without body-fitted remeshing +- weak boundary conditions on the outer box walls, inlet, outlet, and immersed body + +The aim of the example is methodological: it demonstrates how a 3D external flow problem can be solved with the SBM inside the Kratos IGA Application. + +## Physical Parameters + +The material model is defined in the shared file `../FluidMaterials.json`: + +```json +{ + "name": "Newtonian3DLaw", + "Variables": { + "DENSITY": 1.0, + "DYNAMIC_VISCOSITY": 1e-3, + "PENALTY_FACTOR": 1e3 + } +} +``` + +For the immersed SBM boundary, a second property block is assigned: + +```json +{ + "name": "Newtonian3DLaw", + "Variables": { + "PENALTY_FACTOR": 0.0 + } +} +``` + +Therefore: + +- density: `rho = 1.0` +- dynamic viscosity: `mu = 1.0e-3` +- kinematic viscosity: `nu = mu / rho = 1.0e-3` + +## Geometry + +### Immersed car surface + +The immersed geometry is stored in the parent Toyota folder: + +```text +../toyota_immersed_surface_sbm.mdpa +``` + +The imported surface contains: + +- `26939` nodes +- `29860` surface conditions + +From the coordinates in the `.mdpa` file, the axis-aligned extent of the immersed body is: + +```text +x in [-0.111937578, 0.074097578] -> width = 0.186035156 +y in [ 0.055194000, 0.233533669] -> height = 0.178339669 +z in [ 0.244411751, 0.699946535] -> length = 0.455534784 +``` + +In this setup, the imposed inlet velocity is aligned with the `z` direction, so `z` is the streamwise direction. + +### Background box + +The Cartesian background volume is: + +```json +"lower_point_xyz": [-0.2142325, -0.04246225, 0.01640625], +"upper_point_xyz": [ 0.1763925, 0.3481627, 1.5] +``` + +Hence the fluid box dimensions are: + +```text +Lx = 0.390625 +Ly = 0.39062495 +Lz = 1.48359375 +``` + +## Background Mesh and SBM Construction + +The background NURBS volume is created with `NurbsGeometryModelerSbm`. In this case the polynomial order is: + +```json +"polynomial_order": [1, 1, 1] +``` + +Higher-order B-splines can also be used by increasing `polynomial_order`. This usually improves approximation quality, but it also increases the number of integration points required in the elements and boundary conditions, and therefore the computational cost of the analysis. For low-Reynolds-number flows, it can make good sense to use a coarser background mesh together with a higher polynomial order, since the solution is typically smoother and the additional approximation power can be more effective than refining the mesh uniformly. + +The knot-span counts are: + +```json +"number_of_knot_spans": [25, 25, 70] +``` + +This corresponds to a structured Cartesian resolution of: + +```text +dx = 0.390625 / 25 = 0.015625 +dy = 0.39062495 / 25 ≈ 0.015625 +dz = 1.48359375 / 70 ≈ 0.0211942 +``` + +The surrogate-boundary controls are: + +- `lambda_outer = 0.5` +- `lambda_inner = 0.5` +- `number_of_inner_loops = 1` + +## Geometry Pipeline and Model Parts + +`ProjectParameters_3D_fluid.json` follows this three-stage setup: + +1. `import_mdpa_modeler` imports `../toyota_immersed_surface_sbm.mdpa` into `initial_skin_model_part_in`. +2. `NurbsGeometryModelerSbm` creates the background Cartesian NURBS volume and classifies the immersed geometry relative to the background mesh. +3. `IgaModelerSbm` creates the analysis elements and support conditions needed by the SBM fluid formulation. + +The main model parts are: + +- `IgaModelPart.FluidDomain`: background volume where the fluid equations are solved +- `IgaModelPart.SBM_Support_inner`: immersed surrogate boundary associated with the Toyota body +- `IgaModelPart.SBM_Support_outer`: outer wall support surfaces of the box +- `IgaModelPart.Support_inlet`: inlet support surface +- `IgaModelPart.SupportOuterPressure`: outlet pressure support surface + +The immersed skin itself is also kept as a separate model part and its nodal velocity is fixed to zero. The actual fluid boundary condition, however, is imposed weakly on the surrogate boundary `IgaModelPart.SBM_Support_inner`, not by meshing the fluid to the true car surface. + +## Governing Equations and Solver + +This folder solves the transient Navier-Stokes variant: + +- `NavierStokesElement` +- monolithic velocity-pressure solution +- `time_scheme = bdf2_higher_order_vms` +- `analysis_type = non_linear` + +The time data are: + +```json +"start_time": 0.0, +"end_time": 4.0, +"time_step": 0.01 +``` + +So the run advances over `400` nominal time steps. + +The linear solver is: + +```json +"solver_type": "bicgstab", +"preconditioner_type": "ilu0", +"tolerance": 1.0e-14, +"max_iteration": 1000, +"scaling": true +``` + +The nonlinear tolerances are: + +- relative velocity tolerance: `1e-9` +- absolute velocity tolerance: `1e-9` +- relative pressure tolerance: `1e-7` +- absolute pressure tolerance: `1e-7` + +## Boundary Conditions + +The box uses six faces. In the `IgaModelerSbm` setup: + +- `brep_ids [2]` define the inlet support +- `brep_ids [3]` define the outlet pressure support +- `brep_ids [4,5,6,7]` define the outer wall supports + +The immersed Toyota boundary is handled by: + +- `SbmFluidConditionDirichlet` on `IgaModelPart.SBM_Support_inner` + +### Outer box walls + +```text +u = (0, 0, 0) +``` + +is imposed weakly on `IgaModelPart.SBM_Support_outer` through `SupportFluidCondition`. + +### Immersed Toyota body + +```text +u = (0, 0, 0) +``` + +The true surface nodes in `skin_model_part.inner` are assigned zero velocity, and the fluid boundary condition is weakly transferred to the surrogate boundary `IgaModelPart.SBM_Support_inner`. + +### Outlet + +```text +p = 0 +NORMAL_STRESS = (0, 0, 0) +``` + +is imposed on `IgaModelPart.SupportOuterPressure`. + +### Inlet + +The transient case prescribes a streamwise inlet velocity in the `z` direction: + +```json +"VELOCITY": [ + "0.0", + "0.0", + "2.0 * 0.5*(t/2.0 + 1.0 - abs(t/2.0 - 1.0))" +] +``` + +This can be written in piecewise form as: + +```text +u_x = 0 +u_y = 0 +u_z = t for 0 <= t <= 2 +u_z = 2.0 for t >= 2 +``` + +So the inlet is ramped linearly from rest to `2.0` during the first `2` time units and then kept constant. + +## Reynolds Number + +Using: + +```text +U_ref = 2.0 +nu = 1.0e-3 +L_ref = x_max - x_min = 0.186035156 +``` + +gives: + +```text +Re_L = U_ref * L_ref / nu + = 2.0 * 0.186035156 / 1.0e-3 + ≈ 372.07 +``` + +This should be understood as a reference Reynolds number based on the `x` extent of the immersed geometry. + +## Initial Condition + +The JSON file does not prescribe a separate nonzero initial field. In practice, the analysis starts from the default rest state and the inlet ramp reduces the startup shock. + +## Post-Processing Content + +### `plot_conditions.py` + +This helper script rebuilds the geometry/modeler setup and plots the surrogate boundary faces in 3D. It is useful to inspect which background faces have been selected by the SBM construction around the immersed vehicle and around the outer box. + +Its output file is: + +- `surrogate_faces_plot.png` + +### `run_and_post_nurbs_time.py` + +The transient script samples element-centered average velocity and pressure and generates: + +- a 3D cut visualization for elements with `x < 0`, colored with `|v_z|` +- a contour plot of velocity magnitude on the plane closest to `x = 0` +- a contour plot of pressure on the plane closest to `x = 0` +- animated GIFs from the stored frame sequences + +The output files are: + +- `cut_contour_x_lt_0.gif` +- `velocity_magnitude_x_eq_0_contour.gif` +- `pressure_x_eq_0_contour.gif` + +The transient frame caches are written in: + +- `frames_time_cut_contour` +- `frames_time_velocity_magnitude_x0` +- `frames_time_pressure_x0` + +## Files + +- `ProjectParameters_3D_fluid.json`: transient 3D Navier-Stokes setup +- `../FluidMaterials.json`: shared 3D Newtonian material data and penalty settings +- `../toyota_immersed_surface_sbm.mdpa`: shared immersed Toyota surface mesh +- `run_and_post_nurbs_time.py`: transient solver driver and GIF post-processing +- `plot_conditions.py`: surrogate-face visualization utility + +## Python Dependencies + +The post-processing scripts use: + +- `numpy` +- `matplotlib` +- `imageio` + +## What You Can Change + +The standard workflow for another immersed geometry is: + +1. export or generate the external surface elsewhere +2. convert it to a Kratos `.mdpa` skin +3. replace `../toyota_immersed_surface_sbm.mdpa` +4. adapt the background box so the new body is properly enclosed +5. refine `number_of_knot_spans` to obtain adequate wake resolution + +## Run + +Transient Navier-Stokes case: + +```bash +cd /home/nantonelli/Examples/iga/use_cases/toyota_3d_sbm/transient_navier_stokes +python3 run_and_post_nurbs_time.py +``` + +Geometry / surrogate-boundary plot: + +```bash +cd /home/nantonelli/Examples/iga/use_cases/toyota_3d_sbm/transient_navier_stokes +python3 plot_conditions.py +``` + +## Local Refinement Requirement + +Like the 2D Turek SBM example, this case currently uses a globally structured Cartesian background mesh. This is useful for testing, debugging, and demonstrating the workflow, but it is not the most efficient strategy for serious external-aerodynamics simulations around an immersed vehicle. + +For this geometry, the flow features of greatest interest are typically concentrated: + +- near the front stagnation region +- around the roof and side edges +- in the near wake behind the vehicle +- in any separated shear layers generated by the immersed geometry + +Using global refinement to resolve these features increases the total number of degrees of freedom everywhere in the box, including regions where the solution is comparatively smooth. + +More efficient options include: + +- adaptive local refinement with THB-splines [3] +- local multipatch refinement around the vehicle and wake, coupled back to the background patch with Gap-SBM [4] + +## References + +[1] Antonelli, N., Aristio, R., Gorgi, A., Zorrilla, R., Rossi, R., Scovazzi, G., and Wüchner, R., *The Shifted Boundary Method in Isogeometric Analysis*, Computer Methods in Applied Mechanics and Engineering, Volume 430, 2024, 117228. https://doi.org/10.1016/j.cma.2024.117228 + +[2] Antonelli, N., Gorgi, A., Zorrilla, R., and Rossi, R., *Isogeometric analysis for non-Newtonian viscoplastic fluids: challenges for non-smooth solutions*, Computer Methods in Applied Mechanics and Engineering, Volume 447, 2025, 118386. https://doi.org/10.1016/j.cma.2025.118386 + +[3] Giannelli, C., Jüttler, B., and Speleers, H., *THB-splines: The truncated basis for hierarchical splines*, Computer Aided Geometric Design, Volume 29, Issue 7, 2012, pp. 485-498. https://doi.org/10.1016/j.cagd.2012.03.025 + +[4] Antonelli, N., Gorgi, A., Zorrilla, R., and Rossi, R., *Isogeometric multipatch coupling with arbitrary refinement and parametrization using the Gap-Shifted Boundary Method*, Computer Methods in Applied Mechanics and Engineering, Volume 456, 2026, 118913. https://doi.org/10.1016/j.cma.2026.118913 diff --git a/iga/use_cases/toyota_3d_sbm/transient_navier_stokes/plot_conditions.py b/iga/use_cases/toyota_3d_sbm/transient_navier_stokes/plot_conditions.py new file mode 100644 index 00000000..39452423 --- /dev/null +++ b/iga/use_cases/toyota_3d_sbm/transient_navier_stokes/plot_conditions.py @@ -0,0 +1,477 @@ +import importlib +import os + +import matplotlib.pyplot as plt +import numpy as np +import KratosMultiphysics +from matplotlib.patches import Patch +from matplotlib.ticker import MaxNLocator +from mpl_toolkits.mplot3d.art3d import Poly3DCollection + +try: + import KratosMultiphysics.IgaApplication as IGA # noqa: F401 + IGA_IMPORT_ERROR = None +except ImportError as exc: + IGA = None + IGA_IMPORT_ERROR = exc + + +PLOT_SWAP_YZ_FOR_VISUALIZATION = True +PLOT_PROJECTIONS = False + + +def _read_parameters(): + script_dir = os.path.dirname(os.path.abspath(__file__)) + candidate_filenames = ( + "ProjectParameters_3D_fluid.json", + "ProjectParameters_3D_fluid_steady_state.json", + ) + for filename in candidate_filenames: + parameter_path = os.path.join(script_dir, filename) + if not os.path.exists(parameter_path): + continue + with open(parameter_path, "r") as parameter_file: + return KratosMultiphysics.Parameters(parameter_file.read()) + raise FileNotFoundError( + "No local ProjectParameters_3D_fluid*.json file found next to plot_conditions.py" + ) + + +def _map_points(points): + if not PLOT_SWAP_YZ_FOR_VISUALIZATION or points.size == 0: + return points + return points[:, [0, 2, 1]] + + +def _map_bounds(bounds): + if not PLOT_SWAP_YZ_FOR_VISUALIZATION: + return bounds + x_min, x_max, y_min, y_max, z_min, z_max = bounds + return (x_min, x_max, z_min, z_max, y_min, y_max) + + +def _axis_labels(): + if not PLOT_SWAP_YZ_FOR_VISUALIZATION: + return ("X", "Y", "Z") + return ("X", "Z", "Y") + + +def _get_plot_bounds(parameters): + for modeler in parameters["modelers"].values(): + if not modeler.Has("modeler_name"): + continue + if modeler["modeler_name"].GetString() != "NurbsGeometryModelerSbm": + continue + + modeler_params = modeler["Parameters"] + if modeler_params.Has("lower_point_xyz") and modeler_params.Has("upper_point_xyz"): + lower = modeler_params["lower_point_xyz"] + upper = modeler_params["upper_point_xyz"] + elif modeler_params.Has("lower_point_uvw") and modeler_params.Has("upper_point_uvw"): + lower = modeler_params["lower_point_uvw"] + upper = modeler_params["upper_point_uvw"] + else: + continue + + return ( + lower[0].GetDouble(), upper[0].GetDouble(), + lower[1].GetDouble(), upper[1].GetDouble(), + lower[2].GetDouble(), upper[2].GetDouble(), + ) + return None + + +def _collect_surrogate_faces(model_part): + faces = [] + + # The plottable surrogate faces in this setup are stored directly on IgaModelPart + # as Quadrilateral3D4 conditions. The SBM_Support_* submodelparts contain + # quadrature-point wrapper geometries, which are not the face polygons we want here. + for cond in model_part.Conditions: + geom = cond.GetGeometry() + if geom.PointsNumber() < 3 or geom.PointsNumber() > 4: + continue + faces.append(np.array([[node.X, node.Y, node.Z] for node in geom], dtype=float)) + + return faces + + +def _resolve_variable(variable_name): + try: + if KratosMultiphysics.KratosGlobals.HasVariable(variable_name): + return KratosMultiphysics.KratosGlobals.GetVariable(variable_name) + except Exception: + return None + return None + + +def _as_point_array(value): + try: + return np.array([float(value[0]), float(value[1]), float(value[2])], dtype=float) + except Exception: + try: + return np.array([float(value.X), float(value.Y), float(value.Z)], dtype=float) + except Exception: + return None + + +def _read_projection_from_entity(entity, projection_node_variable, projection_coordinates_variable): + if projection_node_variable is not None and entity.Has(projection_node_variable): + try: + projection = _as_point_array(entity.GetValue(projection_node_variable)) + if projection is not None: + return projection, "PROJECTION_NODE" + except Exception: + pass + + if projection_coordinates_variable is not None and entity.Has(projection_coordinates_variable): + try: + projection = _as_point_array(entity.GetValue(projection_coordinates_variable)) + if projection is not None: + return projection, "PROJECTION_NODE_COORDINATES" + except Exception: + pass + + return None, None + + +def _collect_inner_skin_model_parts(model): + candidate_names = [] + preferred_names = ( + "skin_model_part.inner", + "SkinModelPart.inner", + "IgaModelPart.skin_model_part.inner", + ) + + for model_part_name in preferred_names: + if model.HasModelPart(model_part_name): + candidate_names.append(model_part_name) + + model_part_names = model.GetModelPartNames() if hasattr(model, "GetModelPartNames") else [] + for model_part_name in model_part_names: + if "skin_model_part.inner" not in model_part_name.lower(): + continue + if model_part_name in candidate_names: + continue + candidate_names.append(model_part_name) + + model_parts = [] + for model_part_name in candidate_names: + try: + model_parts.append(model[model_part_name]) + except Exception: + continue + return model_parts + + +def _collect_inner_skin_points(model): + points = [] + seen = set() + + for skin_model_part in _collect_inner_skin_model_parts(model): + for node in skin_model_part.Nodes: + point = (float(node.X), float(node.Y), float(node.Z)) + key = tuple(round(value, 12) for value in point) + if key in seen: + continue + seen.add(key) + points.append(point) + + if not points: + return np.empty((0, 3), dtype=float) + + return np.array(points, dtype=float) + + +def _read_projection_from_node_id(entity, projection_node_id_variable, skin_model_parts): + if projection_node_id_variable is None or not entity.Has(projection_node_id_variable): + return None + + try: + projection_node_id = int(entity.GetValue(projection_node_id_variable)) + except Exception: + return None + + if projection_node_id <= 0: + return None + + for skin_model_part in skin_model_parts: + if not skin_model_part.HasNode(projection_node_id): + continue + node = skin_model_part.GetNode(projection_node_id) + return np.array([float(node.X), float(node.Y), float(node.Z)], dtype=float) + + return None + + +def _closest_skin_node_projection(source, skin_points): + if len(skin_points) == 0: + return None + + deltas = skin_points - source + distances = np.einsum("ij,ij->i", deltas, deltas) + return np.array(skin_points[int(np.argmin(distances))], dtype=float) + + +def _execute_modelers(parameters, model): + if not model.HasModelPart("IgaModelPart"): + model.CreateModelPart("IgaModelPart") + + for modeler_data in parameters["modelers"].values(): + if not modeler_data.Has("modeler_name"): + continue + + modeler_name = modeler_data["modeler_name"].GetString() + modeler_params = modeler_data["Parameters"] + print(f"[DEBUG] Creating modeler '{modeler_name}'") + + if KratosMultiphysics.HasModeler(modeler_name): + modeler = KratosMultiphysics.CreateModeler(modeler_name, model, modeler_params) + else: + if not modeler_data.Has("kratos_module"): + raise RuntimeError(f"Python modeler '{modeler_name}' is missing 'kratos_module'.") + + kratos_module_name = modeler_data["kratos_module"].GetString() + if not kratos_module_name.startswith("KratosMultiphysics"): + kratos_module_name = "KratosMultiphysics." + kratos_module_name + python_modeler_name = kratos_module_name + ".modelers." + modeler_name + python_module = importlib.import_module(python_modeler_name) + modeler = python_module.Factory(model, modeler_params) + + setup_geometry = getattr(modeler, "SetupGeometryModel", None) + if callable(setup_geometry): + print(f"[DEBUG] - SetupGeometryModel for '{modeler_name}'") + setup_geometry() + + setup_model_part = getattr(modeler, "SetupModelPart", None) + if callable(setup_model_part): + print(f"[DEBUG] - SetupModelPart for '{modeler_name}'") + setup_model_part() + + +def _style_axis(ax, title, bounds): + x_min, x_max, y_min, y_max, z_min, z_max = bounds + axis_labels = _axis_labels() + ax.set_title(title, pad=10) + ax.set_xlim(x_min, x_max) + ax.set_ylim(y_min, y_max) + ax.set_zlim(z_min, z_max) + ax.set_box_aspect((x_max - x_min, y_max - y_min, z_max - z_min)) + ax.view_init(elev=20, azim=-60) + ax.grid(False) + for axis in (ax.xaxis, ax.yaxis, ax.zaxis): + axis.pane.set_edgecolor((0, 0, 0, 0.08)) + axis.pane.set_facecolor((0.98, 0.98, 0.98, 1.0)) + ax.xaxis.set_major_locator(MaxNLocator(nbins=5)) + ax.set_xlabel(axis_labels[0]) + ax.set_ylabel(axis_labels[1]) + ax.set_zlabel(axis_labels[2]) + + +def _collect_inner_projection_segments(model_part, model): + if not PLOT_PROJECTIONS: + return [] + + projection_node_variable = _resolve_variable("PROJECTION_NODE") + projection_coordinates_variable = _resolve_variable("PROJECTION_NODE_COORDINATES") + projection_node_id_variable = _resolve_variable("PROJECTION_NODE_ID") + skin_model_parts = _collect_inner_skin_model_parts(model) + skin_points = _collect_inner_skin_points(model) + + if not model_part.HasSubModelPart("SBM_Support_inner"): + return [] + + segments = [] + seen = set() + condition_projection_node_count = 0 + condition_projection_coordinates_count = 0 + gauss_projection_node_count = 0 + gauss_projection_coordinates_count = 0 + gauss_projection_node_id_count = 0 + gauss_closest_skin_node_count = 0 + gauss_projection_cache = {} + + def _try_add_segment(source, projection): + key = tuple(round(value, 12) for value in (*source, *projection)) + if key in seen: + return False + seen.add(key) + segments.append((source, projection)) + return True + + for condition in model_part.GetSubModelPart("SBM_Support_inner").Conditions: + center = condition.GetGeometry().Center() + source = np.array([float(center.X), float(center.Y), float(center.Z)], dtype=float) + projection, source_label = _read_projection_from_entity( + condition, + projection_node_variable, + projection_coordinates_variable, + ) + if projection is not None and _try_add_segment(source, projection): + if source_label == "PROJECTION_NODE": + condition_projection_node_count += 1 + elif source_label == "PROJECTION_NODE_COORDINATES": + condition_projection_coordinates_count += 1 + + for gauss_point in condition.GetGeometry(): + gauss_source = np.array([float(gauss_point.X), float(gauss_point.Y), float(gauss_point.Z)], dtype=float) + projection, source_label = _read_projection_from_entity( + gauss_point, + projection_node_variable, + projection_coordinates_variable, + ) + if projection is None: + projection = _read_projection_from_node_id( + gauss_point, + projection_node_id_variable, + skin_model_parts, + ) + source_label = "PROJECTION_NODE_ID" if projection is not None else None + if projection is None: + cache_key = gauss_point.Id + if cache_key not in gauss_projection_cache: + gauss_projection_cache[cache_key] = _closest_skin_node_projection(gauss_source, skin_points) + projection = gauss_projection_cache[cache_key] + source_label = "CLOSEST_SKIN_NODE" if projection is not None else None + + if projection is None or not _try_add_segment(gauss_source, projection): + continue + + if source_label == "PROJECTION_NODE": + gauss_projection_node_count += 1 + elif source_label == "PROJECTION_NODE_COORDINATES": + gauss_projection_coordinates_count += 1 + elif source_label == "PROJECTION_NODE_ID": + gauss_projection_node_id_count += 1 + elif source_label == "CLOSEST_SKIN_NODE": + gauss_closest_skin_node_count += 1 + + print( + "SBM_Support_inner projections: " + f"{condition_projection_node_count} from condition PROJECTION_NODE, " + f"{condition_projection_coordinates_count} from condition PROJECTION_NODE_COORDINATES, " + f"{gauss_projection_node_count} from gauss-point PROJECTION_NODE, " + f"{gauss_projection_coordinates_count} from gauss-point PROJECTION_NODE_COORDINATES, " + f"{gauss_projection_node_id_count} from gauss-point PROJECTION_NODE_ID, " + f"{gauss_closest_skin_node_count} from closest inner skin nodes, " + f"{len(segments)} total" + ) + return segments + + +def _plot_surrogate_faces(ax, faces, bounds): + axis_labels = _axis_labels() + face_groups = { + 0: {"color": (1.0, 0.7, 0.7), "label": f"Faces normal to {axis_labels[0]}"}, + 1: {"color": (0.7, 1.0, 0.7), "label": f"Faces normal to {axis_labels[1]}"}, + 2: {"color": (0.7, 0.7, 1.0), "label": f"Faces normal to {axis_labels[2]}"}, + } + face_count = 0 + + for face_pts in faces: + plot_face_pts = _map_points(face_pts) + normal = np.cross(plot_face_pts[1] - plot_face_pts[0], plot_face_pts[2] - plot_face_pts[0]) + normal_norm = np.linalg.norm(normal) + if normal_norm <= 1e-12: + continue + + normal /= normal_norm + group_id = int(np.argmax(np.abs(normal))) + poly = Poly3DCollection( + [plot_face_pts], + alpha=0.75, + edgecolor="k", + linewidth=0.3, + facecolor=face_groups[group_id]["color"], + ) + ax.add_collection3d(poly) + face_count += 1 + + _style_axis(ax, "Surrogate Boundary Faces", _map_bounds(bounds)) + ax.legend( + handles=[ + Patch(facecolor=group["color"], edgecolor="k", label=group["label"]) + for group in face_groups.values() + ], + loc="upper right", + ) + return face_count + + +def _plot_inner_projection_segments(ax, segments): + if not PLOT_PROJECTIONS: + return + + if not segments: + print("Warning: no inner projection segments found") + return + + for source, projection in segments: + plot_segment = _map_points(np.array([source, projection], dtype=float)) + ax.plot( + plot_segment[:, 0], + plot_segment[:, 1], + plot_segment[:, 2], + color="#111111", + linewidth=0.7, + alpha=0.7, + ) + ax.scatter( + [plot_segment[0, 0]], + [plot_segment[0, 1]], + [plot_segment[0, 2]], + s=6, + color="#111111", + alpha=0.8, + depthshade=False, + ) + + +def main(): + if IGA_IMPORT_ERROR is not None: + raise RuntimeError( + "KratosMultiphysics.IgaApplication failed to import. " + "The current build appears to expose PROJECTION_NODE incorrectly in Python, " + "so the IGA modelers/conditions cannot be created from this interpreter." + ) from IGA_IMPORT_ERROR + + parameters = _read_parameters() + plot_bounds = _get_plot_bounds(parameters) + + global_model = KratosMultiphysics.Model() + _execute_modelers(parameters, global_model) + + iga_model_part = global_model["IgaModelPart"] + surrogate_faces = _collect_surrogate_faces(iga_model_part) + projection_segments = _collect_inner_projection_segments(iga_model_part, global_model) + if not surrogate_faces: + raise RuntimeError("No surrogate face conditions were found on IgaModelPart.") + + if plot_bounds is None: + all_points = [] + for face in surrogate_faces: + for point in face: + all_points.append(tuple(point)) + if not all_points: + raise RuntimeError("No surrogate faces found to plot.") + xs, ys, zs = zip(*all_points) + plot_bounds = (min(xs), max(xs), min(ys), max(ys), min(zs), max(zs)) + + fig = plt.figure(figsize=(10, 8)) + ax = fig.add_subplot(1, 1, 1, projection="3d") + face_count = _plot_surrogate_faces(ax, surrogate_faces, plot_bounds) + if face_count == 0: + plt.close(fig) + raise RuntimeError("No surrogate conditions with 3 or more nodes were found.") + _plot_inner_projection_segments(ax, projection_segments) + + plt.tight_layout() + script_dir = os.path.dirname(os.path.abspath(__file__)) + output_path = os.path.join(script_dir, "surrogate_faces_plot.png") + plt.savefig(output_path, dpi=200) + print(f"Saved surrogate face plot to {output_path}") + plt.show() + + +if __name__ == "__main__": + main() diff --git a/iga/use_cases/toyota_3d_sbm/run_and_post_nurbs_time.py b/iga/use_cases/toyota_3d_sbm/transient_navier_stokes/run_and_post_transient.py similarity index 57% rename from iga/use_cases/toyota_3d_sbm/run_and_post_nurbs_time.py rename to iga/use_cases/toyota_3d_sbm/transient_navier_stokes/run_and_post_transient.py index 8a7648bb..037f6089 100644 --- a/iga/use_cases/toyota_3d_sbm/run_and_post_nurbs_time.py +++ b/iga/use_cases/toyota_3d_sbm/transient_navier_stokes/run_and_post_transient.py @@ -1,11 +1,5 @@ #!/usr/bin/env python3 -"""Transient Toyota post-processing example. - -This script runs the transient SBM/IGA case and then post-processes the -final solution. The particle GIF is generated on the frozen last velocity -field, which is much faster and clearer than plotting particles at every -time step. -""" +"""Transient Toyota post-processing example.""" import importlib import os @@ -25,7 +19,8 @@ import matplotlib.pyplot as plt from matplotlib.colors import Normalize from matplotlib.ticker import MaxNLocator, StrMethodFormatter - from mpl_toolkits.mplot3d.art3d import Line3DCollection, Poly3DCollection + from mpl_toolkits.mplot3d.art3d import Poly3DCollection + from mpl_toolkits.axes_grid1 import make_axes_locatable HAVE_MPL = True @@ -60,7 +55,6 @@ def _font_available(font_name): HAVE_MPL = False plt = None Normalize = None - Line3DCollection = None Poly3DCollection = None try: @@ -73,6 +67,12 @@ def _font_available(font_name): PLOT_SWAP_Z_TO_X = True PLOT_MODEL_PART_NAME = "IgaModelPart" +PROJECT_PARAMETERS_FILENAME = "ProjectParameters_3D_fluid.json" + +WRITE_3D_CUT_GIF = True +WRITE_VELOCITY_PLANE_GIF = True +WRITE_PRESSURE_PLANE_GIF = True + CAMERA_ELEV = 20.0 CAMERA_AZIM = 322.0 COLORBAR_MAX = 4.5 @@ -81,17 +81,23 @@ def _font_available(font_name): GIF_DURATION = 0.05 GIF_PALETTE_SIZE = 64 -PARTICLE_SEED_NX = 10 -PARTICLE_SEED_NY = 5 -PARTICLE_SEED_Z_OFFSET = 0.02 -PARTICLE_STEP_FRACTION = 0.02 -PARTICLE_ALPHA = 0.9 -TRAIL_ALPHA = 0.7 -TRAIL_LINEWIDTH = 1.0 -USE_PHYSICAL_DT = True -PARTICLE_POST_STEPS = 180 -PARTICLE_POST_WAVES = 4 +CUT_FRAME_DIRNAME = "frames_time_cut_contour" +VELOCITY_PLANE_FRAME_DIRNAME = "frames_time_velocity_magnitude_x0" +PRESSURE_PLANE_FRAME_DIRNAME = "frames_time_pressure_x0" +CUT_GIF_FILENAME = "cut_contour_x_lt_0.gif" +VELOCITY_PLANE_GIF_FILENAME = "velocity_magnitude_x_eq_0_contour.gif" +PRESSURE_PLANE_GIF_FILENAME = "pressure_x_eq_0_contour.gif" + +PLANE_FIGSIZE = (8.0, 6.0) +PLANE_TITLE_SIZE = 14 +PLANE_LABEL_SIZE = 13 +PLANE_TICK_SIZE = 11 +PLANE_COLORBAR_LABEL_SIZE = 13 +PLANE_COLORBAR_TICK_SIZE = 11 +PLANE_COLORBAR_SIZE = "4%" +PLANE_COLORBAR_PAD = 0.08 +PLANE_COLORBAR_NBINS = 6 def _map_points(points): if not PLOT_SWAP_Z_TO_X or points.size == 0: @@ -99,10 +105,6 @@ def _map_points(points): return points[:, [0, 2, 1]] -def _axis_labels(): - return ("X", "Z", "Y") if PLOT_SWAP_Z_TO_X else ("X", "Y", "Z") - - def _normalize_solver_settings(parameters): if not parameters.Has("solver_settings"): return @@ -228,31 +230,63 @@ def _spacing(values): return marker_size_pts ** 2 -def _plot_x_plane(ax, plane_points, plane_values, color_norm): +def _plot_x_plane(ax, plane_points, plane_values, color_norm, cmap="jet"): y = plane_points[:, 1] z = plane_points[:, 2] marker_size = _square_marker_size(ax, z, y) + artist = ax.scatter( z, y, c=plane_values, - cmap="jet", + cmap=cmap, norm=color_norm, s=marker_size, marker="s", linewidths=0.0, + rasterized=True, ) + ax.set_facecolor("white") - ax.set_xlabel(r"$Z$") - ax.set_ylabel(r"$Y$") + ax.set_xlabel(r"$Z$", fontsize=PLANE_LABEL_SIZE) + ax.set_ylabel(r"$Y$", fontsize=PLANE_LABEL_SIZE) + + if z.size: + ax.set_xlim(float(z.min()), float(z.max())) + if y.size: + ax.set_ylim(float(y.min()), float(y.max())) + + ax.margins(x=0.0, y=0.0) ax.xaxis.set_major_locator(MaxNLocator(nbins=5)) ax.xaxis.set_major_formatter(StrMethodFormatter("{x:.2f}")) ax.yaxis.set_major_locator(MaxNLocator(nbins=5)) ax.yaxis.set_major_formatter(StrMethodFormatter("{x:.2f}")) - ax.set_aspect("equal") + + ax.tick_params(axis="both", which="major", labelsize=PLANE_TICK_SIZE) + ax.set_aspect("equal", adjustable="box") + return artist +def _add_2d_colorbar(fig, ax, artist, label): + divider = make_axes_locatable(ax) + cax = divider.append_axes( + "right", + size=PLANE_COLORBAR_SIZE, + pad=PLANE_COLORBAR_PAD, + ) + + colorbar = fig.colorbar( + artist, + cax=cax, + ) + + colorbar.set_label(label, fontsize=PLANE_COLORBAR_LABEL_SIZE) + colorbar.ax.tick_params(labelsize=PLANE_COLORBAR_TICK_SIZE) + colorbar.locator = MaxNLocator(nbins=PLANE_COLORBAR_NBINS) + colorbar.update_ticks() + + return colorbar def _collect_surrogate_faces(model): names = [] @@ -309,7 +343,7 @@ def _add_surrogate_faces(ax, faces, alpha): ax.add_collection3d(poly) -def _set_3d_axes(ax, bounds, axis_labels): +def _set_3d_axes(ax, bounds): x_min, x_max, y_min, y_max, z_min, z_max = bounds x_span = x_max - x_min y_span = y_max - y_min @@ -326,11 +360,10 @@ def _set_3d_axes(ax, bounds, axis_labels): ax.set_zlim(z_min, z_max) ax.set_box_aspect((x_span, y_span, z_span)) - ax.set_xlabel(axis_labels[0]) - ax.set_ylabel(axis_labels[1]) - ax.set_zlabel(axis_labels[2]) - ax.xaxis.set_major_locator(MaxNLocator(nbins=3)) - ax.xaxis.set_major_formatter(StrMethodFormatter("{x:.2f}")) + ax.set_xlabel("") + ax.set_ylabel("") + ax.set_zlabel("") + ax.set_xticks([]) ax.view_init(elev=CAMERA_ELEV, azim=CAMERA_AZIM) @@ -366,173 +399,8 @@ def _write_gif(frame_dir, gif_path): print(f"Saved GIF: {gif_path}") -def _seed_particle_points(bounds): - x_min, x_max, y_min, y_max, z_min, z_max = bounds - x_span = x_max - x_min - y_span = y_max - y_min - z_span = z_max - z_min - - seed_xs = np.linspace(x_min + 0.01 * x_span, x_max - 0.01 * x_span, PARTICLE_SEED_NX) - seed_ys = np.linspace(y_min + 0.01 * y_span, y_max - 0.01 * y_span, PARTICLE_SEED_NY) - seed_z = z_min + PARTICLE_SEED_Z_OFFSET * z_span - return [ - np.array([x, y, seed_z], dtype=float) - for x in seed_xs - for y in seed_ys - ] - - -def _nearest_velocity(points, velocity, point): - diff = points - point - nearest = np.argmin(np.einsum("ij,ij->i", diff, diff)) - return velocity[nearest] - - -def _render_particle_frame( - frame_path, - bounds, - axis_labels, - surrogate_faces, - velocity_norm, - step_index, - particle_positions, - particle_paths, - particle_speeds, -): - fig = plt.figure(figsize=(9, 6)) - ax = fig.add_subplot(1, 1, 1, projection="3d") - _set_3d_axes(ax, bounds, axis_labels) - ax.set_title(f"Particle post-process on final field (step {step_index})") - - if surrogate_faces: - faces = [_map_points(np.array(face, dtype=float)).tolist() for face in surrogate_faces] - _add_surrogate_faces(ax, faces, alpha=0.35) - - segments = [] - segment_speeds = [] - for path, speed_history in zip(particle_paths, particle_speeds): - if len(path) < 2: - continue - for index in range(1, len(path)): - segments.append(np.array([path[index - 1], path[index]], dtype=float)) - segment_speeds.append( - speed_history[index - 1] if index - 1 < len(speed_history) else 0.0 - ) - - mappable = None - if segments: - mapped_segments = [_map_points(segment) for segment in segments] - trails = Line3DCollection( - mapped_segments, - cmap="jet", - norm=velocity_norm, - linewidth=TRAIL_LINEWIDTH, - alpha=TRAIL_ALPHA, - ) - trails.set_array(np.array(segment_speeds, dtype=float)) - ax.add_collection3d(trails) - mappable = trails - - if particle_positions: - point_speeds = np.array([history[-1] if history else 0.0 for history in particle_speeds], dtype=float) - plot_positions = _map_points(np.array(particle_positions, dtype=float)) - particles = ax.scatter( - plot_positions[:, 0], - plot_positions[:, 1], - plot_positions[:, 2], - c=point_speeds, - cmap="jet", - norm=velocity_norm, - s=10, - alpha=PARTICLE_ALPHA, - ) - if mappable is None: - mappable = particles - - if mappable is not None: - plt.colorbar(mappable, ax=ax, shrink=0.8, pad=0.03, label=r"$|v_z|$") - - plt.tight_layout() - plt.savefig(frame_path, dpi=FIG_DPI) - plt.close(fig) - - -def _create_particle_postprocess_gif( - frame_dir, - bounds, - axis_labels, - surrogate_faces, - velocity_norm, - points, - velocity, - physical_dt, -): - seed_points = _seed_particle_points(bounds) - particle_positions = [point.copy() for point in seed_points] - particle_paths = [[point.copy()] for point in seed_points] - particle_speeds = [[] for _ in seed_points] - - x_min, x_max, y_min, y_max, z_min, z_max = bounds - x_span = x_max - x_min - y_span = y_max - y_min - z_span = z_max - z_min - step_length = PARTICLE_STEP_FRACTION * min(x_span, y_span, z_span) - field_speed_max = float(np.max(np.abs(velocity[:, 2]))) if velocity.size > 0 else 0.0 - wave_interval = max(1, PARTICLE_POST_STEPS // PARTICLE_POST_WAVES) - - for step_index in range(PARTICLE_POST_STEPS): - if step_index > 0 and step_index % wave_interval == 0: - particle_positions.extend(point.copy() for point in seed_points) - particle_paths.extend([point.copy()] for point in seed_points) - particle_speeds.extend([] for _ in seed_points) - - new_positions = [] - new_paths = [] - new_speeds = [] - for point, path, speed_history in zip(particle_positions, particle_paths, particle_speeds): - local_velocity = _nearest_velocity(points, velocity, point) - speed = float(abs(local_velocity[2])) - if speed > 1e-12: - if physical_dt is not None: - next_point = point + local_velocity * physical_dt - else: - direction = local_velocity / speed - step = step_length - if field_speed_max > 0.0: - step *= speed / field_speed_max - next_point = point + direction * step - else: - next_point = point - - inside_domain = ( - x_min <= next_point[0] <= x_max - and y_min <= next_point[1] <= y_max - and z_min <= next_point[2] <= z_max - ) - if not inside_domain: - continue - - next_point = np.array(next_point, dtype=float) - new_positions.append(next_point) - new_paths.append(path + [next_point.copy()]) - new_speeds.append(speed_history + [speed]) - - particle_positions = new_positions - particle_paths = new_paths - particle_speeds = new_speeds - - frame_path = os.path.join(frame_dir, f"frame_{step_index:06d}.png") - _render_particle_frame( - frame_path, - bounds, - axis_labels, - surrogate_faces, - velocity_norm, - step_index, - particle_positions, - particle_paths, - particle_speeds, - ) +def _save_frame(fig, output_path): + fig.savefig(output_path, dpi=FIG_DPI) def main(): @@ -541,18 +409,22 @@ def main(): return 1 script_dir = os.path.dirname(os.path.abspath(__file__)) - particle_frame_dir = os.path.join(script_dir, "frames_time_particles") - cut_frame_dir = os.path.join(script_dir, "frames_time_cut_contour") - plane_frame_dir = os.path.join(script_dir, "frames_time_plane_x0") - particle_gif_path = os.path.join(script_dir, "particles_trails.gif") - cut_gif_path = os.path.join(script_dir, "cut_contour_x_lt_0.gif") - plane_gif_path = os.path.join(script_dir, "plane_x_eq_0_contour.gif") - - _prepare_output_directory(particle_frame_dir) - _prepare_output_directory(cut_frame_dir) - _prepare_output_directory(plane_frame_dir) - - project_parameters_path = os.path.join(script_dir, "ProjectParameters_3D_fluid.json") + cut_frame_dir = os.path.join(script_dir, CUT_FRAME_DIRNAME) + velocity_plane_frame_dir = os.path.join(script_dir, VELOCITY_PLANE_FRAME_DIRNAME) + pressure_plane_frame_dir = os.path.join(script_dir, PRESSURE_PLANE_FRAME_DIRNAME) + + cut_gif_path = os.path.join(script_dir, CUT_GIF_FILENAME) + velocity_plane_gif_path = os.path.join(script_dir, VELOCITY_PLANE_GIF_FILENAME) + pressure_plane_gif_path = os.path.join(script_dir, PRESSURE_PLANE_GIF_FILENAME) + + if WRITE_3D_CUT_GIF: + _prepare_output_directory(cut_frame_dir) + if WRITE_VELOCITY_PLANE_GIF: + _prepare_output_directory(velocity_plane_frame_dir) + if WRITE_PRESSURE_PLANE_GIF: + _prepare_output_directory(pressure_plane_frame_dir) + + project_parameters_path = os.path.join(script_dir, PROJECT_PARAMETERS_FILENAME) with open(project_parameters_path, "r") as parameter_file: parameters = KratosMultiphysics.Parameters(parameter_file.read()) _normalize_solver_settings(parameters) @@ -563,13 +435,9 @@ def main(): simulation.Initialize() velocity_norm = Normalize(vmin=0.0, vmax=COLORBAR_MAX) - axis_labels = _axis_labels() - surrogate_faces = _collect_surrogate_faces(model) + surrogate_faces = _collect_surrogate_faces(model) if WRITE_3D_CUT_GIF else [] frame_index = 0 - final_points = None - final_velocity = None - final_dt = None while simulation.KeepAdvancingSolutionLoop(): simulation.time = simulation._AdvanceTime() simulation.InitializeSolutionStep() @@ -587,7 +455,7 @@ def main(): break model_part = model[PLOT_MODEL_PART_NAME] - points, velocity, _ = _sample_element_fields(model_part) + points, velocity, pressure = _sample_element_fields(model_part) if points.size == 0: print("No element data found for plotting.") break @@ -599,10 +467,10 @@ def main(): cut_mask = points[:, 0] < 0.0 cut_points = points[cut_mask] - if cut_points.size > 0: + if WRITE_3D_CUT_GIF and cut_points.size > 0: fig_cut = plt.figure(figsize=(9, 6)) ax_cut = fig_cut.add_subplot(1, 1, 1, projection="3d") - _set_3d_axes(ax_cut, bounds, axis_labels) + _set_3d_axes(ax_cut, bounds) ax_cut.set_title(rf"Cut contour $x < 0$ (t = {simulation.time:.5f})") if surrogate_faces: @@ -624,50 +492,108 @@ def main(): ) plt.colorbar(cut_artist, ax=ax_cut, shrink=0.8, pad=0.03, label=r"$|v_z|$") plt.tight_layout() - plt.savefig(os.path.join(cut_frame_dir, f"frame_{frame_index:06d}.png"), dpi=FIG_DPI) + _save_frame(fig_cut, os.path.join(cut_frame_dir, f"frame_{frame_index:06d}.png")) plt.close(fig_cut) - plane_points, plane_values, plane_x = _extract_x_plane_layer(points, np.abs(velocity[:, 2])) - if plane_points is not None and plane_points.size > 0: - fig_plane, ax_plane = plt.subplots(figsize=(8, 6)) - plane_artist = _plot_x_plane(ax_plane, plane_points, plane_values, velocity_norm) - ax_plane.set_title(rf"Plane contour $x \approx {plane_x:.3f}$ (t = {simulation.time:.5f})") - plt.colorbar(plane_artist, ax=ax_plane, shrink=0.9, pad=0.03, label=r"$|v_z|$") - plt.tight_layout() - plt.savefig(os.path.join(plane_frame_dir, f"frame_{frame_index:06d}.png"), dpi=FIG_DPI) - plt.close(fig_plane) + velocity_magnitude = np.linalg.norm(velocity, axis=1) - final_points = points - final_velocity = velocity - final_dt = model_part.ProcessInfo[KratosMultiphysics.DELTA_TIME] if USE_PHYSICAL_DT else None - frame_index += 1 + # ------------------------------------------------------------ + # Velocity magnitude on x ~= 0 plane + # ------------------------------------------------------------ + plane_points, plane_values, plane_x = _extract_x_plane_layer( + points, + velocity_magnitude, + ) - simulation.Finalize() + if WRITE_VELOCITY_PLANE_GIF and plane_points is not None and plane_points.size > 0: + fig_plane, ax_plane = plt.subplots(figsize=PLANE_FIGSIZE) - if final_points is not None and final_velocity is not None: - final_bounds = ( - float(final_points[:, 0].min()), - float(final_points[:, 0].max()), - float(final_points[:, 1].min()), - float(final_points[:, 1].max()), - float(final_points[:, 2].min()), - float(final_points[:, 2].max()), - ) - velocity_norm = Normalize(vmin=0.0, vmax=COLORBAR_MAX) - _create_particle_postprocess_gif( - particle_frame_dir, - final_bounds, - axis_labels, - surrogate_faces, - velocity_norm, - final_points, - final_velocity, - final_dt, + plane_artist = _plot_x_plane( + ax_plane, + plane_points, + plane_values, + velocity_norm, + ) + + ax_plane.set_title( + rf"Velocity magnitude on $x \approx {plane_x:.3f}$, $t = {simulation.time:.5f}$", + fontsize=PLANE_TITLE_SIZE, + ) + + _add_2d_colorbar( + fig_plane, + ax_plane, + plane_artist, + r"$|\mathbf{u}|$", + ) + + fig_plane.tight_layout() + _save_frame( + fig_plane, + os.path.join(velocity_plane_frame_dir, f"frame_{frame_index:06d}.png"), + ) + + plt.close(fig_plane) + + # ------------------------------------------------------------ + # Pressure on x ~= 0 plane + # ------------------------------------------------------------ + pressure_plane_points, pressure_plane_values, pressure_plane_x = _extract_x_plane_layer( + points, + pressure, ) - _write_gif(particle_frame_dir, particle_gif_path) - _write_gif(cut_frame_dir, cut_gif_path) - _write_gif(plane_frame_dir, plane_gif_path) + if WRITE_PRESSURE_PLANE_GIF and pressure_plane_points is not None and pressure_plane_points.size > 0: + pressure_min = float(np.min(pressure_plane_values)) + pressure_max = float(np.max(pressure_plane_values)) + + if abs(pressure_max - pressure_min) < 1.0e-14: + pressure_min -= 0.5 + pressure_max += 0.5 + + pressure_norm = Normalize( + vmin=pressure_min, + vmax=pressure_max, + ) + + fig_pressure, ax_pressure = plt.subplots(figsize=PLANE_FIGSIZE) + + pressure_artist = _plot_x_plane( + ax_pressure, + pressure_plane_points, + pressure_plane_values, + pressure_norm, + ) + + ax_pressure.set_title( + rf"Pressure on $x \approx {pressure_plane_x:.3f}$, $t = {simulation.time:.5f}$", + fontsize=PLANE_TITLE_SIZE, + ) + + _add_2d_colorbar( + fig_pressure, + ax_pressure, + pressure_artist, + r"$p$", + ) + + fig_pressure.tight_layout() + _save_frame( + fig_pressure, + os.path.join(pressure_plane_frame_dir, f"frame_{frame_index:06d}.png"), + ) + + plt.close(fig_pressure) + + frame_index += 1 + + simulation.Finalize() + if WRITE_3D_CUT_GIF: + _write_gif(cut_frame_dir, cut_gif_path) + if WRITE_VELOCITY_PLANE_GIF: + _write_gif(velocity_plane_frame_dir, velocity_plane_gif_path) + if WRITE_PRESSURE_PLANE_GIF: + _write_gif(pressure_plane_frame_dir, pressure_plane_gif_path) return 0