Pre/Post deep learning

Some deep learning workflows applied to physics contexts require the projection of fields defined on an unstructured mesh onto a rectilinear grid, and inversely.

Fig. 2 Example of deep learning pre/post

The mesh and solution associated to a previously computed flow field is read (see also Fig. 2, top-left image):

import numpy as np

#################
# PREPROCESSING #
#################

# Read the solution generated by a physical code
from Muscat.IO import XdmfReader as XR
reader = XR.XdmfReader(filename = "PrePostDeepLearning_Input.xmf")
reader.Read()

dom = reader.xdmf.GetDomain(0)
grid = dom.GetGrid(0)

# Read the mesh
uMesh = grid.GetSupport()

# Read the solution field 'U'
indexU = grid.GetPointFieldsNames().index("U")
U = grid.GetPointFields()[indexU][:,0:2]

# Create a rectilinear mesh of size 48*48 excluding the left part of the mesh (x<0)
Nx = 48
Ny = 48
boundingMin, boundingMax = uMesh.ComputeBoundingBox()
Lx = boundingMax[0] - 0.
Ly = boundingMax[1] - boundingMin[1]

from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
unstructuredRectMesh = CreateConstantRectilinearMesh(dimensions=[Nx,Ny], spacing=[Lx/(Nx-1), Ly/(Ny-1)], origin=[0., boundingMin[1]])

# Compute the projection operator from unstructured mesh to structured mesh
from Muscat.FE.FETools import PrepareFEComputation
from Muscat.Containers.MeshFieldOperations import GetFieldTransferOp
from Muscat.FE.Fields.FEField import FEField

space, numberings, _offset, _NGauss = PrepareFEComputation(uMesh, numberOfComponents = 2)
inputFEField = FEField(name="U",mesh=uMesh, space=space, numbering=numberings[0])
methods = ["Interp/Nearest","Nearest/Nearest","Interp/Clamp","Interp/Extrap","Interp/ZeroFill"]
operator, status, entities  = GetFieldTransferOp(inputFEField, unstructuredRectMesh.nodes, method = methods[2], verbose=True)

# Compute the projected field on the structured mesh
projectedU = operator.dot(U)

# Export the structured mesh and projected field in xdmf format (in ASCII)
# To visualize this xdmf file you can use ParaView (downloadable from https://www.paraview.org/)
from Muscat.IO import XdmfWriter as XW
writer = XW.XdmfWriter('PrePostDeepLearning_OutputI.xdmf')
writer.SetHdf5(False)
writer.Open()
writer.Write(unstructuredRectMesh,PointFields=[projectedU], PointFieldsNames=["U"])
writer.Close()


##################
# POSTPROCESSING #
##################

# Compute the projection operator from the structured mesh to the unstructured mesh (inverse projection)
space, numberings, _offset, _NGauss = PrepareFEComputation(unstructuredRectMesh)
inputFEField = FEField(name="U",mesh=unstructuredRectMesh,space=space,numbering=numberings[0])
operator, status, entities = GetFieldTransferOp(inputFEField, uMesh.nodes, method = methods[2], verbose=True)

# Compute the inverse-projected projected field on the unstructured mesh
inverseProjected_ProjectedU = operator.dot(projectedU)

# Export the unstructured mesh and inverse-projected projected field in xdmf format
# To visualize this xdmf file you can use ParaView (downloadable from https://www.paraview.org/)
writer = XW.XdmfWriter('PrePostDeepLearning_OutputII.xdmf')
writer.SetHdf5(False)
writer.Open()
writer.Write(uMesh, PointFields=[inverseProjected_ProjectedU], PointFieldsNames=["U"])
writer.Close()


##################
# CHECK ACCURACY #
##################

# Create a clippedMesh from the unstructured mesh by removing the left part (x<0)
# (with field, inverse-projected projected field and difference between them)
from Muscat.Containers.Filters.FilterObjects import ElementFilter
from Muscat.FE.FETools import ComputePhiAtIntegPoint, ComputeL2ScalarProductMatrix
from Muscat.Containers.MeshInspectionTools import ExtractElementsByElementFilter
from Muscat.Containers.MeshModificationTools import CleanLonelyNodes
from Muscat.Containers.MeshFieldOperations import CopyFieldsFromOriginalMeshToTargetMesh

uMesh.nodeFields['delta'] = U - inverseProjected_ProjectedU

elFilter = ElementFilter(zone = lambda p: (-p[:,0]))
meshClipped = ExtractElementsByElementFilter(uMesh, elFilter)
CleanLonelyNodes(meshClipped)

CopyFieldsFromOriginalMeshToTargetMesh(uMesh, meshClipped)
nbeOfNodes = meshClipped.GetNumberOfNodes()

deltaClippedMesh = meshClipped.nodeFields['delta']

from Muscat.Helpers.TextFormatHelper import TFormat
from Muscat.Helpers.Timer import Timer

print("Compute the L2(Omega) norm of delta by applying numerical quadrature from Lagrange P1")
print("Finite element integration using three different methods")

#### Method 1
print(TFormat.Center(TFormat.InRed("Method 1: ")+TFormat.InBlue("'by hand'")))

timer = Timer("Duration of method 1")
#compute method 1 three times
for i in range(3):
    timer.Start()
    integrationWeights, phiAtIntegPoint = ComputePhiAtIntegPoint(meshClipped)

    vectDeltaAtIntegPoints = np.empty((2,phiAtIntegPoint.shape[0]))
    for i in range(2):
        vectDeltaAtIntegPoints[i] = phiAtIntegPoint.dot(deltaClippedMesh[:,i])

    normDeltaMethod1 = np.sqrt(np.einsum('ij,ij,j', vectDeltaAtIntegPoints, vectDeltaAtIntegPoints, integrationWeights, optimize = True))
    timer.Stop()

print("norm(Delta) =", normDeltaMethod1)
############

#### Method 2
print(TFormat.Center(TFormat.InRed("Method 2: ")+TFormat.InBlue("using the mass matrix")))

timer = Timer("Duration of method 2").Start()

massMatrix = ComputeL2ScalarProductMatrix(meshClipped, 2)

normDeltaMethod2 = np.sqrt(np.dot(deltaClippedMesh.ravel(order='F'), massMatrix.dot(deltaClippedMesh.ravel(order='F'))))

timer.Stop()
print("norm(Delta) =", normDeltaMethod2)
############

#### Method 3
print(TFormat.Center(TFormat.InRed("Method 3: ")+TFormat.InBlue("using the weak form engine")))

from Muscat.FE.Integration import IntegrateGeneral
from Muscat.FE.SymWeakForm import GetField, GetTestField
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP1, ConstantSpaceGlobal
from Muscat.FE.DofNumbering import ComputeDofNumbering

timer = Timer("Duration of method 3").Start()

U = GetField("U",2)
Tt = GetTestField("T",1)

wf = U.T*U*Tt

#meshClipped
numbering = ComputeDofNumbering(meshClipped,LagrangeSpaceP1)
field1 = FEField("U_0",mesh=meshClipped,space=LagrangeSpaceP1,numbering=numbering, data = deltaClippedMesh[:,0])
field1.ConvertDataForNativeTreatment()
field2 = FEField("U_1",mesh=meshClipped,space=LagrangeSpaceP1,numbering=numbering, data = deltaClippedMesh[:,1])
field2.ConvertDataForNativeTreatment()

numbering = ComputeDofNumbering(meshClipped,ConstantSpaceGlobal)
unknownField = FEField("T",mesh=meshClipped,space=ConstantSpaceGlobal,numbering=numbering)

elFilter = ElementFilter()

K, F = IntegrateGeneral(mesh=meshClipped,
                    wform=wf,
                    constants={},
                    fields=[field1,field2],
                    unknownFields=[unknownField],
                    integrationRuleName="LagrangeP1Quadrature",
                    elementFilter=elFilter)

normDeltaMethod3 = np.sqrt(F[0])

timer.Stop()
print("norm(Delta) =", normDeltaMethod2)
############

print(Timer.PrintTimes())