# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE', which is part of this source code package.
#
#
import os
try:
from mpi4py import MPI
except ImportError:# pragma: no cover
print("MPI capabilities not available")
import numpy as np
import collections
from genericROM.OperatorCompressors import ReducedQuadratureProcedure as RQP
#from genericROM.FE import FETools as FT
from genericROM.FE import FETools as MFT
from Muscat.Helpers.TextFormatHelper import TFormat
import time
import sys
from genericROM.DataCompressors import FusedSnapshotPOD as SP
from genericROM.Containers.OperatorCompressionData import OperatorPreCompressionDataThermal as OPCDT
from genericROM.Containers.OperatorCompressionData import OperatorCompressionDataThermal as OCDT
"""
#### Description of 'operatorPreCompressionData': dict of keys
#PreCompress generated:
- listOfTags
- integrationWeightsRadiation
- nodeToIntegThermalRadiationOperator
- integrationWeights
- nodeToIntegThermalOperator
- nodeToIntegGradThermalOperator
#### Description of 'operatorCompressionData': dict of keys
#Compress generated:
- reducedIntegrationPointsRadiation
- reducedIntegrationWeightsRadiation
- phiRadAtIntegPoints
- phiphiRadAtIntegPoints
- reducedIntegrationPointsCapacity
- reducedIntegrationWeightsCapacity
- phiCapacityAtReducedIntegPoints
- phiphiCapacityAtReducedIntegPoints
- reducedListOTagsCapacity
- reducedIntegrationPointsConductivity
- reducedIntegrationWeightsConductivity
- phiConductivityAtReducedIntegPoints
- gradgradAtReducedIntegPoints
- reducedListOTagsConductivity
"""
[docs]
def ComputeOnline(onlineProblemData, timeSequence, operatorCompressionData, tolerance, callback = None):
"""
Compute the online stage using the method POD and ECM for a thermal problem
The parameters must have been initialized in onlineProblemData
"""
currentFolder = os.getcwd()
folder = currentFolder + os.sep + onlineProblemData.GetDataFolder()
os.chdir(folder)
initialCondition = onlineProblemData.GetInitialCondition()
onlineCompressedSolution = collections.OrderedDict()
onlineCompressedSolution[timeSequence[0]] = initialCondition.GetReducedInitialSnapshot("T")
#nModes = onlineCompressedSolution[timeSequence[0]].shape[0]
#keysConstitutiveLaws = set(onlineProblemData.GetConstitutiveLaws().keys())
constitutiveLawSets = onlineProblemData.GetSetsOfConstitutiveOfType("thermal")
####### Capacity ##############
reducedListOTagsCapacity = operatorCompressionData.GetReducedListOTagsCapacity()
IndicesOfIntegPointsPerMaterialCapacity = MFT.ComputeIndicesOfIntegPointsPerMaterial(reducedListOTagsCapacity, constitutiveLawSets)
####### Conductivity ##########
reducedListOTagsConductivity = operatorCompressionData.GetReducedListOTagsConductivity()
IndicesOfIntegPointsPerMaterialConductivity = MFT.ComputeIndicesOfIntegPointsPerMaterial(reducedListOTagsConductivity, constitutiveLawSets)
for timeStep in range(1, len(timeSequence)):
previousTime = timeSequence[timeStep-1]
time = timeSequence[timeStep]
dtime = time - previousTime
if callback == None:
print("time =", time); sys.stdout.flush()
else:
callback.CurrentTime(timeStep, time)
reducedExternalForcesTemp = ComputeReducedExternalForces(onlineProblemData, time)
reducedExternalForces = np.zeros(reducedExternalForcesTemp.shape)
MPI.COMM_WORLD.Allreduce([reducedExternalForcesTemp, MPI.DOUBLE], [reducedExternalForces, MPI.DOUBLE])
normExt = np.linalg.norm(reducedExternalForces)
before = onlineCompressedSolution[previousTime]# np.copy in niROM
after = before# np.copy in niROM
onlineCompressedSolution[time] = np.copy(before)
count = 0
normRes = 1
while normRes > tolerance:
reducedInternalForcesTemp, reducedTangentMatrixTemp = ComputeReducedInternalForcesAndTangentMatrix(onlineProblemData, operatorCompressionData, IndicesOfIntegPointsPerMaterialCapacity, IndicesOfIntegPointsPerMaterialConductivity, time, dtime, before, after)
reducedInternalForces = np.zeros(reducedInternalForcesTemp.shape)
MPI.COMM_WORLD.Allreduce([reducedInternalForcesTemp, MPI.DOUBLE], [reducedInternalForces, MPI.DOUBLE])
if normExt>0:
normRes = np.linalg.norm(reducedExternalForces-reducedInternalForces)/normExt
else:
normRes = np.linalg.norm(reducedExternalForces-reducedInternalForces)# pragma: no cover
if callback == None:
print("normRes =", normRes)
else:
callback.CurrentNormRes(normRes)
reducedTangentMatrix = np.zeros(reducedTangentMatrixTemp.shape)
MPI.COMM_WORLD.Allreduce([reducedTangentMatrixTemp, MPI.DOUBLE], [reducedTangentMatrix, MPI.DOUBLE])
sol = np.linalg.solve(reducedTangentMatrix, reducedExternalForces-reducedInternalForces)
onlineCompressedSolution[time] += sol
after = onlineCompressedSolution[time]# np.copy in niROM
count += 1
if count == 20:
raise RuntimeError("problem could not converge after 20 iterations") # pragma: no cover
if callback == None:
print("=== Newton iterations:", count)
else:
callback.CurrentNewtonIterations(count)
os.chdir(currentFolder)
return onlineCompressedSolution
[docs]
def ComputeReducedExternalForces(problemData, time):
reducedExternalForces = 0.
for loading in problemData.GetLoadings().values():
if loading.GetSolutionName() == "T":
reducedExternalForces += loading.ComputeContributionToReducedExternalForces(time)
return reducedExternalForces
[docs]
def ComputeReducedInternalForcesAndTangentMatrix(onlineProblemData, opCompDat, IndicesOfIntegPointsPerMaterialCapacity, IndicesOfIntegPointsPerMaterialConductivity, time, dtime, before, after):
constLaws = onlineProblemData.GetConstitutiveLaws()
reducedTangentMatrix = 0.
reducedInternalForces = 0.
#### capacity terms ############################
reducedIntegrationWeightsCapacity = opCompDat.GetReducedIntegrationWeightsCapacity()
reducedPhiAtCapacityReducedIntegPoint = opCompDat.GetReducedPhiAtCapacityReducedIntegPoint()
reducedPhiPhiAtCapacityReducedIntegPoint = opCompDat.GetReducedPhiPhiAtCapacityReducedIntegPoint()
TredIntegBefore = np.dot(before, reducedPhiAtCapacityReducedIntegPoint)#before on reduced integration points
TredIntegAfter = np.dot(after, reducedPhiAtCapacityReducedIntegPoint)#after on reduced integration points
res = 0.
mat = 0.
deltaEnergyInteg = np.empty(len(reducedIntegrationWeightsCapacity))
capacityInteg = np.empty(len(reducedIntegrationWeightsCapacity))
for tag, intPoints in IndicesOfIntegPointsPerMaterialCapacity.items():
lawTag = ('thermal', tag)
deltaEnergyInteg[intPoints] = constLaws[lawTag].ComputeInternalEnergyVectorized(TredIntegAfter[intPoints])-constLaws[lawTag].ComputeInternalEnergyVectorized(TredIntegBefore[intPoints])
capacityInteg[intPoints] = constLaws[lawTag].ComputeCapacity(TredIntegAfter[intPoints])
res += np.dot(reducedPhiAtCapacityReducedIntegPoint, np.multiply(reducedIntegrationWeightsCapacity, deltaEnergyInteg))
mat += np.dot(reducedPhiPhiAtCapacityReducedIntegPoint, np.multiply(reducedIntegrationWeightsCapacity, capacityInteg))
reducedInternalForces += res/dtime
reducedTangentMatrix += mat/dtime
################################################
#### conductivity terms ########################
reducedIntegrationWeightsConductivity = opCompDat.GetReducedIntegrationWeightsConductivity()
reducedPhiAtConductivityReducedIntegPoint = opCompDat.GetReducedPhiAtConductivityReducedIntegPoint()
reducedGradGradAtConductivityReducedIntegPoint = opCompDat.GetReducedGradGradAtConductivityReducedIntegPoint()
TredIntegAfter = np.dot(after, reducedPhiAtConductivityReducedIntegPoint)#after on reduced integration points
conductivityInteg = np.empty(len(reducedIntegrationWeightsConductivity))
for tag, intPoints in IndicesOfIntegPointsPerMaterialConductivity.items():
lawTag = ('thermal', tag)
conductivityInteg[intPoints] = constLaws[lawTag].ComputeConductivity(TredIntegAfter[intPoints])
reducedInternalForces += np.dot(np.multiply(np.dot(after, reducedGradGradAtConductivityReducedIntegPoint),reducedIntegrationWeightsConductivity), conductivityInteg)
reducedTangentMatrix += np.dot(np.multiply(reducedGradGradAtConductivityReducedIntegPoint,reducedIntegrationWeightsConductivity), conductivityInteg)
################################################
#### radiation terms ###########################
bcRad = onlineProblemData.GetLoadingsOfType('radiation')
assert len(bcRad) < 2 #only one radiation
if len(bcRad) == 1 and opCompDat.GetReducedIntegrationWeightsRadiation() is not None:
reducedIntegrationWeightsRadiation = opCompDat.GetReducedIntegrationWeightsRadiation()
reducedPhiAtRadReducedIntegPoint = opCompDat.GetReducedPhiAtRadReducedIntegPoint()
reducedPhiPhiAtRadReducedIntegPoint = opCompDat.GetReducedPhiPhiAtRadReducedIntegPoint()
TredInteg = np.dot(after, reducedPhiAtRadReducedIntegPoint)#T on reduced integration points
T4redInteg = np.power(TredInteg,4)#T^4 on reduced integration points
T4BoundRedRad = np.einsum('kl,l->kl', reducedPhiAtRadReducedIntegPoint, T4redInteg, optimize = True)
# cf https://stackoverflow.com/questions/19388152/numpy-element-wise-multiplication-of-an-array-and-a-vector
reducedInternalForces += bcRad[0].stefanBoltzmannConstant*np.dot(T4BoundRedRad, reducedIntegrationWeightsRadiation)
T3redInteg = np.power(TredInteg,3)#T^3 on reduced integration points
T3BoundRedRadMat = np.einsum('klm,m->klm', reducedPhiPhiAtRadReducedIntegPoint, T3redInteg, optimize = True)
reducedTangentMatrix += 4*bcRad[0].stefanBoltzmannConstant*np.dot(T3BoundRedRadMat, reducedIntegrationWeightsRadiation)
################################################
#### convection terms ##########################
bcConv = onlineProblemData.GetLoadingsOfType('convection_heat_flux')
for bc in bcConv:
convTensor = bc.reducedPhiTPhiT
h,_ = bc.GetCoefficientsAtTime(time)
reducedTangentMatrix += h*convTensor
reducedInternalForces += h*np.dot(convTensor,after)
################################################
return reducedInternalForces, reducedTangentMatrix
[docs]
def PreCompressOperator(mesh, radiationSet = None):
listOfTags = MFT.ComputeIntegrationPointsTags(mesh, mesh.GetDimensionality())
if radiationSet != None:
integrationWeightsRadiation, phiAtIntegPointRadiationSet = MFT.ComputePhiAtIntegPoint(mesh, [radiationSet], -1)
else:
integrationWeightsRadiation, phiAtIntegPointRadiationSet = None, None
"""integrationWeightsRadiation0, phiAtIntegPointRadiationSet0 = MFT.ComputePhiAtIntegPoint(mesh, [radiationSet], -1)
integrationWeightsRadiation, phiAtIntegPointRadiationSet = FT.ComputeNodeToIntegThermalOperatorRadiation(mesh, radiationSet)
phiAtIntegPointRadiationSet = phiAtIntegPointRadiationSet.T
print("rel diff integrationWeightsRadiation =", np.linalg.norm(integrationWeightsRadiation-integrationWeightsRadiation0)/np.linalg.norm(integrationWeightsRadiation0))
print("rel diff phiAtIntegPointRadiationSet0 =", np.linalg.norm(phiAtIntegPointRadiationSet0.todense()-phiAtIntegPointRadiationSet.todense())/np.linalg.norm(phiAtIntegPointRadiationSet0.todense()))
1./0."""
integrationWeights, phiAtIntegPoint = MFT.ComputePhiAtIntegPoint(mesh)
_, gradPhiAtIntegPoint = MFT.ComputeGradPhiAtIntegPoint(mesh)
operatorPreCompressionData = OPCDT.OperatorPreCompressionDataThermal("T", \
gradPhiAtIntegPoint, phiAtIntegPoint, integrationWeights, listOfTags, integrationWeightsRadiation, phiAtIntegPointRadiationSet)
return operatorPreCompressionData
[docs]
def CompressOperator(
collectionProblemData, operatorPreCompressionData, mesh, tolerance, toleranceCompressSnapshotsForRedQuad = 0.
):
"""
Operator Compression for the POD and ECM for a thermal problem
Parameters
----------
collectionProblemData : CollectionProblemData
definition of the training data in a CollectionProblemData object
mesh : MeshBase
mesh
tolerance : float
tolerance for Empirical Cubature Method
Returns
-------
operatorCompressionOutputData : (regressor, scaler, scaler)
(fitted regressor, fitted scaler on the coefficients, fitted scaler on the parameters)
"""
print("CompressOperator starting...")
if toleranceCompressSnapshotsForRedQuad > 0:
collectionProblemData.DefineQuantity("Cond")
collectionProblemData.DefineQuantity("IntEnergy")
collectionProblemData.DefineQuantity("T4Rad")
reducedOrderBasis = collectionProblemData.GetReducedOrderBasis("T")
snapshots = collectionProblemData.GetSnapshots("T")
operatorCompressionData = collectionProblemData.GetOperatorCompressionData("T")
if operatorCompressionData == None:
operatorCompressionData = OCDT.OperatorCompressionDataThermal("T")
operatorCompressionData.SetOperatorPreCompressionData(operatorPreCompressionData)
collectionProblemData.AddOperatorCompressionData(operatorCompressionData)
listOfTags = operatorCompressionData.GetListOfTags()
ProblemDataIndicesOfIntegPointsPerMaterial = {}
for key, problemData in collectionProblemData.GetProblemDatas().items():
#laws = problemData.GetConstitutiveLaws()
#IndicesOfIntegPointsPerMaterial = MFT.ComputeIndicesOfIntegPointsPerMaterial(listOfTags, set(laws.keys()))
constitutiveLawSets = problemData.GetSetsOfConstitutiveOfType("thermal")
IndicesOfIntegPointsPerMaterial = MFT.ComputeIndicesOfIntegPointsPerMaterial(listOfTags, constitutiveLawSets)
ProblemDataIndicesOfIntegPointsPerMaterial[key] = IndicesOfIntegPointsPerMaterial
integrationWeightsRadiation = operatorCompressionData.GetIntegrationWeightsRadiation()
phiAtIntegPointRadiationSet = operatorCompressionData.GetPhiAtIntegPointRadiationSet()
integrationWeights = operatorCompressionData.GetIntegrationWeights()
phiAtIntegPoint = operatorCompressionData.GetPhiAtIntegPoint()
gradPhiAtIntegPoint = operatorCompressionData.GetGradPhiAtIntegPoint()
#### Radiation integrator #####
if integrationWeightsRadiation is not None:
start = time.time()
print(TFormat.InGreen("Starting computation reduced quadrature for the Radiation term"))
#reduced integrators
reducedPhiRadAtIntegPoint = phiAtIntegPointRadiationSet.dot(reducedOrderBasis.T).T
reducedPhiPhiRadAtIntegPoint = np.einsum('kl,ml->kml', reducedPhiRadAtIntegPoint, reducedPhiRadAtIntegPoint, optimize = True)
T4BoundRad = ComputeT4OnRadiationSet(collectionProblemData, phiAtIntegPointRadiationSet, reducedPhiRadAtIntegPoint, snapshots, toleranceCompressSnapshotsForRedQuad)
#rank = MPI.COMM_WORLD.Get_rank()
#T4BoundRadnew = np.load("T4BoundRad"+str(rank)+".npy")
#T4BoundRad2 = np.dot(T4BoundRad, integrationWeightsRadiation)
#print("np.linalg.norm(T4BoundRad2) = ", np.linalg.norm(T4BoundRad2))
#print("np.linalg.norm(T4BoundRadnew) = ", np.linalg.norm(T4BoundRadnew))
#print("rel diff T4BoundRadnew =", np.linalg.norm(T4BoundRadnew - T4BoundRad2)/np.linalg.norm(T4BoundRad2))
reducedIntegrationPointsInitSet = operatorCompressionData.GetReducedIntegrationPointsRadiation()
reducedIntegrationPointsRadiation, reducedIntegrationWeightsRadiation = RQP.ComputeReducedIntegrationScheme(integrationWeightsRadiation, T4BoundRad, tolerance, imposedIndices = None, reducedIntegrationPointsInitSet = reducedIntegrationPointsInitSet)
#reducedIntegrationPointsRadiation = np.arange(len(integrationWeightsRadiation))
#reducedIntegrationWeightsRadiation = integrationWeightsRadiation
operatorCompressionData.SetReducedIntegrationPointsRadiation(reducedIntegrationPointsRadiation)
operatorCompressionData.SetReducedIntegrationWeightsRadiation(reducedIntegrationWeightsRadiation)
#hyperreduced integrators
operatorCompressionData.SetReducedPhiAtRadReducedIntegPoint(reducedPhiRadAtIntegPoint[:,reducedIntegrationPointsRadiation])
operatorCompressionData.SetReducedPhiPhiAtRadReducedIntegPoint(reducedPhiPhiRadAtIntegPoint[:,:,reducedIntegrationPointsRadiation])
print(TFormat.InGreen("Duration computation reduced quadrature for the Radiation term: "+str(time.time()-start)))
################################
tempAtIntegPoints = phiAtIntegPoint.dot(snapshots.T).T
#reduced integrators
reducedPhiAtIntegPoint = phiAtIntegPoint.dot(reducedOrderBasis.T).T
reducedPhiPhiAtIntegPoint = np.einsum('kl,ml->kml', reducedPhiAtIntegPoint, reducedPhiAtIntegPoint, optimize = True)
reducedGradAtIntegPoint = np.empty((reducedOrderBasis.shape[0], len(integrationWeights), mesh.GetDimensionality()))
for i in range(mesh.GetDimensionality()):
reducedGradAtIntegPoint[:,:,i] = gradPhiAtIntegPoint[i].dot(reducedOrderBasis.T).T
reducedGradGradAtIntegPoint = np.einsum('klt,mlt->kml', reducedGradAtIntegPoint , reducedGradAtIntegPoint , optimize = True)
#### Capacity integrator #########
start = time.time()
print(TFormat.InGreen("Starting computation reduced quadrature for the Capacity term"))
internalEnergy = ComputeInternalEnergy(collectionProblemData, tempAtIntegPoints, reducedPhiAtIntegPoint, ProblemDataIndicesOfIntegPointsPerMaterial, snapshots, toleranceCompressSnapshotsForRedQuad)
#internalEnergynew = np.load("internalEnergy"+str(rank)+".npy")
#internalEnergy2 = np.dot(internalEnergy, integrationWeights)
#print("rel diff internalEnergy =", np.linalg.norm(internalEnergynew - internalEnergy2)/np.linalg.norm(internalEnergy2))
#1./0.
reducedIntegrationPointsInitSet = operatorCompressionData.GetReducedIntegrationPointsCapacity()
reducedIntegrationPointsCapacity, reducedIntegrationWeightsCapacity = RQP.ComputeReducedIntegrationScheme(integrationWeights, internalEnergy, tolerance, imposedIndices = None, reducedIntegrationPointsInitSet = reducedIntegrationPointsInitSet)
#reducedIntegrationPointsCapacity = np.arange(len(integrationWeights))
#reducedIntegrationWeightsCapacity = integrationWeights
operatorCompressionData.SetReducedIntegrationPointsCapacity(reducedIntegrationPointsCapacity)
operatorCompressionData.SetReducedIntegrationWeightsCapacity(reducedIntegrationWeightsCapacity)
#hyperreduced integrators
operatorCompressionData.SetReducedPhiAtCapacityReducedIntegPoint(reducedPhiAtIntegPoint[:,reducedIntegrationPointsCapacity])
operatorCompressionData.SetReducedPhiPhiAtCapacityReducedIntegPoint(reducedPhiPhiAtIntegPoint[:,:,reducedIntegrationPointsCapacity])
reducedListOTagsCapacity = PrecomputeReducedMaterialVariable(mesh, listOfTags, reducedIntegrationPointsCapacity)
operatorCompressionData.SetReducedListOTagsCapacity(reducedListOTagsCapacity)
print(TFormat.InGreen("Duration computation reduced quadrature for the Capacity term: "+str(time.time()-start)))
##################################
#### Conductivity integrator #####
start = time.time()
print(TFormat.InGreen("Starting computation reduced quadrature for the Conductivity term"))
gradGradT = ComputeGradGradT(collectionProblemData, tempAtIntegPoints, reducedPhiAtIntegPoint, reducedGradGradAtIntegPoint, ProblemDataIndicesOfIntegPointsPerMaterial, snapshots, toleranceCompressSnapshotsForRedQuad)
#gradGradTnew = np.load("gradGradT"+str(rank)+".npy")
#gradGradT2 = np.dot(gradGradT, integrationWeights)
#print("rel diff gradGradT =", np.linalg.norm(gradGradTnew - gradGradT2)/np.linalg.norm(gradGradT2))
#1./0.
reducedIntegrationPointsInitSet = operatorCompressionData.GetReducedIntegrationPointsConductivity()
reducedIntegrationPointsConductivity, reducedIntegrationWeightsConductivity = RQP.ComputeReducedIntegrationScheme(integrationWeights, gradGradT, tolerance, imposedIndices = None, reducedIntegrationPointsInitSet = reducedIntegrationPointsInitSet)
#reducedIntegrationPointsConductivity = np.arange(len(integrationWeights))
#reducedIntegrationWeightsConductivity = integrationWeights
operatorCompressionData.SetReducedIntegrationPointsConductivity(reducedIntegrationPointsConductivity)
operatorCompressionData.SetReducedIntegrationWeightsConductivity(reducedIntegrationWeightsConductivity)
#hyperreduced integrators
operatorCompressionData.SetReducedPhiAtConductivityReducedIntegPoint(reducedPhiAtIntegPoint[:,reducedIntegrationPointsConductivity])
operatorCompressionData.SetReducedGradGradAtConductivityReducedIntegPoint(reducedGradGradAtIntegPoint[:,:,reducedIntegrationPointsConductivity])
reducedListOTagsConductivity = PrecomputeReducedMaterialVariable(mesh, listOfTags, reducedIntegrationPointsConductivity)
operatorCompressionData.SetReducedListOTagsConductivity(reducedListOTagsConductivity)
print(TFormat.InGreen("Duration computation reduced quadrature for the Conductivity term: "+str(time.time()-start)))
##################################
[docs]
def ComputeGradGradT(collectionProblemData, tempAtIntegPoints, phi, gradgrad, ProblemDataIndicesOfIntegPointsPerMaterial, snapshots, toleranceCompressSnapshotsForRedQuad):
"""
computes lambda(T_{i+1})(x_k) sum_l hat{T}_{i+1}_l grad{Psi_j} cdot grad{Psi_l} cdot(x_k)
"""
#first update Cond reducedOrderBasisCond
TNumberOfSnapshots = tempAtIntegPoints.shape[0]
numberOfModes = collectionProblemData.GetReducedOrderBasisNumberOfModes("T")
reducedCoordinates = np.empty((TNumberOfSnapshots,numberOfModes))
arrayOfCond = np.empty(tempAtIntegPoints.shape)
count = 0
for key, problemData in collectionProblemData.GetProblemDatas().items():
laws = problemData.GetConstitutiveLaws()
localNbeSnapshots = problemData.GetSolution("T").GetNumberOfSnapshots()
IndicesOfIntegPointsPerMaterial = ProblemDataIndicesOfIntegPointsPerMaterial[key]
reducedCoordinates[count:count+localNbeSnapshots,:] = np.array(
problemData.GetSolution("T").GetReducedCoordinatesAtTimes(problemData.GetSolution("T").GetTimeSequenceFromSnapshots()),
dtype = object)
for key2, value in IndicesOfIntegPointsPerMaterial.items():
numberOfMaterialIntegPoints = len(IndicesOfIntegPointsPerMaterial[key2])
arrayOfCond[count:count+localNbeSnapshots,IndicesOfIntegPointsPerMaterial[key2]] = laws[('thermal', key2)].ComputeConductivity(tempAtIntegPoints[count:count+localNbeSnapshots,IndicesOfIntegPointsPerMaterial[key2]].ravel()).reshape((localNbeSnapshots,numberOfMaterialIntegPoints))
count += localNbeSnapshots
numberOfIntegrationPoints = phi.shape[1]
if toleranceCompressSnapshotsForRedQuad > 0.:
SP.CompressData(collectionProblemData, "Cond", toleranceCompressSnapshotsForRedQuad, snapshots = arrayOfCond)
reducedOrderBasisCond = collectionProblemData.GetReducedOrderBasis("Cond")
numberOfModesCond = collectionProblemData.GetReducedOrderBasisNumberOfModes("Cond")
basisProjection = np.dot(reducedOrderBasisCond, phi.T)
gradBasisProjection = np.dot(basisProjection,gradgrad)
return np.einsum('km,klm->klm', reducedOrderBasisCond, gradBasisProjection, optimize = True).reshape(numberOfModesCond*numberOfModes,numberOfIntegrationPoints)
else:
gradgradhatT = np.dot(reducedCoordinates, gradgrad)
return np.einsum('km,klm->klm', arrayOfCond, gradgradhatT, optimize = True).reshape(TNumberOfSnapshots*numberOfModes,numberOfIntegrationPoints)
[docs]
def ComputeInternalEnergy(collectionProblemData, tempAtIntegPoints, phi, ProblemDataIndicesOfIntegPointsPerMaterial, snapshots, toleranceCompressSnapshotsForRedQuad):
"""
computes U(T_{i})(x_k) * Psi_j(x_k)
"""
#first update IntEnergy reducedOrderBasisIntEnergy
arrayOfIntEnergy = np.empty(tempAtIntegPoints.shape)
count = 0
for key, problemData in collectionProblemData.GetProblemDatas().items():
laws = problemData.GetConstitutiveLaws()
localNbeSnapshots = problemData.GetSolution("T").GetNumberOfSnapshots()
IndicesOfIntegPointsPerMaterial = ProblemDataIndicesOfIntegPointsPerMaterial[key]
for key2, value in IndicesOfIntegPointsPerMaterial.items():
numberOfMaterialIntegPoints = len(IndicesOfIntegPointsPerMaterial[key2])
arrayOfIntEnergy[count:count+localNbeSnapshots,IndicesOfIntegPointsPerMaterial[key2]] = laws[('thermal', key2)].ComputeInternalEnergyVectorized(tempAtIntegPoints[count:count+localNbeSnapshots,IndicesOfIntegPointsPerMaterial[key2]].ravel()).reshape((localNbeSnapshots,numberOfMaterialIntegPoints))
count += localNbeSnapshots
numberOfIntegrationPoints = phi.shape[1]
numberOfModes = collectionProblemData.GetReducedOrderBasisNumberOfModes("T")
if toleranceCompressSnapshotsForRedQuad > 0.:
SP.CompressData(collectionProblemData, "IntEnergy", toleranceCompressSnapshotsForRedQuad, snapshots = arrayOfIntEnergy)
reducedOrderBasisIntEnergy = collectionProblemData.GetReducedOrderBasis("IntEnergy")
numberOfModesIntEnergy = collectionProblemData.GetReducedOrderBasisNumberOfModes("IntEnergy")
return np.einsum('kl,ml->kml', phi, reducedOrderBasisIntEnergy, optimize = True).reshape(numberOfModes*numberOfModesIntEnergy,numberOfIntegrationPoints)
else:
TNumberOfSnapshots = tempAtIntegPoints.shape[0]
return np.einsum('kl,ml->kml', phi, arrayOfIntEnergy, optimize = True).reshape(numberOfModes*TNumberOfSnapshots,numberOfIntegrationPoints)
[docs]
def ComputeT4OnRadiationSet(collectionProblemData, phiAtIntegPointRadiationSet, phiRad, snapshots, toleranceCompressSnapshotsForRedQuad):
"""
computes T_i^4(x_k) Psi_j(x_k)
"""
#first update T4rad reducedOrderBasisT4Rad
arrayOfT4RadSnapshots = np.power(phiAtIntegPointRadiationSet.dot(snapshots.T), 4).T
numberOfModes = collectionProblemData.GetReducedOrderBasisNumberOfModes("T")
numberOfIntegrationPoints = phiRad.shape[1]
if toleranceCompressSnapshotsForRedQuad > 0.:
SP.CompressData(collectionProblemData, "T4Rad", toleranceCompressSnapshotsForRedQuad, snapshots = arrayOfT4RadSnapshots)
reducedOrderBasisT4Rad = collectionProblemData.GetReducedOrderBasis("T4Rad")
numberOfModesT4Rad = collectionProblemData.GetReducedOrderBasisNumberOfModes("T4Rad")
return np.einsum('kl,ml->kml', phiRad, reducedOrderBasisT4Rad, optimize = True).reshape(numberOfModes*numberOfModesT4Rad,numberOfIntegrationPoints)
else:
TNumberOfSnapshots = arrayOfT4RadSnapshots.shape[0]
return np.einsum('kl,ml->kml', phiRad, arrayOfT4RadSnapshots, optimize = True).reshape(numberOfModes*TNumberOfSnapshots,numberOfIntegrationPoints)
[docs]
def PrecomputeReducedMaterialVariable(mesh, listOfTags, reducedIntegrationPoints):
"""
1.
"""
reducedListOTags = [list(listOfTags[intPoint]) for intPoint in reducedIntegrationPoints]
for i in range(len(reducedIntegrationPoints)):
reducedListOTags[i].append("ALLELEMENT")
return reducedListOTags
if __name__ == "__main__":# pragma: no cover
from genericROM import RunTestFile
RunTestFile(__file__)