Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
for ebNE in range(self.mesh.nExteriorElementBoundaries_global):
ebN = self.mesh.exteriorElementBoundariesArray[ebNE]
eN_global = self.mesh.elementBoundaryElementsArray[ebN, 0]
ebN_element = self.mesh.elementBoundaryLocalElementBoundariesArray[ebN, 0]
for i in range(self.mesh.nNodes_element):
if i != ebN_element:
I = self.mesh.elementNodesArray[eN_global, i]
self.internalNodes -= set([I])
self.nNodes_internal = len(self.internalNodes)
self.internalNodesArray = np.zeros((self.nNodes_internal,), 'i')
for nI, n in enumerate(self.internalNodes):
self.internalNodesArray[nI] = n
#
del self.internalNodes
self.internalNodes = None
logEvent("Updating local to global mappings", 2)
self.updateLocal2Global()
logEvent("Building time integration object", 2)
logEvent(memory("inflowBC, internalNodes,updateLocal2Global", "OneLevelTransport"), level=4)
# mwf for interpolating subgrid error for gradients etc
if self.stabilization and self.stabilization.usesGradientStabilization:
self.timeIntegration = TimeIntegrationClass(self, integrateInterpolationPoints=True)
else:
self.timeIntegration = TimeIntegrationClass(self)
if options is not None:
self.timeIntegration.setFromOptions(options)
logEvent(memory("TimeIntegration", "OneLevelTransport"), level=4)
logEvent("Calculating numerical quadrature formulas", 2)
self.calculateQuadrature()
self.setupFieldStrides()
argsDict["isDiffusiveFluxBoundary_u"] = self.ebqe[('diffusiveFlux_bc_flag', 0, 0)]
argsDict["ebqe_bc_diffusiveFlux_u_ext"] = self.ebqe[('diffusiveFlux_bc', 0, 0)]
argsDict["ebqe_phi"] = self.coefficients.ebqe_phi
argsDict["epsFact"] = self.coefficients.epsFact
argsDict["ebqe_dissipation"] = self.coefficients.ebqe_dissipation
argsDict["ebqe_porosity"] = self.coefficients.ebqe_porosity
argsDict["ebqe_u"] = self.ebqe[('u', 0)]
argsDict["ebqe_flux"] = self.ebqe[('advectiveFlux', 0)]
self.kappa.calculateResidual(argsDict)
if self.forceStrongConditions:
for dofN, g in list(self.dirichletConditionsForceDOF.DOFBoundaryConditionsDict.items()):
r[dofN] = 0
if self.stabilization:
self.stabilization.accumulateSubgridMassHistory(self.q)
prof.logEvent("Global residual", level=9, data=r)
# mwf decide if this is reasonable for keeping solver statistics
self.nonlinear_function_evaluations += 1
if self.globalResidualDummy is None:
self.globalResidualDummy = np.zeros(r.shape, 'd')
def outputRow(self, time):
""" Outputs a single row of currently calculated gauge data to self.file"""
assert self.isGaugeOwner
if self.isPointGauge or self.isLineGauge:
self.localQuantitiesBuf = np.concatenate([gaugesVec.getArray() for gaugesVec in
self.pointGaugeVecs]).astype(np.double)
logEvent("Sending local array of type %s and shape %s to root on comm %s" % (
str(self.localQuantitiesBuf.dtype), str(self.localQuantitiesBuf.shape), str(self.gaugeComm)), 9)
if self.gaugeComm.rank == 0:
logEvent("Receiving global array of type %s and shape %s on comm %s" % (
str(self.localQuantitiesBuf.dtype), str(self.globalQuantitiesBuf.shape), str(self.gaugeComm)), 9)
self.gaugeComm.Gatherv(sendbuf=[self.localQuantitiesBuf, MPI.DOUBLE],
recvbuf=[self.globalQuantitiesBuf, (self.globalQuantitiesCounts, None),
MPI.DOUBLE], root=0)
self.gaugeComm.Barrier()
if self.isLineIntegralGauge:
lineIntegralGaugeBuf = self.lineIntegralGaugesVec.getArray()
globalLineIntegralGaugeBuf = lineIntegralGaugeBuf.copy()
self.gaugeComm.Reduce(lineIntegralGaugeBuf, globalLineIntegralGaugeBuf, op=MPI.SUM)
else:
globalLineIntegralGaugeBuf = []
if self.gaugeComm.rank == 0:
self.packFraction=closure.packFraction
self.packMargin=closure.packMargin
self.maxFraction=closure.maxFraction
self.frFraction=closure.frFraction
self.sigmaC=closure.sigmaC
self.C3e=closure.C3e
self.C4e=closure.C4e
self.eR=closure.eR
self.fContact=closure.fContact
self.mContact=closure.mContact
self.nContact=closure.nContact
self.angFriction=closure.angFriction
self.vos_limiter = closure.vos_limiter
self.mu_fr_limiter = closure.mu_fr_limiter
self.sedFlag = 1
prof.logEvent("INFO: Loading parameters for sediment closure",2)
except:
self.aDarcy=-1.
self.betaForch=-1.
self.grain=-1.
self.packFraction=-1.
self.packMargin=-1.
self.maxFraction=-1.
self.frFraction=-1.
self.sigmaC=-1.
self.C3e=-1.
self.C4e=-1.
self.eR=-1.
self.fContact=-1.
self.mContact=-1.
self.nContact=-1.
self.angFriction=-1.
def globalConstantSolve(self, u, r):
U = 0.0
R = 0.0
J = 0.0
(R, J) = self.globalConstantRJ(u, r, U)
its = 0
log(" Mass Conservation Residual 0 ", level=3, data=R)
RNORM_OLD = fabs(R)
while ((fabs(R) > self.atol and its < self.maxIts) or its < 1):
U -= old_div(R, (J + 1.0e-8))
(R, J) = self.globalConstantRJ(u, r, U)
lsits = 0
while(fabs(R) > 0.99 * RNORM_OLD and lsits < self.maxLSits):
lsits += 1
U += (0.5)**lsits * (old_div(R, (J + 1.0e-8)))
(R, J) = self.globalConstantRJ(u, r, U)
its += 1
log(" Mass Conservation Residual " + repr(its)+" ", level=3, data=R)
self.u[0].dof.flat[:] = U
sfConfig = p0.domain.PUMIMesh.size_field_config()
logEvent("h-adapt mesh by calling AdaptPUMIMesh")
if(sfConfig=="pseudo"):
logEvent("Testing solution transfer and restart feature of adaptation. No actual mesh adaptation!")
else:
p0.domain.PUMIMesh.adaptPUMIMesh()
#code to suggest adapting until error is reduced;
#not fully baked and can lead to infinite loops of adaptation
#if(sfConfig=="ERM"):
# p0.domain.PUMIMesh.get_local_error()
# while(p0.domain.PUMIMesh.willAdapt()):
# p0.domain.PUMIMesh.adaptPUMIMesh()
# p0.domain.PUMIMesh.get_local_error()
logEvent("Converting PUMI mesh to Proteus")
#ibaned: PUMI conversion #2
#TODO: this code is nearly identical to
#PUMI conversion #1, they should be merged
#into a function
if p0.domain.nd == 3:
mesh = MeshTools.TetrahedralMesh()
else:
mesh = MeshTools.TriangularMesh()
mesh.convertFromPUMI(p0.domain.PUMIMesh,
p0.domain.faceList,
p0.domain.regList,
parallel = self.comm.size() > 1,
dim = p0.domain.nd)
self.PUMI2Proteus(mesh)
logEvent("Updating local to global mappings", 2)
self.updateLocal2Global()
logEvent("Building time integration object", 2)
logEvent(memory("inflowBC, internalNodes,updateLocal2Global",
"OneLevelTransport"), level=4)
# mwf for interpolating subgrid error for gradients etc
if self.stabilization and self.stabilization.usesGradientStabilization:
self.timeIntegration = TimeIntegrationClass(
self, integrateInterpolationPoints=True)
else:
self.timeIntegration = TimeIntegrationClass(self)
if options is not None:
self.timeIntegration.setFromOptions(options)
logEvent(memory("TimeIntegration", "OneLevelTransport"), level=4)
logEvent("Calculating numerical quadrature formulas", 2)
self.calculateQuadrature()
self.setupFieldStrides()
comm = Comm.get()
self.comm = comm
if comm.size() > 1:
assert numericalFluxType is not None and numericalFluxType.useWeakDirichletConditions, "You must use a numerical flux to apply weak boundary conditions for parallel runs"
logEvent(memory("stride+offset", "OneLevelTransport"), level=4)
if numericalFluxType is not None:
if options is None or options.periodicDirichletConditions is None:
self.numericalFlux = numericalFluxType(
self,
dofBoundaryConditionsSetterDict,
advectiveFluxBoundaryConditionsSetterDict,
diffusiveFluxBoundaryConditionsSetterDictDict)
self.elementBoundaryQuadraturePoints, self.ebqe['x'])
self.u[0].femSpace.elementMaps.getJacobianValuesGlobalExteriorTrace(
self.elementBoundaryQuadraturePoints,
self.ebqe['inverse(J)'],
self.ebqe['g'],
self.ebqe['sqrt(det(g))'],
self.ebqe['n'])
cfemIntegrals.calculateIntegrationWeights(
self.ebqe['sqrt(det(g))'], self.elementBoundaryQuadratureWeights[
('u', 0)], self.ebqe[
('dS_u', 0)])
#
# get physical locations of element boundary quadrature points
#
# assume all components live on the same mesh
log("initalizing basis info")
self.u[0].femSpace.elementMaps.getBasisValuesTraceRef(
self.elementBoundaryQuadraturePoints)
self.u[0].femSpace.elementMaps.getBasisGradientValuesTraceRef(
self.elementBoundaryQuadraturePoints)
self.u[0].femSpace.getBasisValuesTraceRef(
self.elementBoundaryQuadraturePoints)
self.u[0].femSpace.getBasisGradientValuesTraceRef(
self.elementBoundaryQuadraturePoints)
self.u[1].femSpace.getBasisValuesTraceRef(
self.elementBoundaryQuadraturePoints)
self.u[1].femSpace.getBasisGradientValuesTraceRef(
self.elementBoundaryQuadraturePoints)
self.u[0].femSpace.elementMaps.getValuesGlobalExteriorTrace(
self.elementBoundaryQuadraturePoints, self.ebqe['x'])
log("setting flux boundary conditions")
if not domainMoved:
ksp.buildResidual(r_work)
truenorm = r_work.norm()
if its == 0:
self.rnorm0 = truenorm
logEvent("NumericalAnalytics KSPSchurResidual: %12.5e" %(truenorm) )
logEvent("NumericalAnalytics KSPSchurResidual(relative): %12.5e" %(truenorm / self.rnorm0) )
logEvent(" KSP it %i norm(r) = %e norm(r)/|b| = %e ; atol=%e rtol=%e " % (its,
truenorm,
(truenorm/ self.rnorm0),
ksp.atol,
ksp.rtol))
return False
else:
logEvent("NumericalAnalytics KSPSchurResidual: %12.5e" %(truenorm) )
logEvent("NumericalAnalytics KSPSchurResidual(relative): %12.5e" %(truenorm / self.rnorm0) )
logEvent(" KSP it %i norm(r) = %e norm(r)/|b| = %e ; atol=%e rtol=%e " % (its,
truenorm,
(truenorm/ self.rnorm0),
ksp.atol,
ksp.rtol))
if truenorm < self.rnorm0*ksp.rtol:
return p4pyPETSc.KSP.ConvergedReason.CONVERGED_RTOL
if truenorm < ksp.atol:
return p4pyPETSc.KSP.ConvergedReason.CONVERGED_ATOL
return False
else:
mlMesh.generatePartitionedMeshFromPUMI(
mesh,n0.nLevels,
nLayersOfOverlap=n0.nLayersOfOverlapForParallel)
self.mlMesh_nList=[]
for p in self.pList:
self.mlMesh_nList.append(mlMesh)
if (p0.domain.PUMIMesh.size_field_config() == "isotropicProteus"):
mlMesh.meshList[0].subdomainMesh.size_field = numpy.ones((mlMesh.meshList[0].subdomainMesh.nNodes_global,1),'d')*1.0e-1
if (p0.domain.PUMIMesh.size_field_config() == 'anisotropicProteus'):
mlMesh.meshList[0].subdomainMesh.size_scale = numpy.ones((mlMesh.meshList[0].subdomainMesh.nNodes_global,3),'d')
mlMesh.meshList[0].subdomainMesh.size_frame = numpy.ones((mlMesh.meshList[0].subdomainMesh.nNodes_global,9),'d')
#may want to trigger garbage collection here
modelListOld = self.modelList
logEvent("Allocating models on new mesh")
self.allocateModels()
logEvent("Attach auxiliary variables to new models")
#(cut and pasted from init, need to cleanup)
self.simOutputList = []
self.auxiliaryVariables = {}
self.newAuxiliaryVariables = {}
if self.simFlagsList is not None:
for p, n, simFlags, model, index in zip(
self.pList,
self.nList,
self.simFlagsList,
self.modelList,
list(range(len(self.pList)))):
self.simOutputList.append(
SimTools.SimulationProcessor(
flags=simFlags,