.. _MecaSequential:

MecaSequential
==============

This use case aims to validate the methodology POD + ECM for quasistatic elastoviscoplastic computations with Z-set, with a reconstruction of
the dual quantities using the gappy POD, see Section :ref:`Publications`, article 1.

Optional prerequisite: Zmat.

The data needed to execute the tutorials below can be loaded here: :download:`exampleMecaSequential.zip <exampleMecaSequential.zip>`

Description of the physics problem
----------------------------------


Geometry, loading and boundary conditions
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Consider a 1m-wide cube, illustrated in :numref:`label-mesh`. This cube is subjected to

   * a variable thermal loading, in the form of time- and space-dependant scalar field obtained from a previous thermal computation, see :numref:`label-mesh` (right)

   * a centrifugal effect generated by a rotation along the :math:`Z`-axis

   * a time- and space-dependant pressure filed applied on the :math:`X_1` face

   * fixed displacement imposed along the :math:`X`-axis on the :math:`X_0` face, along the :math:`Y`-axis on the :math:`X_0Y_0` edge and along the :math:`Z`-axis on the :math:`X_0Y_1Z_1` vertex


.. figure:: cube.jpg
    :name: label-mesh
    :align: center
    :width: 75%

    (left) Dual mesh of the test case with the accumulated plasticity filed at the end of the simulation, (right) mesh and maximal temperature loading field.


Modeling hypothesis
~~~~~~~~~~~~~~~~~~~

Quasistatic equilibrium equations for the deformable solide in small perturbations.


Material behavior
~~~~~~~~~~~~~~~~~

The material is elastoviscoplastic with nonlinear hardening and a Norton flow.
The constitutive law is specified by the .mat file.


Reduction strategy
------------------

Inputs and outputs
~~~~~~~~~~~~~~~~~~

This is a simple "reproducting use case", where a reduced-order model is constructed to reproduce a high-fidelity solution and check its quality.
This use case contains the modeling complexity of the high-pressure turbine blade cyclic extrapolation, except for the
parallel computation in distributed memory and the fact that the turbine features two material (elastic and elastoviscoplastic).

The number of simulation cycles can be seen as the input, while the outputs are the displacement field ``U`` and the accumulated plasticity field :math:`p` (named ``evrcum``).

Data compression
~~~~~~~~~~~~~~~~

The reduced-order basis is generated by the snapshot-POD method ``CompressData`` from the ``DataCompressors.FusedSnapshotPOD`` module.


Operator compression
~~~~~~~~~~~~~~~~~~~~

The operator compression step is computed using the ECM (Empirical Cubature Method) ``CompressOperator`` from the ``OperatorCompressors.Mechanical`` module.


Algorithm
---------

Results formats and simulation codes
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Snapshots are generated by Z-set and read from the Z-set format.


.. _MecaSequential_verification_qualite_approx:


Quality approximation verification
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


The quality of the *data compression* is evaluated by computed the relative projection error of the high-fidelity solutions on the reduced-order basis:

.. math::
    \dfrac{\left\|u - \sum_{i=1}^N \left( u, \Psi_i^u\right) \Psi_i^u\right\|}{\left\|u\right\|} \le 1 \times 10^{-5},

where :math:`\Psi_i^u` denote the POD modes (or vectors of the reduced-order basis).
The same quantity is considered from :math:`p`.

The quality of the complete procedure is evaluated by compute the relative :math:`\mathcal{l}_2`-norm error between the high-fidelity solutions and the reduced ones:

.. math::
    \dfrac{\left\|u - \sum_{i=1}^N \gamma_i^u \Psi_i^u\right\|}{\left\|u\right\|} \le 1 \times 10^{-5},

with :math:`\gamma_i^u` the coefficients of the reduced solution (or ``reducedCoordinates``) computed by the reduced-order model.


Code - offline stage
--------------------

List of imports required for the *offline* stage of this example:

.. code-block:: python

   from genericROM.IO import ZsetMeshReader as ZMR
   from genericROM.IO import ZsetSolutionReader as ZSR
   from Mordicus.Containers import ProblemData as PD
   from Mordicus.Containers import CollectionProblemData as CPD
   from Mordicus.Containers import Solution as S
   from genericROM.FE import FETools as FT
   from genericROM.DataCompressors import FusedSnapshotPOD as SP
   from genericROM.OperatorCompressors import Mechanical
   from Mordicus.IO import StateIO as SIO
   import numpy as np


Then, filename and dimensions related to the mesh and the solutions have to be declared, readers are initalized, and the mesh is read:

.. code-block:: python

   folder = GetTestDataPath()+"Zset"+os.sep+"MecaSequential"+os.sep

   meshFileName = folder + "cube.geof"
   solutionFileName = folder + "cube.ut"
   meshReader = ZMR.ZsetMeshReader(meshFileName)
   solutionReader = ZSR.ZsetSolutionReader(solutionFileName)
   mesh = meshReader.ReadMesh()
   numberOfNodes = mesh.GetNumberOfNodes()
   numberOfIntegrationPoints = FT.ComputeNumberOfIntegrationPoints(mesh)
   nbeOfComponentsPrimal = 3
   nbeOfComponentsDual = 6


Then, the part of the ECM algorithm depending only on the mesh is carried out:

.. code-block:: python

   operatorPreCompressionData = Mechanical.PreCompressOperator(mesh)

Then, the objects ``Solution`` are declared and populated with data from the precomputed Z-set solutions:

.. code-block:: python

   outputTimeSequence = solutionReader.ReadTimeSequenceFromSolutionFile()
   solutionU = S.Solution("U", nbeOfComponentsPrimal, numberOfNodes, primality = True)
   solutionSigma = S.Solution("sigma", nbeOfComponentsDual, numberOfIntegrationPoints, primality = False)
   solutionEvrcum = S.Solution("evrcum", 1, numberOfIntegrationPoints, primality = False)
   for time in outputTimeSequence:
      solutionU.AddSnapshot(solutionReader.ReadSnapshot("U", time, nbeOfComponentsPrimal, primality=True), time)
      solutionSigma.AddSnapshot(solutionReader.ReadSnapshot("sig", time, nbeOfComponentsDual, primality=False), time)
      solutionEvrcum.AddSnapshot(solutionReader.ReadSnapshotComponent("evrcum", time, primality=False), time)

Then, the objects ``CollectionProblemData`` and ``ProblemData`` are declared, which will enable to agregate the ``Solution`` objects previously constructed in a standard fashion in Mordicus:

.. code-block:: python

   problemData = PD.ProblemData("MecaSequential")
   problemData.AddSolution(solutionU)
   problemData.AddSolution(solutionSigma)
   problemData.AddSolution(solutionEvrcum)
   collectionProblemData = CPD.CollectionProblemData()
   collectionProblemData.AddVariabilityAxis('config', str, description="dummy variability")
   collectionProblemData.DefineQuantity("U", "displacement", "m")
   collectionProblemData.DefineQuantity("sigma", "stress", "Pa")
   collectionProblemData.DefineQuantity("evrcum", "accumulated plasticity", "")
   collectionProblemData.AddProblemData(problemData, config="case-1")

Then, the :math:`L_2 (\Omega)` correlation operator between snapshots is computed (identified by "U"):

.. code-block:: python

   snapshotCorrelationOperator = {"U":FT.ComputeL2ScalarProducMatrix(mesh, 3)}

Then, using the snapshots-POD method, we compute the reduced-order basis for the solutions ``U`` with the :math:`L_2 (\Omega)` correlation operator, and for the solutions :math:`p`
without correlation operator (the default operator is the identity) as a precomputing step for the Gappy-POD reconstruction method on :math:`p`.

.. code-block:: python

   SP.CompressData(collectionProblemData, "U", 1.e-6, snapshotCorrelationOperator["U"])
   SP.CompressData(collectionProblemData, "evrcum", 1.e-6)

Then, we compute the reduced coefficients (or ``reducedCoordinates``) by projecting the high-fidelity snapshots onto the reduced-order basis:

.. code-block:: python

   collectionProblemData.CompressSolutions("U", snapshotCorrelationOperator["U"])

Notice that the two previous steps can be done in one by setting the attribute ``compressSolutions = True`` in the function ``SP.CompressData``.
Then, we verify the quality of the *data compression* on ``U`` as explained in MecaSequential_verification_qualite_approx_ :

.. code-block:: python

   reducedOrderBasisU = collectionProblemData.GetReducedOrderBasis("U")
   CompressedSolutionU = solutionU.GetReducedCoordinates()
   compressionErrors = []
   for t in outputTimeSequence:
      reconstructedCompressedSolution = np.dot(CompressedSolutionU[t], reducedOrderBasisU)
      exactSolution = solutionU.GetSnapshot(t)
      norml2ExSol = np.linalg.norm(exactSolution)
      if norml2ExSol != 0:
         relError = np.linalg.norm(reconstructedCompressedSolution-exactSolution)/norml2ExSol
      else:
         relError = np.linalg.norm(reconstructedCompressedSolution-exactSolution)
         compressionErrors.append(relError)

Then, we carry out the ECM algorithm to determine the reduced quadrature scheme:

.. code-block:: python

   Mechanical.CompressOperator(collectionProblemData, operatorPreCompressionData, mesh, 1.e-5,
   listNameDualVarOutput = ["evrcum"], listNameDualVarGappyIndicesforECM = ["evrcum"])

Finally, at the end of the *offline*, the Modicus data model containing the results of this stage, is saved on disk in order to used them during the *online* stage.

.. code-block:: python

   SIO.SaveState("collectionProblemData", collectionProblemData)
   SIO.SaveState("snapshotCorrelationOperator", snapshotCorrelationOperator)


Code - online stage
-------------------

List of imports required for the *offline* stage of this example:

.. code-block:: python

   from genericROM.IO import ZsetInputReader as ZIR
   from genericROM.IO import ZsetMeshReader as ZMR
   from genericROM.IO import ZsetSolutionReader as ZSR
   from genericROM.IO import ZsetSolutionWriter as ZSW
   from Mordicus.Containers import ProblemData as PD
   from Mordicus.Containers import Solution as S
   from genericROM.FE import FETools as FT
   from genericROM.IO import PXDMFWriter as PW
   from genericROM.OperatorCompressors import Mechanical as Meca
   from Mordicus.IO import StateIO as SIO
   import numpy as np


First, data saved on disk at the end of the *offline* stage is read:

.. code-block:: python

   collectionProblemData = SIO.LoadState("collectionProblemData")
   operatorCompressionDataMechanical = collectionProblemData.GetOperatorCompressionData("U")
   snapshotCorrelationOperator = SIO.LoadState("snapshotCorrelationOperator")
   reducedOrderBases = collectionProblemData.GetReducedOrderBases()

Then, filename and dimensions related to the mesh and the solutions have to be declared and readers are initalized, in the same fashion as the *offline* stage.
We mention here that the physical setting for the *online* stage has been taken identical to the one used in the
*offline* stage (the folder "MecaSequential/"). In meaningful workflow, the physical setting for the *online* stage
would be different.

.. code-block:: python

   folder = GetTestDataPath()+"Zset"+os.sep+"MecaSequential"+os.sep
   inputFileName = folder + "cube.inp"
   inputReader = ZIR.ZsetInputReader(inputFileName)
   meshFileName = folder + "cube.geof"

Then, the mesh is read (which is required when the variability is not parametrized):

.. code-block:: python

   mesh = ZMR.ReadMesh(meshFileName)

Then, an object ``ProblemData`` is defined, which will store the data computed during the *online* stage:

.. code-block:: python

   onlineProblemData = PD.ProblemData("Online")
   onlineProblemData.SetDataFolder(os.path.relpath(folder, folderHandler.scriptFolder))

Then, the temporal sequence and the constitutive law are read from the Z-Set input file. These are "fixed data" for the *online* resolution:

.. code-block:: python

   timeSequence = inputReader.ReadInputTimeSequence()
   constitutiveLawsList = inputReader.ConstructConstitutiveLawsList()
   onlineProblemData.AddConstitutiveLaw(constitutiveLawsList)

Then, the loadings and initial condition are read from the Z-Set input file and are reduced by projecting them onto the reduced-order basis:

.. code-block:: python

   loadingList = inputReader.ConstructLoadingsList()
   onlineProblemData.AddLoading(loadingList)
   for loading in onlineProblemData.GetLoadingsForSolution("U"):
      loading.ReduceLoading(mesh, onlineProblemData, reducedOrderBases, operatorCompressionData)

   initialCondition = inputReader.ConstructInitialCondition()
   onlineProblemData.SetInitialCondition(initialCondition)
   initialCondition.ReduceInitialSnapshot(reducedOrderBases, snapshotCorrelationOperator)

Then, the reduced solution is computed in a nonintrusive fashion using a reduced Newton iterative algorithm for solving the reduced nonlinear system of equations at each time-step:

.. code-block:: python

   onlineCompressedSolution = Meca.ComputeOnline(onlineProblemData, timeSequence, operatorCompressionDataMechanical, 1.e-8)

Then, the reduced coefficients (or ``reducedCoordinates``) for the dual quantified of interest, here :math:`p`, are computed using the *online* part of the Gappy-POD:

.. code-block:: python

   onlineData = onlineProblemData.GetOnlineData("U")
   onlineEvrcumCompressedSolution, errorGappy = Meca.ReconstructDualQuantity("evrcum", operatorCompressionDataMechanical, \
   onlineData, timeSequence = list(onlineCompressedSolution.keys()))

In order to compare the reduced solutions to the high-fidelity reference ones,
``Solution`` objects are created and populated with precomputed Z-set solutions:

.. code-block:: python

   numberOfIntegrationPoints = FT.ComputeNumberOfIntegrationPoints(mesh)
   nbeOfComponentsPrimal = 3
   numberOfNodes = mesh.GetNumberOfNodes()
   solutionFileName = folder + "cube.ut"
   solutionReader = ZSR.ZsetSolutionReader(solutionFileName)
   outputTimeSequence = solutionReader.ReadTimeSequenceFromSolutionFile()

   solutionEvrcumExact = S.Solution("evrcum", 1, numberOfIntegrationPoints, primality = False)
   solutionUExact = S.Solution("U", nbeOfComponentsPrimal, numberOfNodes, primality = True)
   for t in outputTimeSequence:
      evrcum = solutionReader.ReadSnapshotComponent("evrcum", t, primality=False)
      solutionEvrcumExact.AddSnapshot(evrcum, t)
      U = solutionReader.ReadSnapshot("U", t, nbeOfComponentsPrimal, primality=True)
      solutionUExact.AddSnapshot(U, t)

``Solution`` objects corresponding to the reduced solutions are constructed, populated with the reduced coefficients (or ``reducedCoordinates``) computed by the *online* stage, and
reconstruted on the complete mesh:

.. code-block:: python

   solutionEvrcumApprox = S.Solution("evrcum", 1, numberOfIntegrationPoints, primality = False)
   solutionEvrcumApprox.SetReducedCoordinates(onlineEvrcumCompressedSolution)
   solutionEvrcumApprox.UncompressSnapshots(reducedOrderBases["evrcum"])
   solutionUApprox = S.Solution("U", nbeOfComponentsPrimal, numberOfNodes, primality = True)
   solutionUApprox.SetReducedCoordinates(onlineCompressedSolution)
   solutionUApprox.UncompressSnapshots(reducedOrderBases["U"])

Then, we verify the quality of the reduced solutions ``U`` and ``evrcum`` as explained in MecaSequential_verification_qualite_approx_ :

.. code-block:: python

   ROMErrorsU = []
   ROMErrorsEvrcum = []
   for t in outputTimeSequence:
      exactSolution = solutionEvrcumExact.GetSnapshotAtTime(t)
      approxSolution = solutionEvrcumApprox.GetSnapshotAtTime(t)
      norml2ExactSolution = np.linalg.norm(exactSolution)
      if norml2ExactSolution > 1.e-10:
         relError = np.linalg.norm(approxSolution-exactSolution)/norml2ExactSolution
      else:
         relError = np.linalg.norm(approxSolution-exactSolution)
      ROMErrorsEvrcum.append(relError)

      exactSolution = solutionUExact.GetSnapshotAtTime(t)
      approxSolution = solutionUApprox.GetSnapshotAtTime(t)
      norml2ExactSolution = np.linalg.norm(exactSolution)
      if norml2ExactSolution > 1.e-10:
         relError = np.linalg.norm(approxSolution-exactSolution)/norml2ExactSolution
      else:
         relError = np.linalg.norm(approxSolution-exactSolution)
      ROMErrorsU.append(relError)

Finally, reduced predictions for ``U`` and ``evrcum`` are exported in the Z-set format:

.. code-block:: python

   onlineProblemData.AddSolution(solutionUApprox)
   onlineProblemData.AddSolution(solutionEvrcumApprox)
   ZSW.WriteZsetSolution(mesh, meshFileName, "reduced", collectionProblemData, onlineProblemData, "U")



Results
-------

A comparison between the reduced and reference high-fidelity solutions is illustrated in :numref:`label-res`.

.. figure:: res.jpg
    :name: label-res
    :align: center
    :width: 75%

    Comparison between the reduced and reference high-fidelity solutions: (top) on the displacement ``U``, (bottom) on the accumulated plasticity ``p``.