<commit_before># -*- coding: utf-8 -*-
#
# TODO update
'''
Numerical solution schemes for the Navier--Stokes equation in cylindrical
coordinates,

.. math::
    \\DeclareMathOperator{\\div}{div}

.. math::
    \\begin{align*}
        &\\rho \\left(\\frac{du}{dt} + (u\\cdot\\nabla)u\\right)
          = -\\nabla p
            + \\mu \\left(
                \\frac{1}{r} \\div(r \\nabla u)
                - e_r \\frac{u_r}{r^2}
                \\right)
            + f,\\\\
        &\\frac{1}{r} \\div(r u) = 0,
    \\end{align*}

cf.
https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations#Representations_in_3D.
In the weak formulation, we consider integrals in pseudo 3D, resulting in a
weighting with :math:`2\\pi r` of the equations. (The volume element is
:math:`2\\pi r\\,\\text{d}x`.)

The order of the variables is taken to be :math:`(r, z, \\theta)`. This makes
sure that for planar domains, the :math:`x`- and :math:`y`-coordinates are
interpreted as :math:`r`, :math:`z`.
'''

from dolfin import (
    TestFunction, Function, Constant, dot, grad, inner, pi, dx, div, solve,
    derivative, TrialFunction, PETScPreconditioner, PETScKrylovSolver,
    as_backend_type, info, assemble, norm, FacetNormal, sqrt, ds, as_vector,
    NonlinearProblem, NewtonSolver, SpatialCoordinate, project
    )

from . import stabilization as stab
from .message import Message


def _momentum_equation(u, v, p, f, rho, mu, stabilization, my_dx):
    '''Weak form of the momentum equation.
    '''
    # rho and my are Constant() functions
    assert rho.values()[0] > 0.0
    assert mu.values()[0] > 0.0

    # Skew-symmetric formulation.
    # Don't include the boundary term
    #
    #   - mu *inner(r*grad(u2)*n  , v2) * 2*pi*ds.
    #
    # This effectively means that at all boundaries where no sufficient
    # Dirichlet-conditions are posed, we assume grad(u)*n to vanish.
    #
    # The original term
    #    u2[0]/(r*r) * v2[0]
    # doesn't explode iff u2[0]~r and v2[0]~r at r=0. Hence, we need to enforce
    # homogeneous Dirichlet-conditions for n.u at r=0. This corresponds to no
    # flow in normal direction through the symmetry axis -- makes sense.
    # When using the 2*pi*r weighting, we can even be a bit lax on the
    # enforcement on u2[0].
    #
    # For this to be well defined, u[0]/r and u[2]/r must be bounded for r=0,
    # so u[0]~u[2]~r must hold. This either needs to be enforced in the
    # boundary conditions (homogeneous Dirichlet for u[0], u[2] at r=0) or must
    # follow from the dynamics of the system.
    #
    # TODO some more explanation for the following lines of code
    mesh = v.function_space().mesh()
    r = SpatialCoordinate(mesh)[0]
    F = rho * 0.5 * (dot(grad(u) * u, v) - dot(grad(v) * u, u)) \
            * 2*pi*r*my_dx \
        + mu * inner(r * grad(u), grad(v)) * 2 * pi * my_dx  \
        + mu * u[0] / r * v[0] * 2 * pi * my_dx  \
        - dot(f, v) * 2*pi*r * my_dx
    if p:
        F += (p.dx(0) * v[0] + p.dx(1) * v[1]) * 2*pi*r * my_dx
    if len(u) == 3:
        F += rho * (-u[2] * u[2] * v[0] + u[0] * u[2] * v[2]) * 2*pi * my_dx
        F += mu * u[2] / r * v[2] * 2 * pi * my_dx

    if stabilization == 'SUPG':
        # TODO check this part of the code
        #
        # SUPG stabilization has the form
        #
        #     <R, tau*grad(v)*u[0]>
        #
        # with R being the residual in strong form. The choice of tau is
        # subject to research.
        rho_val = rho.values()[0]
        mu_val = mu.values()[0]
        tau = stab.supg2(
                u.function_space().mesh(),
                u,
                mu_val / rho_val,
                u.function_space().ufl_element().degree()
                )
        # Strong residual:
        R = + rho * grad(u) * u * 2*pi*r \
            - mu * div(r * grad(u)) * 2*pi \
            - f * 2*pi*r
        if p:
            R += grad(p) * 2*pi*r

        gv = tau * grad(v) * u
        F += dot(R, gv) * my_dx

        # We need to deal with the term
        #
        #     \int mu * (u2[0]/r**2, 0) * dot(R, grad(v2)*b_tau) 2*pi*r*dx
        #
        # somehow. Unfortunately, it's not easy to construct (u2[0]/r**2,
        # 0), cf.  <https://answers.launchpad.net/dolfin/+question/228353>.
        F += mu * u[0] / r * 2*pi * gv[0] * my_dx
        if u.function_space().num_sub_spaces() == 3:
            F += rho * (-u[2] * u[2] * gv[0] + u[0] * u[2] * gv[2]) \
                    * 2*pi*my_dx<commit_msg>fix for navier-stokes test<commit_after># -*- coding: utf-8 -*-
#
# TODO update
'''
Numerical solution schemes for the Navier--Stokes equation in cylindrical
coordinates,

.. math::
    \\DeclareMathOperator{\\div}{div}

.. math::
    \\begin{align*}
        &\\rho \\left(\\frac{du}{dt} + (u\\cdot\\nabla)u\\right)
          = -\\nabla p
            + \\mu \\left(
                \\frac{1}{r} \\div(r \\nabla u)
                - e_r \\frac{u_r}{r^2}
                \\right)
            + f,\\\\
        &\\frac{1}{r} \\div(r u) = 0,
    \\end{align*}

cf.
https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations#Representations_in_3D.
In the weak formulation, we consider integrals in pseudo 3D, resulting in a
weighting with :math:`2\\pi r` of the equations. (The volume element is
:math:`2\\pi r\\,\\text{d}x`.)

The order of the variables is taken to be :math:`(r, z, \\theta)`. This makes
sure that for planar domains, the :math:`x`- and :math:`y`-coordinates are
interpreted as :math:`r`, :math:`z`.
'''

from dolfin import (
    TestFunction, Function, Constant, dot, grad, inner, pi, dx, div, solve,
    derivative, TrialFunction, PETScPreconditioner, PETScKrylovSolver,
    as_backend_type, info, assemble, norm, FacetNormal, sqrt, ds, as_vector,
    NonlinearProblem, NewtonSolver, SpatialCoordinate, project
    )

from . import stabilization as stab
from .message import Message


def _momentum_equation(u, v, p, f, rho, mu, stabilization, my_dx):
    '''Weak form of the momentum equation.
    '''
    # rho and my are Constant() functions
    assert rho.values()[0] > 0.0
    assert mu.values()[0] > 0.0

    # Skew-symmetric formulation.
    # Don't include the boundary term
    #
    #   - mu *inner(r*grad(u2)*n  , v2) * 2*pi*ds.
    #
    # This effectively means that at all boundaries where no sufficient
    # Dirichlet-conditions are posed, we assume grad(u)*n to vanish.
    #
    # The original term
    #    u2[0]/(r*r) * v2[0]
    # doesn't explode iff u2[0]~r and v2[0]~r at r=0. Hence, we need to enforce
    # homogeneous Dirichlet-conditions for n.u at r=0. This corresponds to no
    # flow in normal direction through the symmetry axis -- makes sense.
    # When using the 2*pi*r weighting, we can even be a bit lax on the
    # enforcement on u2[0].
    #
    # For this to be well defined, u[0]/r and u[2]/r must be bounded for r=0,
    # so u[0]~u[2]~r must hold. This either needs to be enforced in the
    # boundary conditions (homogeneous Dirichlet for u[0], u[2] at r=0) or must
    # follow from the dynamics of the system.
    #
    # TODO some more explanation for the following lines of code
    mesh = v.function_space().mesh()
    r = SpatialCoordinate(mesh)[0]
    F = rho * 0.5 * (dot(grad(u) * u, v) - dot(grad(v) * u, u)) \
            * 2*pi*r*my_dx \
        + mu * inner(r * grad(u), grad(v)) * 2 * pi * my_dx  \
        + mu * u[0] / r * v[0] * 2 * pi * my_dx  \
        - dot(f, v) * 2*pi*r * my_dx
    if p:
        F += (p.dx(0) * v[0] + p.dx(1) * v[1]) * 2*pi*r * my_dx
    if len(u) == 3:
        F += rho * (-u[2] * u[2] * v[0] + u[0] * u[2] * v[2]) * 2*pi * my_dx
        F += mu * u[2] / r * v[2] * 2 * pi * my_dx

    if stabilization == 'SUPG':
        # TODO check this part of the code
        #
        # SUPG stabilization has the form
        #
        #     <R, tau*grad(v)*u[0]>
        #
        # with R being the residual in strong form. The choice of tau is
        # subject to research.
        rho_val = rho.values()[0]
        mu_val = mu.values()[0]
        tau = stab.supg(
                u.function_space().mesh(),
                u,
                mu_val / rho_val,
                u.function_space().ufl_element().degree()
                )
        # Strong residual:
        R = + rho * grad(u) * u * 2*pi*r \
            - mu * div(r * grad(u)) * 2*pi \
            - f * 2*pi*r
        if p:
            R += grad(p) * 2*pi*r

        gv = tau * grad(v) * u
        F += dot(R, gv) * my_dx

        # We need to deal with the term
        #
        #     \int mu * (u2[0]/r**2, 0) * dot(R, grad(v2)*b_tau) 2*pi*r*dx
        #
        # somehow. Unfortunately, it's not easy to construct (u2[0]/r**2,
        # 0), cf.  <https://answers.launchpad.net/dolfin/+question/228353>.
        F += mu * u[0] / r * 2*pi * gv[0] * my_dx
        if u.function_space().num_sub_spaces() == 3:
            F += rho * (-u[2] * u[2] * gv[0] + u[0] * u[2] * gv[2]) \
                    * 2*pi*my_dx
{'fraction_non_alphanumeric': 0.15567724631202504, 'fraction_numerical': 0.02313366115333035, 'mean_word_length': 2.0243325447786416, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 11, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1330299', 'n_tokens_mistral': 3377, 'n_tokens_neox': 3178, 'n_words': 950}