utils.py

# Copyright (C) 2021-2023 Jørgen S. Dokken and Igor A. Baratta
#
# SPDX-License-Identifier:    MIT

from typing import Dict, Tuple

from mpi4py import MPI

import basix.ufl
import numpy as np
import ufl
from dolfinx import cpp, default_scalar_type, fem

from generate_team30_meshes import mesh_parameters, model_parameters, surface_map

__all__ = ["DerivedQuantities2D", "update_current_density"]


def _cross_2D(A, B):
    """Compute cross of two 2D vectors"""
    return A[0] * B[1] - A[1] * B[0]


class DerivedQuantities2D:
    """
    Collection of methods for computing derived quantities used in the TEAM 30 benchmark including:
    - Torque of rotor (using classical surface calculation and Arkkio's method)
    - Loss in the rotor (steel and aluminium component separately)
    - Induced voltage in one copper winding
    """

    def __init__(
        self,
        Az: fem.Function,
        Azn: fem.Function,
        u,
        sigma: fem.Function,
        domains: dict,
        ct: cpp.mesh.MeshTags_int32,
        ft: cpp.mesh.MeshTags_int32,
        form_compiler_options: dict = {},
        jit_parameters: dict = {},
    ):
        """
        Parameters
        ==========
        Az            The mixed function of the magnetic vector potential Az

        An
            The mixed function of the magnetic vector potential Az

        u
            Rotational velocity (Expressed as an ufl expression)

        sigma
            Conductivity

        domains
            dictonary were each key indicates a material in the problem. Each item is a tuple of
            indices relating to the volume tags ct and facet tags

        ct
            Meshtag containing cell indices
        ft
            Meshtag containing facet indices

        form_compiler_options
            Parameters used in FFCx compilation of this form. Run `ffcx --help` at
            the commandline to see all available options. Takes priority over all
            other parameter values, except for `scalar_type` which is determined by
            DOLFINx.

        jit_parameters
            Parameters used in CFFI JIT compilation of C code generated by FFCx.
            See `python/dolfinx/jit.py` for all available parameters.
            Takes priority over all other parameter values.

        """
        self.mesh = Az.function_space.mesh
        self.comm = self.mesh.comm

        # Functions
        self.sigma = sigma

        # Constants
        self.dt = fem.Constant(self.mesh, default_scalar_type(0))
        self.L = 1  # Depth of domain (for torque and voltage calculations)

        # Integration quantities
        x = ufl.SpatialCoordinate(self.mesh)
        r = ufl.sqrt(x[0] ** 2 + x[1] ** 2)
        self.domains = domains
        self.dx = ufl.Measure("dx", domain=self.mesh, subdomain_data=ct)
        self.dS = ufl.Measure("dS", domain=self.mesh, subdomain_data=ft)

        # Interior facet restriction (1 on interior airgap, 0 on exterior airgap)
        V_c = fem.functionspace(self.mesh, ("DG", 0))
        gap_markers = domains["AirGap"]
        self._restriction = fem.Function(V_c)
        self._restriction.interpolate(
            lambda x: np.ones(x.shape[1], dtype=default_scalar_type), cells0=ct.find(gap_markers[1])
        )
        self._restriction.interpolate(
            lambda x: np.zeros(x.shape[1], dtype=default_scalar_type),
            cells0=ct.find(gap_markers[0]),
        )
        self._restriction.x.scatter_forward()

        # Derived quantities
        self.B = ufl.as_vector((Az.dx(1), -Az.dx(0)))  # Electromagnetic field
        self.Bphi = ufl.inner(self.B, ufl.as_vector((-x[1], x[0]))) / r
        self.Br = ufl.inner(self.B, x) / r
        self.E = -(Az - Azn) / self.dt  # NOTE: as grad(V)=dV/dz=0 in 2D (-ufl.grad(V)) is excluded
        self.Ep = self.E + _cross_2D(u, self.B)

        # Parameters
        self.fp = form_compiler_options
        self.jp = jit_parameters

        self._init_voltage()
        self._init_loss()
        self._init_torque()

    def _init_voltage(self):
        """
        Initializer for computation of induced voltage in for each the copper winding
        (phase A and -A)
        """
        N = 1  # Number of turns in winding
        if len(self.domains["Cu"]) == 2:
            windings = self.domains["Cu"]
        elif len(self.domains["Cu"]) == 6:
            windings = [
                self.domains["Cu"][0],
                self.domains["Cu"][2],
            ]  # NOTE: assumption on ordering of input windings
        else:
            raise RuntimeError("Only single or three phase computations implemented")
        self._C = []
        self._voltage = []
        for winding in windings:
            self._C.append(
                N
                * self.L
                / self.comm.allreduce(
                    fem.assemble_scalar(fem.form(1 * self.dx(winding))), op=MPI.SUM
                )
            )
            self._voltage.append(
                fem.form(
                    self.E * self.dx(winding), form_compiler_options=self.fp, jit_options=self.jp
                )
            )

    def compute_voltage(self, dt: float):
        """
        Compute induced voltage between two time steps of distance dt
        """
        self.dt.value = dt
        voltages = [self.comm.allreduce(fem.assemble_scalar(voltage)) for voltage in self._voltage]
        return [voltages[i] * self._C[i] for i in range(len(voltages))]

    def _init_loss(self):
        """
        Compute the Loss in the rotor, total and steel component.
        """
        # Induced voltage
        q = self.sigma * ufl.inner(self.Ep, self.Ep)
        al = q * self.dx(self.domains["Al"])  # Loss in rotor
        steel = q * self.dx(self.domains["Rotor"])  # Loss in only steel
        self._loss_al = fem.form(al, form_compiler_options=self.fp, jit_options=self.jp)
        self._loss_steel = fem.form(steel, form_compiler_options=self.fp, jit_options=self.jp)

    def compute_loss(self, dt: float) -> Tuple[default_scalar_type, default_scalar_type]:
        """
        Compute loss between two time steps of distance dt
        """
        self.dt.value = dt
        al = self.comm.allreduce(fem.assemble_scalar(self._loss_al), op=MPI.SUM)
        steel = self.comm.allreduce(fem.assemble_scalar(self._loss_steel), op=MPI.SUM)
        return (al, steel)

    def _init_torque(self):
        """
        Compute torque induced by magnetic field on the TEAM 30 engine using the surface formulation
        (with Maxwell's stress tensor) or Akkio's method.
        """
        mu_0 = model_parameters["mu_0"]

        dS_air = dS_air = self.dS(surface_map["MidAir"])

        # Create variational form for Electromagnetic torque
        x = ufl.SpatialCoordinate(self.mesh)
        r = ufl.sqrt(x[0] ** 2 + x[1] ** 2)
        dF = 1 / mu_0 * ufl.dot(self.B, x / r) * self.B
        dF -= 1 / mu_0 * 0.5 * ufl.dot(self.B, self.B) * x / r
        torque = self.L * self._restriction * _cross_2D(x, dF)
        torque_surface = (torque("+") + torque("-")) * dS_air
        self._surface_torque = fem.form(
            torque_surface, form_compiler_options=self.fp, jit_options=self.jp
        )

        # Volume formulation of torque (Arkkio's method)
        torque_vol = (
            r
            * self.L
            / (mu_0 * (mesh_parameters["r3"] - mesh_parameters["r2"]))
            * self.Br
            * self.Bphi
        ) * self.dx(self.domains["AirGap"])
        self._volume_torque = fem.form(
            torque_vol, form_compiler_options=self.fp, jit_options=self.jp
        )

    def torque_surface(self) -> float:
        """
        Compute torque using surface integration in air gap and Maxwell's stress tensor
        """
        return self.comm.allreduce(fem.assemble_scalar(self._surface_torque), op=MPI.SUM)

    def torque_volume(self) -> float:
        """
        Compute torque using Arkkio's method, derived on Page 55 of:
        "Analysis of induction motors based on the numerical solution of the magnetic field
        and circuit equations", Antero Arkkio, 1987.
        """
        return self.comm.allreduce(fem.assemble_scalar(self._volume_torque), op=MPI.SUM)


class MagneticField2D:
    def __init__(
        self, Az: fem.Function, form_compiler_options: dict = {}, jit_parameters: dict = {}
    ):
        """
        Class for interpolate the magnetic vector potential to the magnetic flux intensity B=curl(A)

        Parameters
        ==========
        Az
            The magnetic vector potential
        form_compiler_options
            Parameters used in FFCx compilation of this form. Run `ffcx --help` at
            the commandline to see all available options. Takes priority over all
            other parameter values, except for `scalar_type` which is determined by
            DOLFINx.

        jit_parameters
            Parameters used in CFFI JIT compilation of C code generated by FFCx.
            See `python/dolfinx/jit.py` for all available parameters.
            Takes priority over all other parameter values.
        """
        degree = Az.function_space.ufl_element().degree
        mesh = Az.function_space.mesh
        cell = mesh.ufl_cell()

        # Create dolfinx Expression for electromagnetic field B (post processing)
        # Use minimum DG 1 as VTXFile only supports CG/DG>=1
        el_B = basix.ufl.element(
            "DG", cell.cellname(), max(degree - 1, 1), shape=(mesh.geometry.dim,)
        )
        VB = fem.functionspace(mesh, el_B)
        self.B = fem.Function(VB)
        B_2D = ufl.as_vector((Az.dx(1), -Az.dx(0)))
        self.Bexpr = fem.Expression(
            B_2D,
            VB.element.interpolation_points,
            form_compiler_options=form_compiler_options,
            jit_options=jit_parameters,
        )

    def interpolate(self):
        """
        Interpolate magnetic field
        """
        self.B.interpolate(self.Bexpr)


def update_current_density(
    J_0: fem.Function,
    omega: float,
    t: float,
    ct: cpp.mesh.MeshTags_int32,
    currents: Dict[np.int32, Dict[str, float]],
):
    """
    Given a DG-0 scalar field J_0, update it to be alpha*J*cos(omega*t + beta)
    in the domains with copper windings
    """
    J_0.x.array[:] = 0
    for domain, values in currents.items():
        _cells = ct.find(domain)
        J_0.x.array[_cells] = np.full(
            len(_cells),
            model_parameters["J"] * values["alpha"] * np.cos(omega * t + values["beta"]),
        )