# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE.txt', which is part of this source code package.
#
import time
import numpy as np
from scipy.sparse import coo_matrix
from Muscat.Types import MuscatIndex, MuscatFloat
from Muscat.Helpers.Decorators import froze_it
from Muscat.Helpers.Logger import Debug, Error
import Muscat.MeshContainers.ElementsDescription as ED
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
from Muscat.MeshContainers.Filters.FilterOperators import FrozenFilter
from Muscat.MeshContainers.Filters.FilterTools import GetListOfPartialElementFilter, GetFrozenFilter
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo
from Muscat.FE.Fields.FEField import FEField
from Muscat.FE.Fields.IPField import IPField
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.FETools import PrepareFEComputation
from Muscat.FE.SymWeakForm import GetField, GetTestField, GetNormal
import Muscat.FE.WeakForms.NumericalWeakForm as WeakForm
from Muscat.FE.IntegrationRules import LagrangeIsoParamQuadrature, IntegrationRulesFactory, GetIntegrationRuleByName
UseKokkos = True
UseCpp = True
UseMultiThread = True
MultiThreadThreshold = 1000
[docs]
def IntegrateGeneral(mesh, wform, constants, fields, unknownFields=None, testFields=None,
integrationRuleName=None, onlyEvaluation=False, elementFilter=None,
userIntegrator=None, integrationRule=None, nbThreads = None ):
"""Integration of a weak formulation
For more information about the argument please refer to IntegrationClass
Returns
-------
K : coo_matrix
the assembled matrix
rhs : ndarray
Array with the values of the right hand side term
"""
if userIntegrator is None:
intClass = GetIntegrator()
else:
intClass = userIntegrator
intClass.SetMesh(mesh)
intClass.SetOnlyEvaluation(onlyEvaluation)
intClass.SetElementFilter(elementFilter)
intClass.SetIntegrationRule(integrationRuleName=integrationRuleName, integrationRule=integrationRule)
intClass.SetConstants(constants)
intClass.SetUnknownFields(unknownFields)
intClass.SetTestFields(testFields)
intClass.SetExtraFields(fields)
intClass.SetWeakForm(wform)
intClass.PreStartCheck()
intClass.Allocate()
intClass.Compute(nbThreads=nbThreads)
return intClass.GetK(), intClass.GetRhs()
[docs]
def GetIntegrator():
import Muscat.FE.Integration as BTFEI
if BTFEI.UseKokkos:
from Muscat.FE.Integrators.NativeKokkosIntegration import IntegrationKokkos
return IntegrationKokkos()
elif BTFEI.UseCpp:
from Muscat.FE.Integrators.NativeIntegration import IntegratorCpp
return IntegratorCpp()
else:
from Muscat.FE.Integrators.PythonIntegration import IntegratorPython
return IntegratorPython()
[docs]
def IntegrateGeneralMonoThread(mesh, wform, constants, fields, unknownFields=None, testFields=None, integrationRuleName=None, onlyEvaluation=False, elementFilter=None, userIntegrator=None, integrationRule=None):
"""Integration of a weak formulation using only one thread
For more information about the argument please refer to IntegrationClass
Returns
-------
K : coo_matrix
the assembled matrix
rhs : ndarray
Array with the values of the right hand side term
"""
return IntegrateGeneral(mesh, wform, constants, fields, unknownFields=unknownFields, testFields=testFields,
integrationRuleName=integrationRuleName, onlyEvaluation=onlyEvaluation, elementFilter=elementFilter,
integrationRule=integrationRule, nbThreads = 1 )
##################### CheckIntegrity #####################
[docs]
def CheckIntegrityNormalFlux(GUI: bool = False):
from Muscat.MeshTools.MeshCreationTools import CreateMeshOf
points = [[0, 0, 0], [1, 0, 1], [0, 1, 1]]
mesh = CreateMeshOf(points, [[0, 1, 2]], ED.Triangle_3)
mesh.ConvertDataForNativeTreatment()
sdim = 3
dofs = ["u_" + str(i) for i in range(sdim)]
ut = GetTestField("u", sdim,sdim)
p = GetField("p", 1,sdim)
normal = GetNormal(3)
wformflux = p*normal.T*ut
constants = {"alpha": 1.0}
pressField = FEField("p", mesh)
pressField.Allocate(1)
print(mesh)
elemfilt = ElementFilter() # all elements
unknownFields = [FEField(dofs[n], mesh=mesh, space=pressField.space, numbering=pressField.numbering) for n in range(len(dofs))]
_matrix, F = IntegrateGeneral(mesh=mesh, wform=wformflux, constants=constants, fields=[pressField],
unknownFields=unknownFields,
elementFilter=elemfilt)
np.testing.assert_allclose(F, [-0.16666667, -0.16666667, -0.16666667, -0.16666667, -0.16666667, -0.16666667, 0.16666667, 0.16666667, 0.16666667])
if GUI:
from Muscat.Actions.OpenInParaView import OpenInParaView
F.shape = (3, 3)
F = F.T
mesh.nodeFields["normalflux"] = F
OpenInParaView(mesh)
return "ok"
[docs]
def CheckIntegrityKF(edim, sdim, testCase):
from Muscat.MeshTools.MeshCreationTools import CreateMeshOfTriangles, CreateMeshOf
from Muscat.FE.MaterialHelp import HookeIso
if edim == 1:
nu = 0.25
# https://www.colorado.edu/engineering/CAS/courses.d/IFEM.d/IFEM.Ch20.d/IFEM.Ch20.pdf
if sdim == 1:
E = 1
A = 1
points = np.array([[0, ], [2, ], ])
L = np.sqrt(np.sum((points[1, :] - points[0, :])**2))
K = np.array([[E], ])
KValidation = (E*A/L)*np.array([[1, -1], [-1, 1]])
permut = None
elif sdim == 2:
E = 1000.
A = 5.
points = np.array([[0, 0], [30, 40], ])
L = np.sqrt(np.sum((points[1, :] - points[0, :])**2))
K = A*E*np.array([[1, 1, 0], [1, 1, 0], [0, 0, 0]])
KValidation = np.array([[36, 48, -36, -48],
[48, 64, -48, -64],
[-36, -48, 36, 48],
[-48, -64, 48, 64]])
permut = [0, 2, 1, 3]
KValidation = KValidation[permut, :][:, permut]
elif sdim == 3:
points = np.array([[0, 0, 0], [2., 3., 6.], ])
L = np.sqrt(np.sum((points[1, :] - points[0, :])**2))
E = 343.
A = 10.
K = A*E*np.array([[1, 1, 1, 0, 0, 0],
[1, 1, 1, 0, 0, 0],
[1, 1, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0], ])
permut = [0, 3, 1, 4, 2, 5]
KValidation = np.array([[40, 60, 120, -40, -60, -120],
[60, 90, 180, -60, -90, -180],
[120, 180, 360, -120, -180, -360],
[-40, -60, -120, 40, 60, 120],
[-60, -90, -180, 60, 90, 180],
[-120, -180, -360, 120, 180, 360]])
KValidation = KValidation[permut, :][:, permut]
mesh = CreateMeshOf(points, [[0, 1], ], ED.Bar_2)
elif edim == 2:
if sdim == 2:
if testCase[0] == "A":
# http://www.unm.edu/~bgreen/ME360/2D%20Triangular%20Elements.pdf
points = [[3, 0], [3, 2], [0, 0], [0, 2]]
# ,[3,2,1][0,1,2]
mesh = CreateMeshOfTriangles(points, [[0, 1, 2], [3, 2, 1]])
E = 3.0/2
nu = 0.25
K = HookeIso(E, nu, dim=2)
permut = [0, 2, 6, 4, 1, 3, 7, 5]
KValidation = np.array([[0.9833333333333333333, -0.5, -.45, .2, 0, 0, -0.5333333333333333333, 0.3],
[-0.5, 1.4, .3, -1.2, 0, 0, .2, -.2],
[-0.45, .3, 0.9833333333333333333, 0, -0.5333333333333333333, 0.2, 0, -0.5],
[.2, -1.2, 0, 1.4, .3, -.2, -0.5, 0],
[0, 0, -0.5333333333333333333, .3, 0.983333333333333333, -0.5, -0.45, .2],
[0, 0, 0.2, -0.2, -0.5, 1.4, .3, -1.2],
[-0.5333333333333333333, 0.2, 0.0, -0.5, -0.45, 0.3, 0.9833333333333333333, 0],
[0.3, -0.2, -0.5, 0, 0.2, -1.2, 0, 1.4]])
# ^^^
# attention il y a un erreur dans le pdf (ce termet est possitive dans le pdf)
KValidation = KValidation[permut, :][:, permut]
elif testCase[0] == "B":
# https://www.colorado.edu/engineering/CAS/courses.d/IFEM.d/IFEM.Ch15.d/IFEM.Ch15.pdf
# pages 15-11
points = [[0, 0], [3, 1], [2, 2]]
mesh = CreateMeshOfTriangles(points, [[0, 1, 2], ])
#mesh.GetElementsOfType(ED.Triangle_3).tags.CreateTag("2D").SetIds(np.arange(mesh.GetElementsOfType(ED.Triangle_3).GetNumberOfElements() ) )
K = np.array([[8, 2, 0], [2, 8, 0], [0, 0, 3]])*8
permut = [0, 2, 4, 1, 3, 5]
KValidation = np.array([[11, 5, -10, -2, -1, -3],
[5, 11, 2, 10, -7, -21],
[-10, 2, 44, -20, -34, 18],
[-2, 10, -20, 44, 22, -54],
[-1, -7, -34, 22, 35, -15],
[-3, -21, 18, -54, -15, 75]])
KValidation = KValidation[permut, :][:, permut]
else:
raise RuntimeError
elif sdim == 3:
if testCase[0] == "A":
points = [[0, 0, 0], [6, -8, 5], [6, 2, 3], [0, 5, 2]]
mesh = CreateMeshOf(points, [[0, 1, 2, 3]], ED.Quadrangle_4)
E = 3.0/2
nu = 0.25
K = HookeIso(E, nu, dim=sdim)
mesh.GetElementsOfType(ED.Quadrangle_4).tags.CreateTag("2D").SetIds(np.arange(mesh.GetElementsOfType(ED.Quadrangle_4).GetNumberOfElements()))
permut = [0, 3, 6, 1, 4, 7, 2, 5, 8]
else:
raise RuntimeError
elif edim == 3:
# https://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch09.d/AFEM.Ch09.pdf
# page 9-17
points = [[2., 3., 4], [6, 3, 2], [2, 5, 1], [4, 3, 6]]
mesh = CreateMeshOf(points, [[0, 1, 2, 3]], ED.Tetrahedron_4)
E = 480
nu = 1./3.
K = HookeIso(E, nu, dim=sdim)
#mesh.GetElementsOfType(ED.Tetrahedron_4).tags.CreateTag("3D").SetIds(np.arange(mesh.GetElementsOfType(ED.Tetrahedron_4).GetNumberOfElements() ) )
permut = [0, 3, 6, 9, 1, 4, 7, 10, 2, 5, 8, 11]
KValidation = np.array([[745, 540, 120, -5, 30, 60, -270, -240, 0, -470, -330, -180],
[540, 1720, 270, -120, 520, 210, -120, -1080, -60, -300, -1160, -420],
[120, 270, 565, 0, 150, 175, 0, -120, -270, -120, -300, -470],
[-5, -120, 0, 145, -90, -60, -90, 120, 0, -50, 90, 60],
[30, 520, 150, -90, 220, 90, 60, -360, -60, 0, -380, -180],
[60, 210, 175, -60, 90, 145, 0, -120, -90, 0, -180, -230],
[-270, -120, 0, -90, 60, 0, 180, 0, 0, 180, 60, 0],
[-240, -1080, -120, 120, -360, -120, 0, 720, 0, 120, 720, 240],
[0, -60, -270, 0, -60, -90, 0, 0, 180, 0, 120, 180],
[-470, -300, -120, -50, 0, 0, 180, 120, 0, 340, 180, 120],
[-330, -1160, -300, 90, -380, -180, 60, 720, 120, 180, 820, 360],
[-180, -420, -470, 60, -180, -230, 0, 240, 180, 120, 360, 520]])
KValidation = KValidation[permut, :][:, permut]
else:
raise RuntimeError
# CompureVolume(mesh)
mesh.ConvertDataForNativeTreatment()
space = LagrangeSpaceGeo
numbering = ComputeDofNumbering(mesh, space)
if sdim == 1:
dofs = ["u"]
else:
dofs = ["u_" + str(i) for i in range(sdim)]
spaces = [space]*sdim
numberings = [numbering]*sdim
from Muscat.FE.SymWeakForm import Strain
from Muscat.FE.SymWeakForm import ToVoigtEpsilon
from Muscat.FE.WeakForms.NumericalWeakForm import SymWeakToNumWeak
u = GetField("u", sdim, sdim=sdim)
ut = GetTestField("u", sdim, sdim=sdim)
weak = ToVoigtEpsilon(Strain(u, sdim)).T@K@ToVoigtEpsilon(Strain(ut, sdim))
constants = {}
fields = []
startt = time.time()
elemFilt = ElementFilter(eTag=(str(edim)+"D"))
unknownFields = [FEField(dofs[n], mesh=mesh, space=spaces[n], numbering=numberings[n]) for n in range(len(dofs))]
kInteg, _rhs = IntegrateGeneral(mesh=mesh, wform=weak, constants=constants, fields=fields,
unknownFields=unknownFields,
elementFilter=elemFilt)
stopt = time.time() - startt
print(stopt)
kDense = kInteg.todense()
permut = []
offset = 0
for f in unknownFields:
for p in range(mesh.GetNumberOfNodes()):
permut.append(f.numbering.GetDofOfPoint(p)+offset)
offset += f.numbering.size
if permut is not None:
kDense = kDense[:, permut][permut, :]
error = np.sum(abs(kDense-KValidation))/(np.sum(abs(KValidation)))
if error > 1e-14 or error is np.nan:
print((edim, sdim, testCase))
print("K Validation")
print(KValidation)
print("K Calcul")
print(kDense)
print("ERROR : ")
print(error)
print(kDense-KValidation)
print(numbering)
return "KO"
return "ok"
[docs]
def CheckIntegrityIntegrationWithIntegrationPointField(GUI: bool = False):
from Muscat.FE.IntegrationRules import LagrangeP1Quadrature
from Muscat.MeshTools.MeshCreationTools import CreateCube
mesh = CreateCube([2., 3., 4.], [-1.0, -1.0, -1.0], [2./10, 2./10, 2./10])
mesh.ConvertDataForNativeTreatment()
space, numberings, _offset, _NGauss = PrepareFEComputation(mesh, numberOfComponents=1)
rhoField = IPField("rho", mesh=mesh, rule=LagrangeP1Quadrature)
factor: MuscatFloat = .1
rhoField.Allocate(factor)
tField = GetField("T", 1,sdim=3)
rho = GetField("rho", 1,sdim=3)
tFieldTest = GetTestField("T", 1,sdim=3)
weakForm = tField.T*tFieldTest + rho.T*tFieldTest
constants = {}
fields = [rhoField]
unknownFields = [FEField("T", mesh=mesh, space=space, numbering=numberings[0])]
startt = time.time()
elemFilt = ElementFilter(eTag="3D")
K, F = IntegrateGeneral(mesh=mesh,
wform=weakForm,
constants=constants,
fields=fields,
unknownFields=unknownFields,
integrationRuleName="LagrangeP1Quadrature",
elementFilter=elemFilt)
intClass = GetIntegrator()
intClass.SetMesh(mesh)
intClass.SetIntegrationRule(integrationRuleName="LagrangeP1Quadrature")
intClass.SetOnlyEvaluation(False)
intClass.SetElementFilter(elemFilt)
intClass.SetConstants(constants)
intClass.SetUnknownFields(unknownFields)
intClass.SetTestFields(None)
intClass.SetExtraFields(fields)
intClass.SetWeakForm(weakForm)
intClass.PreStartCheck()
intClass.Allocate()
intClass.Compute()
matrixKII = intClass.GetK()
vectorFII = intClass.GetRhs()
_stopt = time.time() - startt
volk = np.sum(K)
print("volume (k): " + str(volk))
print("volume class (k): " + str(np.sum(matrixKII)))
volf = np.sum(F)
print("volume (f): " + str(volf))
print("volume class (f): " + str(np.sum(vectorFII)))
if abs(volk*factor - volf) >= volf/100000000:
return "KO"
volkII = np.sum(intClass.GetK())
volfII = np.sum(intClass.GetRhs())
if abs(volkII*factor - volfII) >= volfII/100000000:
return "KO"
return "ok"
[docs]
def CheckIntegrity(GUI: bool = False):
import Muscat.FE.Integration as BTFEI
saveOldeStateOfUseCpp = BTFEI.UseCpp
saveOldeStateOfUseKokkos = BTFEI.UseKokkos
toCheck = [CheckIntegrityNormalFlux,
CheckIntegrityIntegrationWithIntegrationPointField,
CheckIntegrityTangentAndRHS,
CheckIntegrityConstantRectilinearIntegrationII,
CheckIntegrityConstantRectilinearIntegration]
from Muscat.Helpers.CheckTools import RunListOfCheckIntegrities
#(False,False),
for k,c in [(False,True),(True,True)]:
BTFEI.UseKokkos = k
BTFEI.UseCpp = c
print(f" -------- UseCpp={BTFEI.UseCpp} -------- UseKokkos ={BTFEI.UseKokkos} -------- ")
res = RunListOfCheckIntegrities(toCheck)
if res.lower() != "ok":
BTFEI.UseCpp = saveOldeStateOfUseCpp
BTFEI.UseKokkos = saveOldeStateOfUseKokkos
return res
BTFEI.UseCpp = saveOldeStateOfUseCpp
BTFEI.UseKokkos = saveOldeStateOfUseKokkos
res = RunListOfCheckIntegrities([CheckIntegritySpeedUp])
if res.lower() != "ok":
return res
return "ok"
from Muscat.Helpers.Cache import CachedResultDecorator
[docs]
def GetCheckIntegrityMesh(dimensions, spacing, origin):
from Muscat.MeshTools.MeshCreationTools import CreateConstantRectilinearMesh
mesh = CreateConstantRectilinearMesh(dimensions=dimensions, spacing=spacing, origin=origin)
#force the generation of the connectivity
mesh.GetElementsOfType(ED.Hexahedron_8).connectivity
return mesh
[docs]
def CheckIntegritySpeedUpCore(s=20, toTests = [True, True, True]):
if toTests[0] == False and toTests[1] == False and toTests[2] == False:
return [None, None, None]
from Muscat.MeshTools.MeshCreationTools import CreateConstantRectilinearMesh
nNodesX = s
nNodesY = s
nNodesZ = s
mesh = GetCheckIntegrityMesh(dimensions=[nNodesX, nNodesY, nNodesZ], spacing=[1./(nNodesX-1), 1./(nNodesY-1), 1./(nNodesZ-1)], origin=[0, 0, 0])
#force the generation of the connectivity
mesh.GetElementsOfType(ED.Hexahedron_8).connectivity
from Muscat.FE.SymWeakForm import testcharacter
tFieldTest = GetTestField("T", 1, sdim=3)
Test = GetField("T", 1, sdim=3)
weakForm = 1*Test*tFieldTest
constants = {}
fields = []
uFields= [FEField("T", mesh=mesh)]
testFields = [FEField("T"+testcharacter, mesh=mesh)]
elemFilter = ElementFilter(dimensionality=3)
import Muscat.FE.Integration as BTFEI
saveOldeStateOfUseCpp = BTFEI.UseCpp
saveOldeStateOfUseKokkos = BTFEI.UseKokkos
dtime = [ ]
for testid, (k,c) in enumerate([(False,False),(False,True),(True,True)]):
if toTests[testid] == False:
dtime.append(None)
continue
BTFEI.UseKokkos = k
BTFEI.UseCpp = c
startTime = time.time()
f = IntegrateGeneral(mesh=mesh,
wform=weakForm,
constants=constants,
fields=fields,
unknownFields=uFields,
testFields=testFields,
integrationRuleName="LagrangeP1Quadrature",
elementFilter=elemFilter)[1]
volume = np.sum(f)
dtime.append(time.time() - startTime)
BTFEI.UseCpp = saveOldeStateOfUseCpp
BTFEI.UseKokkos = saveOldeStateOfUseKokkos
return dtime
[docs]
def CheckIntegritySpeedUp(GUI: bool = False):
dtime = CheckIntegritySpeedUpCore(30)
from Muscat.Helpers.TextFormatHelper import TFormat
def PrintRow(datarow):
res = "|"
for d in datarow:
if type(d) is str:
res += TFormat.Center(str(d), fill=" ", width=20) + "|"
else:
res += TFormat.Center("%.4f" % d, fill=" ", width=20) + "|"
print(res)
PrintRow(["Python [s}", "CPP [s]", "Speedup CPP", "Kokkos [s]", "speedup Kokkos"])
if dtime[0] != None and dtime[1] != None:
cppSpeedUp = dtime[0]/dtime[1]
else:
cppSpeedUp = "NA"
if dtime[0] != None and dtime[2] != None:
cppKokkosSpeedUp = dtime[0]/dtime[2]
else:
cppKokkosSpeedUp = "NA"
PrintRow([dtime[0], dtime[1], cppSpeedUp, dtime[2], cppKokkosSpeedUp ])
return "ok"
[docs]
def CheckIntegrityTangentAndRHS(GUI: bool = False):
problems = [(1, 1, "A bars 1D"),
(1, 2, "A bars 2D"),
(1, 3, "A bars 3D"),
(2, 2, "A Triangles in 2D"),
(2, 2, "B Triangle in 2D"),
#(2,3,"A Triangle in 3D"),
(3, 3, "A tet in 3D"),
]
startt = time.time()
for edim, sdim, testCase in problems:
print("*-"*40)
print((edim, sdim, testCase))
res = CheckIntegrityKF(edim=edim, sdim=sdim, testCase=testCase)
if res.lower() != "ok":
return res
else:
print(" --- CheckIntegrityKF integration -- OK ")
stopt = time.time() - startt
print("Total time : ")
print(stopt)
print("ALL ok")
return "ok"
[docs]
def CheckIntegrityConstantRectilinearIntegrationII(GUI: bool = False):
print("Test integrator on Constant rectilinear mesh")
from Muscat.MeshTools.MeshCreationTools import CreateConstantRectilinearMesh
nNodesX = 3
nNodesY = 4
nNodesZ = 5
mesh = CreateConstantRectilinearMesh(dimensions=[nNodesX, nNodesY, nNodesZ], spacing=[1./(nNodesX-1), 1./(nNodesY-1), 1./(nNodesZ-1)], origin=[0, 0, 0])
from Muscat.FE.Spaces.FESpaces import ConstantSpaceGlobal
testNumbering = ComputeDofNumbering(mesh, ConstantSpaceGlobal)
from Muscat.FE.SymWeakForm import testcharacter
tFieldTest = GetTestField("T", 1, sdim=3)
weakForm = 1*tFieldTest
constants = {}
fields = []
testFields = [FEField("T"+testcharacter, mesh=mesh, space=ConstantSpaceGlobal, numbering=testNumbering)]
startt = time.time()
elemFilter = ElementFilter()
print(startt)
K, F = IntegrateGeneral(mesh=mesh,
wform=weakForm,
constants=constants,
fields=fields,
testFields=testFields,
integrationRuleName="LagrangeP1Quadrature",
elementFilter=elemFilter)
print("K ", K)
print("F ", F)
return "ok"
[docs]
def CheckIntegrityConstantRectilinearIntegration(GUI: bool = False):
print("Test integrator on Constant rectilinear mesh")
from Muscat.MeshTools.MeshCreationTools import CreateConstantRectilinearMesh
nNodesY = 3
nNodesY = 4
nNodesZ = 5
myCRMesh = CreateConstantRectilinearMesh(dimensions=[nNodesY, nNodesY, nNodesZ], spacing=[1./(nNodesY-1), 1./(nNodesY-1), 1./(nNodesZ-1)], origin=[0, 0, 0])
space, numberings, _offset, _NGauss = PrepareFEComputation(myCRMesh, numberOfComponents=1)
tField = GetField("T", 1, sdim=3)
tFieldTest = GetTestField("T", 1, sdim=3)
weakForm = tField.T*tFieldTest + 1*tFieldTest
constants = {}
fields = []
unknownFields = [FEField("T", mesh=myCRMesh, space=space, numbering=numberings[0])]
startt = time.time()
elemFilter = ElementFilter()
print(startt)
K, F = IntegrateGeneral(mesh=myCRMesh,
wform=weakForm,
constants=constants,
fields=fields,
unknownFields=unknownFields,
integrationRuleName="LagrangeP1Quadrature",
elementFilter=elemFilter)
stopt = time.time() - startt
volk = np.sum(K)
print("volume (k): " + str(volk))
volf = np.sum(F)
print("volume (f): " + str(volf))
print(stopt)
if volk*1 - volf < volf/100000000:
return "ok"
print("volk: "+str(volk))
print("volf: "+str(volf))
return "Not Ok"
[docs]
def RunHeavySpeedUpCheck():
#from Muscat.Helpers.Kokkos.KokkosHelper import PrintConfigKokkos
#PrintConfigKokkos()
CheckIntegritySpeedUpCore(2)
datas = []
data = [0,0,0]
import os
OMP_NUM_THREADS = os.environ.get("OMP_NUM_THREADS",None)
from Muscat.Helpers.CPU import GetNumberOfAvailableCores
cores = GetNumberOfAvailableCores()
#print(f"{cores=}")
import matplotlib.pyplot as plt
import numpy as np
for t in np.logspace(0.4,3,num=25):
t = int(t)
toTest = [(d!=None and d<6 ) for d in data]
if np.all(np.array(toTest)==False):
print(", ".join(map(str,[cores, OMP_NUM_THREADS,t,(t-1)**3, None, None, None])))
break
data = [1e10,1e10,1e10]
for _try in range(2):
ldata = CheckIntegritySpeedUpCore(t, toTest)
ldata = [(d if d != None else 1e10 ) for d in ldata ]
data = np.minimum(data,ldata)
data = [ (d if d < 1e10 else None ) for d in data]
datas.append([cores, OMP_NUM_THREADS,t,(t-1)**3]+list(data))
print(", ".join(map(str,datas[-1])))
nbelements = [ d[ 3] for d in datas]
python = [ d[ 4] for d in datas]
cpp = [ d[ 5] for d in datas]
kokkos = [ d[ 6] for d in datas]
if __name__ == '__main__':
print(CheckIntegrity(False)) # pragma: no cover
#print(CheckIntegritySpeedUpCore(2, [False,False, True]))