.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_how_to/example_fvc_fvm_operators.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_how_to_example_fvc_fvm_operators.py: Compute finite-volume operators from Python =========================================== :mod:`pybFoam.fvc` and :mod:`pybFoam.fvm` mirror OpenFOAM's C++ operator API. ``fvc`` evaluates explicit operators and returns a ``tmp_*`` field wrapper; ``fvm`` produces implicit matrix terms that can be summed into an equation and solved. This how-to evaluates every ``fvc`` operator the solver API exposes against the shipped VOF case, assembles a momentum matrix with ``fvm``, and plots ``|∇p_rgh|`` on the cell centres. Prerequisite: finished :doc:`/auto_tutorials/example_02_scalar_fields`. .. GENERATED FROM PYTHON SOURCE LINES 19-24 Set up case + mesh ------------------ :func:`pybFoam.clone_example` copies the named example into a tmp dir and restores ``0.orig`` → ``0``. The mesh is generated from ``system/blockMeshDict``. .. GENERATED FROM PYTHON SOURCE LINES 24-34 .. code-block:: Python from pybFoam import Time, argList, clone_example, dictionary, fvMesh from pybFoam.meshing import generate_blockmesh case = clone_example("case") time = Time(argList([str(case), "-case", str(case)])) generate_blockmesh(time, dictionary.read(str(case / "system" / "blockMeshDict"))) mesh = fvMesh(time) print(f"mesh: {mesh.nCells()} cells") .. rst-class:: sphx-glr-script-out .. code-block:: none mesh: 2268 cells .. GENERATED FROM PYTHON SOURCE LINES 35-37 Read fields ----------- .. GENERATED FROM PYTHON SOURCE LINES 37-64 .. code-block:: Python import numpy as np from pybFoam import createPhi, volScalarField, volVectorField, write p_rgh = volScalarField.read_field(mesh, "p_rgh") U = volVectorField.read_field(mesh, "U") # The shipped IC is uniform zero, which would leave grad / div / lap # all at numerical noise. Seed an analytic sinusoidal pressure and a # matching divergence-free velocity so the operator outputs are # meaningful (and the gradient plot has visible structure). C = np.asarray(mesh.C()["internalField"]) L = 0.584 # case extent (scale 0.146 × 4) np.asarray(p_rgh["internalField"])[:] = ( 1000.0 * np.sin(np.pi * C[:, 0] / L) * np.sin(np.pi * C[:, 1] / L) ) np.asarray(U["internalField"])[:] = np.column_stack( [np.sin(np.pi * C[:, 1] / L), -np.sin(np.pi * C[:, 0] / L), np.zeros_like(C[:, 0])] ) p_rgh.correctBoundaryConditions() U.correctBoundaryConditions() write(p_rgh) write(U) phi = createPhi(U) .. GENERATED FROM PYTHON SOURCE LINES 65-70 Explicit operators (``fvc``) ---------------------------- Each ``fvc`` call returns a ``tmp_*`` reference so OpenFOAM can reuse the underlying buffer — calling it (``()``) materialises the concrete field. Pass ``tmp_*`` directly into another operator to chain. .. GENERATED FROM PYTHON SOURCE LINES 70-85 .. code-block:: Python from pybFoam import fvc grad_p = fvc.grad(p_rgh)() div_U = fvc.div(U)() lap_p = fvc.laplacian(p_rgh)() print(f"grad(p_rgh) shape : {np.asarray(grad_p['internalField']).shape}") print( f"div(U) min/max: " f"{np.asarray(div_U['internalField']).min():.3e} / " f"{np.asarray(div_U['internalField']).max():.3e}" ) print(f"lap(p_rgh) sum : {np.asarray(lap_p['internalField']).sum():.3e}") .. rst-class:: sphx-glr-script-out .. code-block:: none grad(p_rgh) shape : (2268, 3) div(U) min/max: -1.678e+02 / 7.874e+01 lap(p_rgh) sum : 9.790e+06 .. GENERATED FROM PYTHON SOURCE LINES 86-88 Flux + surface-normal gradient + reconstruct -------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 88-102 .. code-block:: Python flux = fvc.flux(U)() # surfaceScalarField div_phiU = fvc.div(flux, U)() snGrad_p = fvc.snGrad(p_rgh)() reconstructed = fvc.reconstruct(flux)() flux_vals = np.asarray(flux["internalField"]) snGrad_vals = np.asarray(snGrad_p["internalField"]) recon_vals = np.asarray(reconstructed["internalField"]) print(f"flux : surfaceScalarField, {flux_vals.shape[0]} internal faces") print(f"snGrad_p : {snGrad_vals.shape[0]} internal face values") print(f"reconstructed: volVectorField, {recon_vals.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none flux : surfaceScalarField, 4432 internal faces snGrad_p : 4432 internal face values reconstructed: volVectorField, (2268, 3) .. GENERATED FROM PYTHON SOURCE LINES 103-108 Implicit/matrix operators (``fvm``) ----------------------------------- ``fvm`` operators produce terms of an equation matrix rather than field values. Combine them with ``+`` / ``-`` and wrap into an ``fvScalarMatrix`` or ``fvVectorMatrix``. .. GENERATED FROM PYTHON SOURCE LINES 108-127 .. code-block:: Python from pybFoam import ( Word, dimensionedScalar, dimLength, dimTime, fvm, fvScalarMatrix, fvVectorMatrix, ) pEqn = fvScalarMatrix(fvm.laplacian(p_rgh)) nu = dimensionedScalar(Word("nu"), dimLength * dimLength / dimTime, 1.0e-5) UEqn = fvVectorMatrix(fvm.ddt(U) + fvm.div(phi, U) - fvm.laplacian(nu, U)) print(f"pEqn: {type(pEqn).__name__}") print(f"UEqn: {type(UEqn).__name__}") .. rst-class:: sphx-glr-script-out .. code-block:: none pEqn: fvScalarMatrix UEqn: fvVectorMatrix .. GENERATED FROM PYTHON SOURCE LINES 128-135 Visualize ``∇p_rgh`` on the cell centres ---------------------------------------- ``mesh.C()`` gives the cell-centre ``volVectorField``. Combined with ``np.asarray(grad_p['internalField'])`` we get a (n_cells, 3) array of the gradient and (n_cells, 3) cell coordinates. The case is 2-D (single layer in z), so a matplotlib ``tricontourf`` of the in-plane gradient magnitude is enough. .. GENERATED FROM PYTHON SOURCE LINES 135-166 .. code-block:: Python import matplotlib.pyplot as plt C = np.asarray(mesh.C()["internalField"]) g = np.asarray(grad_p["internalField"]) mag_g = np.linalg.norm(g[:, :2], axis=1) fig, ax = plt.subplots(figsize=(6, 4)) tcf = ax.tricontourf(C[:, 0], C[:, 1], mag_g, levels=20, cmap="magma") fig.colorbar(tcf, ax=ax, label=r"$|\nabla p_{rgh}|$") # Quiver on a coarse subsample so arrows are readable. n_cells = C.shape[0] step = max(1, n_cells // 400) ax.quiver( C[::step, 0], C[::step, 1], g[::step, 0], g[::step, 1], color="white", scale=None, scale_units="xy", width=0.002, ) ax.set_aspect("equal") ax.set_xlabel("x [m]") ax.set_ylabel("y [m]") ax.set_title(r"$\nabla p_{rgh}$ from fvc.grad") fig.tight_layout() plt.show() .. image-sg:: /auto_how_to/images/sphx_glr_example_fvc_fvm_operators_001.png :alt: $\nabla p_{rgh}$ from fvc.grad :srcset: /auto_how_to/images/sphx_glr_example_fvc_fvm_operators_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 167-181 Mixing implicit and explicit terms ---------------------------------- A typical PISO/SIMPLE loop builds the momentum matrix with ``fvm``, adds an explicit pressure gradient from ``fvc``, and hands the system to ``solve``: .. code-block:: python if piso.momentumPredictor(): solve(UEqn + fvc.grad(p)) The shipped case uses ``p_rgh`` (full pressure) whose dimensions do not match ``UEqn``, so we stop at showing the matrix shapes. For a worked solver see ``examples/cavity/icoFoam.py``. .. GENERATED FROM PYTHON SOURCE LINES 183-188 See also -------- - :doc:`/auto_tutorials/example_04_run_and_sample` — feed ``UEqn`` into the icoFoam PISO loop and sample the result. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.203 seconds) .. _sphx_glr_download_auto_how_to_example_fvc_fvm_operators.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_fvc_fvm_operators.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_fvc_fvm_operators.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_fvc_fvm_operators.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_