.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_tutorials/example_03_build_a_solver.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_tutorials_example_03_build_a_solver.py: Build a scalar-transport solver from scratch ============================================ Build a complete NeoFOAM solver for unsteady scalar transport with a prescribed velocity field: .. math:: \frac{\partial T}{\partial t} + \nabla \cdot (U\, T) - \nabla \cdot (D \nabla T) = 0 The smallest PDE that exercises every major framework piece — time loop, field reads, dependency graph, ``fvm`` solve, topological ordering — without pressure-velocity coupling. Complete :doc:`example_02_passive_scalar_plugin` first; this tutorial reuses ``field()``, ``depends_on``, and parameter injection from there. .. GENERATED FROM PYTHON SOURCE LINES 24-26 Solver imports -------------- .. GENERATED FROM PYTHON SOURCE LINES 26-67 .. code-block:: Python import contextlib import subprocess from typing import Annotated, Any, Optional import pybFoam as pyf import pyvista as pv from pybFoam import ( Info, fvm, fvScalarMatrix, surfaceScalarField, volScalarField, volVectorField, ) from neofoam import Depends, FieldUpdates, Solver, field from neofoam.framework.context import Context from neofoam.framework.graph import DAGResolver from neofoam.framework.initialization import ( ConfigContext, InitializerBuilder, InitStep, LoadResult, StagedInitRunner, StagedInitSpec, lazy, ) from neofoam.framework.initialization.execution.ordering import ( _topological_sort as topological_sort, ) from neofoam.framework.operations import ( IterativeOp, Operation, Operations, StepBuilder, ) from neofoam.framework.types import OperationMetadata from neofoam.tutorial import clone_case .. GENERATED FROM PYTHON SOURCE LINES 68-74 Define the spec and the time-loop predicate ------------------------------------------- ``Solver("scalar_transport")`` returns an empty ``SolverSpec`` — the analogue of ``Model`` from tutorial 2. ``TimeLoop`` is the iteration condition for the outer loop; ``IterativeOp`` runs it until it returns False. .. GENERATED FROM PYTHON SOURCE LINES 74-84 .. code-block:: Python class TimeLoop: def __call__(self, ctx: Context) -> bool: return bool(ctx.runtime.run()) scalar_transport = Solver("scalar_transport") .. GENERATED FROM PYTHON SOURCE LINES 85-98 Build the staged initializer ---------------------------- ``create_init`` returns a ``StagedInitRunner``. Its three stages are registered as decorators on a per-call ``StagedInitSpec`` builder and defined as closures so ``build`` can read ``runner.argv`` (the framework sets ``runner.argv`` after ``create_init`` returns but before ``runner.run()``). This solver has no plugin models, so LOAD and RESOLVE are trivial. BUILD creates the runtime + mesh, then registers one ``field(...)`` per object with ``depends_on`` listing prerequisite step names. ``create_phi`` reads ``ctx["fields.U"]`` (the prefixed name) because ``field("U", ...)`` produces an ``InitStep`` named ``fields.U``. .. GENERATED FROM PYTHON SOURCE LINES 98-145 .. code-block:: Python def create_init(case_dir: Optional[Any] = None) -> StagedInitRunner: spec_builder = StagedInitSpec.build("scalar_transport") @spec_builder.load def load_config() -> LoadResult: return LoadResult(core_models=[], optional_models=[]) @spec_builder.resolve def resolve_models(config: ConfigContext) -> None: pass @spec_builder.build def build_lazy( core_models: list[Any], optional_models: list[Any] ) -> list[InitStep]: argv = runner.argv builder = InitializerBuilder() builder.add(lazy("runtime", lambda _ctx: pyf.Time(pyf.argList(argv)))) builder.add( lazy("mesh", lambda ctx: pyf.fvMesh(ctx["runtime"]), depends_on=["runtime"]) ) def create_T(ctx: dict[str, Any]) -> volScalarField: return volScalarField.read_field(ctx["mesh"], "T") def create_U(ctx: dict[str, Any]) -> volVectorField: return volVectorField.read_field(ctx["mesh"], "U") def create_D(ctx: dict[str, Any]) -> volScalarField: return volScalarField.read_field(ctx["mesh"], "D") def create_phi(ctx: dict[str, Any]) -> surfaceScalarField: return pyf.createPhi(ctx["fields.U"]) builder.add(field("T", create_T, depends_on=["mesh"])) builder.add(field("U", create_U, depends_on=["mesh"])) builder.add(field("D", create_D, depends_on=["mesh"])) builder.add(field("phi", create_phi, depends_on=["fields.U"])) return builder.build() runner = StagedInitRunner(spec_builder.finalize()) return runner .. GENERATED FROM PYTHON SOURCE LINES 146-151 Register the initializer ------------------------ Every solver has exactly one ``@initializer``. The ``Annotated[..., Depends(create_init)]`` marker tells the framework to call ``create_init()`` and pass its return value as ``init_obj``. .. GENERATED FROM PYTHON SOURCE LINES 151-161 .. code-block:: Python @scalar_transport.initializer def initialize( self: Any, init_obj: Annotated[StagedInitRunner, Depends(create_init)], ) -> Context: return init_obj.run() .. GENERATED FROM PYTHON SOURCE LINES 162-167 Define the execution graph -------------------------- ``builder.loop(...)`` returns a fresh ``StepBuilder`` scoped to the new loop's body. We return an empty ``Operations()`` because no plugin contributes operations to this solver. .. GENERATED FROM PYTHON SOURCE LINES 167-189 .. code-block:: Python @scalar_transport.execution_graph_step def execution_graph( self: Any, domain_name: Optional[str] = None, ) -> tuple[StepBuilder, Operations]: ops = self.operations builder = StepBuilder() time_loop_op = Operation( func=IterativeOp(TimeLoop()), metadata=OperationMetadata(op_name="time_loop"), ) with builder.loop(time_loop_op) as time_builder: time_builder.step(ops["increment_time"]) time_builder.step(ops["solve_T"]) time_builder.step(ops["write_output"]) return builder, Operations() .. GENERATED FROM PYTHON SOURCE LINES 190-195 Add the solver's operations --------------------------- Three operations: increment time, solve ``T``, write output. Parameter-name injection on ``solve_T`` pulls ``T``, ``U``, ``phi``, ``D`` from ``ctx.fields``. .. GENERATED FROM PYTHON SOURCE LINES 195-222 .. code-block:: Python @scalar_transport.operation() def increment_time(self: Any, ctx: Context) -> None: Info(f"Time = {ctx.runtime.timeName()}") ctx.runtime.increment() @scalar_transport.operation(depends_on=["increment_time"]) def solve_T( self: Any, T: volScalarField, U: volVectorField, phi: surfaceScalarField, D: volScalarField, ) -> FieldUpdates: TEqn = fvScalarMatrix(fvm.ddt(T) + fvm.div(phi, T) - fvm.laplacian(D, T)) TEqn.solve() return FieldUpdates({"T": T}) @scalar_transport.operation(depends_on=["solve_T"]) def write_output(self: Any, ctx: Context) -> None: ctx.runtime.write(True) ctx.runtime.printExecutionTime() .. GENERATED FROM PYTHON SOURCE LINES 223-230 Read a graph error ------------------ Worth seeing the resolver fail once. If you accidentally make ``U`` depend on ``phi`` while ``phi`` already depends on ``U``, the cycle detector refuses to run. The fix is to delete the bogus edge: ``U`` is read from disk and depends only on ``mesh``; ``phi`` *derives from* ``U``, not the other way around. .. GENERATED FROM PYTHON SOURCE LINES 230-240 .. code-block:: Python cyclic_steps = [ field("U", lambda ctx: None, depends_on=["fields.phi"]), field("phi", lambda ctx: None, depends_on=["fields.U"]), ] try: topological_sort(cyclic_steps) except Exception as exc: print(f"caught: {type(exc).__name__}: {exc}") .. rst-class:: sphx-glr-script-out .. code-block:: none caught: InitializationGraphError: Circular dependency: fields.U -> fields.phi .. GENERATED FROM PYTHON SOURCE LINES 241-246 Run the new solver on the bundled case -------------------------------------- ``tutorials/scalar_transport_min`` is a 50×50 unit-square mesh with uniform inflow ``U=(0.5, 0, 0)`` and a tracer-1 inlet driving advection-diffusion of ``T``. .. GENERATED FROM PYTHON SOURCE LINES 246-259 .. code-block:: Python case = clone_case("scalar_transport_min") (case / "scalar_transport_min.foam").touch() subprocess.run(["blockMesh", "-case", str(case)], check=True) with contextlib.chdir(case): solver = scalar_transport.instantiate(argv=["scalar_transport"]) ctx = solver.initialize() builder, model_ops = solver.execution_graph() resolver = DAGResolver() resolved = resolver.resolve(builder, model_ops) resolved.operations.run(ctx) .. GENERATED FROM PYTHON SOURCE LINES 260-265 Animate the diffusion-advection front ------------------------------------- Watch ``T`` sweep across the domain as advection pushes the inlet's ``T = 1`` boundary condition downstream and diffusion smears the front. A static frame would mean ``solve_T`` never ran. .. GENERATED FROM PYTHON SOURCE LINES 265-295 .. code-block:: Python reader = pv.OpenFOAMReader(str(case / "scalar_transport_min.foam")) pl = pv.Plotter(off_screen=True, window_size=(700, 700)) pl.open_gif(str(case / "scalar_transport.gif"), fps=8) reader.set_active_time_value(reader.time_values[0]) mesh = reader.read()["internalMesh"] pl.add_mesh( mesh, scalars="T", cmap="inferno", clim=[0, 1], scalar_bar_args={ "title": "T [-]", "vertical": False, "position_x": 0.2, "position_y": 0.05, "width": 0.6, "height": 0.04, }, ) pl.view_xy() for t in reader.time_values: reader.set_active_time_value(t) mesh.copy_from(reader.read()["internalMesh"]) pl.write_frame() pl.show() .. image-sg:: /auto_tutorials/images/sphx_glr_example_03_build_a_solver_001.gif :alt: example 03 build a solver :srcset: /auto_tutorials/images/sphx_glr_example_03_build_a_solver_001.gif :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.328 seconds) .. _sphx_glr_download_auto_tutorials_example_03_build_a_solver.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_03_build_a_solver.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_03_build_a_solver.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_03_build_a_solver.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_