utils.py

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

from typing import Dict, Tuple

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

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, AzV: fem.Function, AnVn: 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
        ==========
        AzV
            The mixed function of the magnetic vector potential Az and the Scalar electric potential V

        AnVn
            The mixed function of the magnetic vector potential Az and the Scalar electric potential V
            from the previous time step

        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 = AzV.function_space.mesh
        self.comm = self.mesh.comm

        # Functions
        Az = AzV[0]
        Azn = AnVn[0]
        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), cells=ct.find(gap_markers[1]))
        self._restriction.interpolate(lambda x: np.zeros(
            x.shape[1], dtype=default_scalar_type), cells=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, AzV: fem.Function,
                 form_compiler_options: dict = {}, jit_parameters: dict = {}):
        """
        Class for interpolate the magnetic vector potential (here as the first part of the mixed function AvZ)
        to the magnetic flux intensity B=curl(A)

        Parameters
        ==========
        AzV
            The mixed function of the magnetic vector potential Az and the Scalar electric potential V

        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 = AzV.function_space.ufl_element().degree()
        mesh = AzV.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,),
                                 gdim=mesh.geometry.dim)
        VB = fem.FunctionSpace(mesh, el_B)
        self.B = fem.Function(VB)
        B_2D = ufl.as_vector((AzV[0].dx(1), -AzV[0].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"]))