.. _UserDefinedPostprocessing:

.. currentmodule:: flow360

User Defined Postprocessing via :class:`UserDefinedField`
===================================================================

In Flow360, users can specify custom expressions to calculate outputs based on some solver variables. These variables can be used in the :code:`output_fields` of :class:`~VolumeOutput`, :class:`~SurfaceOutput`, :class:`~SliceOutput`, :class:`~IsosurfaceOutput` and :class:`~ProbeOutput`.

We will illustrate how to use this feature using several examples.

.. note::

   We provide a complete list of available functions and some guidelines in our :ref:`knowledge base<userDefinedExpressionsKnowledgeBase>`.

.. _UDFTotPressure:

Custom Output Variable
--------------------------------------------------------------------------

In this example, we compute the total pressure coefficient for the entire domain of an :ref:`om6Wing<quick_start_geometry>` simulation. For the definition of the total pressure coefficient please see :ref:`here<totalPressureCoefficientDefinition>`.

First we need to define the expression for total pressure coefficient.

.. code-block:: python
   :linenos:

   fl.UserDefinedField(
      name="TotalPressureCoeff",
      expression="double gamma = 1.40; double pow1 = gamma/(gamma-1); double pow2 = (gamma-1) / 2; double MachRefSq = MachRef * MachRef; "
      + "double Mach = sqrt(primitiveVars[1] * primitiveVars[1] + primitiveVars[2] * primitiveVars[2] + primitiveVars[3] * primitiveVars[3]) / sqrt(gamma * primitiveVars[4] / primitiveVars[0]);  double MachSq = Mach * Mach; "
      + "TotalPressureCoeff = (gamma*primitiveVars[4]*pow((1+pow2*MachSq),pow1)-pow((1+pow2*MachRefSq),pow1))/(gamma/2*MachRefSq);",
   )


A statement-by-statement breakdown of the above expression is shown in following table:

.. _UserDefinedFields_:

.. table::
   :widths: 40 50

   +---------------------------------------------------+-------------------------------------------------------------------+
   | Expression                                        | Explanation                                                       |
   +===================================================+===================================================================+
   |  :code:`double gamma = 1.40;....`                 | A declaration of local variables that can be used \               |
   |                                                   | in subsequent expression.                                         |
   +---------------------------------------------------+-------------------------------------------------------------------+
   | :code:`double Mach = sqrt(primit....`\            | Still a declaration. But the Mach number is calculated locally  \ |
   |                                                   | even though it is a valid option for :code:`outputFields`.       \|
   |                                                   | This is an example that not all the pre-defined output variables \|
   |                                                   | are available to be used in the expression.                       |
   +---------------------------------------------------+-------------------------------------------------------------------+
   | :code:`TotalPressureCoeff = (gamma*pri....`\      | This assigns the final result to the total pressure coefficient.\ |
   |                                                   | Note that the output variable name in the expression            \ |
   |                                                   | (:code:`TotalPressureCoeff`) is exactly                         \ |
   |                                                   | the same as its :code:`name` entry in the JSON.                 \ |
   +---------------------------------------------------+-------------------------------------------------------------------+

Then we need to add :code:`TotalPressureCoeff` in one of the outputs. Here we use :class:`~VolumeOutput`.

.. code-block:: python
   :linenos:

   fl.VolumeOutput(
      output_format="paraview",
      output_fields=["TotalPerssureCoeff"],
   )

.. _UDFSurfIntegral:

Custom surface integral
--------------------------------------------------------------------------

In Flow360 we have added a new type of monitor :class:`~SurfaceIntegralOutput` in addition to the probe monitor which monitors quantities at a set of given coordinates. The surface-integral monitor will try to compute the surface integral of the given variable on given patches. Here we show an example of calculating the pressure force on a no-slip wall as an example. We first set up the :py:attr:`~SimulationParams.user_defined_fields` as 

.. code-block:: python
   :linenos:
   
   pressure_force_field = fl.UserDefinedField(
      name="PressureForce",
      expression="double prel = primitiveVars[4] - pressureFreestream; "
      + "PressureForce[0] = prel * nodeNormals[0]; "
      + "PressureForce[1] = prel * nodeNormals[1]; "
      + "PressureForce[2] = prel * nodeNormals[2];",
   )

   fl.SimulationParams(
      ...
      user_defined_fields=[pressure_force_field]
   )

Note that the :code:`nodeNormals` is a vector whose direction is the normal of the given patch at each node and the magnitude is the area associated with the node. Therefore the right hand side of the :code:`PressureForce` corresponds to the red part of the following equation for pressure force:

.. math::
   :label: def_pressureForce

   \overrightarrow{\mathbf{f_p}} = \int {\color{Red} \left(p - p_{ref}\right )d\overrightarrow{\mathbf{A}}}

where the :math:`\overrightarrow{\mathbf{f_p}}` is the pressure force, :math:`p_{ref}` is the free steam pressure. Surface integrals for scalar outputs should use :code:`magnitude(nodeNormals)` instead to include the area (:math:`dA`) in integral. Otherwise the result will be a summation of the value over the nodes instead of surface area integral.

And we append to the :py:attr:`~SimulationParams.outputs` to include:

.. code-block:: python
   :linenos:

   fl.SurfaceIntegralOutput(
      surfaces=[volume_mesh["wing1"], volume_mesh["wing2"]], output_fields="PressureForce"
   )

which means that we want :code:`PressureForce` integral computed on patch :code:`wing1` and patch :code:`wing2` as a whole.