Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
argsDict["q_n"] = self.q[('grad(u)', 0)]
argsDict["ebqe_u"] = self.ebqe[('u', 0)]
argsDict["ebqe_n"] = self.ebqe[('grad(u)', 0)]
argsDict["q_r"] = self.q[('r', 0)]
argsDict["q_porosity"] = self.coefficients.q_porosity
argsDict["offset_u"] = self.offset[0]
argsDict["stride_u"] = self.stride[0]
argsDict["globalResidual"] = r
argsDict["nExteriorElementBoundaries_global"] = self.mesh.nExteriorElementBoundaries_global
argsDict["exteriorElementBoundariesArray"] = self.mesh.exteriorElementBoundariesArray
argsDict["elementBoundaryElementsArray"] = self.mesh.elementBoundaryElementsArray
argsDict["elementBoundaryLocalElementBoundariesArray"] = self.mesh.elementBoundaryLocalElementBoundariesArray
self.mcorr.calculateResidual(argsDict,
self.coefficients.useExact)
logEvent("Global residual", level=9, data=r)
self.coefficients.massConservationError = fabs(globalSum(r[:self.mesh.nNodes_owned].sum()))
logEvent(" Mass Conservation Error: ", level=3, data=self.coefficients.massConservationError)
self.nonlinear_function_evaluations += 1
if self.globalResidualDummy is None:
self.globalResidualDummy = np.zeros(r.shape, 'd')
self.transport.hv_dof_old[:] = self.u_dof_last[2]
self.transport.heta_dof_old[:] = self.u_dof_last[3]
self.transport.hw_dof_old[:] = self.u_dof_last[4]
elif self.timeOrder == 2:
if self.lstage == 1:
logEvent("First stage of SSP22 method finished", level=4)
for ci in range(self.nc):
self.u_dof_lstage[ci][:] = self.transport.u[ci].dof
# Update u_dof_old
self.transport.h_dof_old[:] = self.u_dof_lstage[0]
self.transport.hu_dof_old[:] = self.u_dof_lstage[1]
self.transport.hv_dof_old[:] = self.u_dof_lstage[2]
self.transport.heta_dof_old[:] = self.u_dof_lstage[3]
self.transport.hw_dof_old[:] = self.u_dof_lstage[4]
else:
logEvent("Second stage of SSP22 method finished", level=4)
for ci in range(self.nc):
self.u_dof_lstage[ci][:] = self.transport.u[ci].dof
self.u_dof_lstage[ci][:] *= old_div(1., 2.)
self.u_dof_lstage[ci][:] += 1. / 2. * self.u_dof_last[ci]
# update solution to u[0].dof
self.transport.u[ci].dof[:] = self.u_dof_lstage[ci]
# Update u_dof_old
self.transport.h_dof_old[:] = self.u_dof_last[0]
self.transport.hu_dof_old[:] = self.u_dof_last[1]
self.transport.hv_dof_old[:] = self.u_dof_last[2]
self.transport.heta_dof_old[:] = self.u_dof_last[3]
self.transport.hw_dof_old[:] = self.u_dof_last[4]
else:
assert self.timeOrder == 1
logEvent("FE method finished", level=4)
def updateStage(self):
"""
Need to switch to use coefficients
"""
self.lstage += 1
assert self.timeOrder in [1, 2, 3]
assert self.lstage > 0 and self.lstage <= self.timeOrder
# print "within update stage...: ", self.lstage
if self.timeOrder == 3:
if self.lstage == 1:
logEvent("First stage of SSP33 method finished", level=4)
for ci in range(self.nc):
self.u_dof_lstage[ci][:] = self.transport.u[ci].dof
# update u_dof_old
self.transport.h_dof_old[:] = self.u_dof_lstage[0]
self.transport.hu_dof_old[:] = self.u_dof_lstage[1]
self.transport.hv_dof_old[:] = self.u_dof_lstage[2]
elif self.lstage == 2:
logEvent("Second stage of SSP33 method finished", level=4)
for ci in range(self.nc):
self.u_dof_lstage[ci][:] = self.transport.u[ci].dof
self.u_dof_lstage[ci] *= old_div(1., 4.)
self.u_dof_lstage[ci] += 3. / 4. * self.u_dof_last[ci]
# update u_dof_old
self.transport.h_dof_old[:] = self.u_dof_lstage[0]
self.transport.hu_dof_old[:] = self.u_dof_lstage[1]
self.transport.hv_dof_old[:] = self.u_dof_lstage[2]
self.vectors_elementBoundaryQuadrature = set()
self.tensors_elementBoundaryQuadrature = set()
#
# show quadrature
#
logEvent("Dumping quadrature shapes for model %s" % self.name, level=9)
logEvent("Element quadrature array (q)", level=9)
for (k, v) in list(self.q.items()):
logEvent(str((k, v.shape)), level=9)
logEvent("Element boundary quadrature (ebq)", level=9)
for (k, v) in list(self.ebq.items()):
logEvent(str((k, v.shape)), level=9)
logEvent("Global element boundary quadrature (ebq_global)", level=9)
for (k, v) in list(self.ebq_global.items()):
logEvent(str((k, v.shape)), level=9)
logEvent("Exterior element boundary quadrature (ebqe)", level=9)
for (k, v) in list(self.ebqe.items()):
logEvent(str((k, v.shape)), level=9)
logEvent("Interpolation points for nonlinear diffusion potential (phi_ip)", level=9)
for (k, v) in list(self.phi_ip.items()):
logEvent(str((k, v.shape)), level=9)
#
# allocate residual and Jacobian storage
#
#
# allocate residual and Jacobian storage
#
self.elementResidual = [np.zeros(
(self.mesh.nElements_global,
self.nDOF_test_element[ci]),
'd')]
self.inflowBoundaryBC = {}
self.numericalFlux.isDOFBoundary[0],
self.numericalFlux.isDOFBoundary[1],
self.numericalFlux.isDOFBoundary[2],
self.ebqe[('stressFlux_bc_flag',0)],
self.ebqe[('stressFlux_bc_flag',1)],
self.ebqe[('stressFlux_bc_flag',2)],
self.csrColumnOffsets_eb[(0,0)],
self.csrColumnOffsets_eb[(0,1)],
self.csrColumnOffsets_eb[(0,2)],
self.csrColumnOffsets_eb[(1,0)],
self.csrColumnOffsets_eb[(1,1)],
self.csrColumnOffsets_eb[(1,2)],
self.csrColumnOffsets_eb[(2,0)],
self.csrColumnOffsets_eb[(2,1)],
self.csrColumnOffsets_eb[(2,2)])
logEvent("Jacobian ",level=10,data=jacobian)
if self.forceStrongConditions:
scaling = 1.0#probably want to add some scaling to match non-dirichlet diagonals in linear system
for cj in range(self.nc):
for dofN in list(self.dirichletConditionsForceDOF[cj].DOFBoundaryConditionsDict.keys()):
global_dofN = self.offset[cj]+self.stride[cj]*dofN
for i in range(self.rowptr[global_dofN],self.rowptr[global_dofN+1]):
if (self.colind[i] == global_dofN):
self.nzval[i] = scaling
else:
self.nzval[i] = 0.0
#mwf decide if this is reasonable for solver statistics
self.nonlinear_function_jacobian_evaluations += 1
#jacobian.fwrite("jacobian_p"+`self.nonlinear_function_jacobian_evaluations`)
return jacobian
def calculateElementQuadrature(self):
self.dt_model = m.timeIntegration.dt
if self.nSteps >= self.nStepsOsher: # start ramping up the time step
self.dt_model = self.dt_model * self.red_ratio
#logEvent("Osher-PsiTC dt %12.5e" %(self.dt_model),level=1)
self.set_dt_allLevels()
# physical time step
self.t_model = self.substeps[0]
self.substeps.append(self.substeps[0])
logEvent("Osher-PsiTC iteration %d dt = %12.5e |res| = %12.5e %g " % (self.nSteps, self.dt_model, res, (old_div(res, self.res0)) * 100.0), level=1)
elif self.nSteps >= self.nStepsMax:
logEvent("Osher-PsiTC DID NOT Converge |res| = %12.5e but quitting anyway" % (res,))
logEvent("Osher-PsiTC tolerance %12.5e " % (self.res0 * self.rtol + self.atol,))
self.nSteps = 0
else:
logEvent("Osher-PsiTC converged |res| = %12.5e %12.5e" % (res, ssError * 100.0))
logEvent("Osher-PsiTC tolerance %12.5e " % (self.res0 * self.rtol + self.atol,))
self.nSteps = 0
self.ebqe[('stressFlux_bc_flag',1)] = np.zeros((self.mesh.nExteriorElementBoundaries_global,self.nElementBoundaryQuadraturePoints_elementBoundary),'i')
self.ebqe[('stressFlux_bc_flag',2)] = np.zeros((self.mesh.nExteriorElementBoundaries_global,self.nElementBoundaryQuadraturePoints_elementBoundary),'i')
self.ebqe[('stressFlux_bc',0)] = np.zeros((self.mesh.nExteriorElementBoundaries_global,self.nElementBoundaryQuadraturePoints_elementBoundary),'d')
self.ebqe[('stressFlux_bc',1)] = np.zeros((self.mesh.nExteriorElementBoundaries_global,self.nElementBoundaryQuadraturePoints_elementBoundary),'d')
self.ebqe[('stressFlux_bc',2)] = np.zeros((self.mesh.nExteriorElementBoundaries_global,self.nElementBoundaryQuadraturePoints_elementBoundary),'d')
self.points_elementBoundaryQuadrature= set()
self.scalars_elementBoundaryQuadrature= set([('u',ci) for ci in range(self.nc)])
self.vectors_elementBoundaryQuadrature= set()
self.tensors_elementBoundaryQuadrature= set()
#
#show quadrature
#
logEvent("Dumping quadrature shapes for model %s" % self.name,level=9)
logEvent("Element quadrature array (q)", level=9)
for (k,v) in list(self.q.items()): logEvent(str((k,v.shape)),level=9)
logEvent("Element boundary quadrature (ebq)",level=9)
for (k,v) in list(self.ebq.items()): logEvent(str((k,v.shape)),level=9)
logEvent("Global element boundary quadrature (ebq_global)",level=9)
for (k,v) in list(self.ebq_global.items()): logEvent(str((k,v.shape)),level=9)
logEvent("Exterior element boundary quadrature (ebqe)",level=9)
for (k,v) in list(self.ebqe.items()): logEvent(str((k,v.shape)),level=9)
logEvent("Interpolation points for nonlinear diffusion potential (phi_ip)",level=9)
for (k,v) in list(self.phi_ip.items()): logEvent(str((k,v.shape)),level=9)
#
# allocate residual and Jacobian storage
#
self.elementResidual = [np.zeros(
(self.mesh.nElements_global,
self.nDOF_test_element[ci]),
'd') for ci in range(self.nc)]
self.elementSpatialResidual = [np.zeros(
(self.mesh.nElements_global,
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)
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()
# hEps: this is use to regularize the flux and re-define the dry states
self.hEps = None
self.hReg = None
self.ML = None # lumped mass matrix
self.MC_global = None # consistent mass matrix
# Global C Matrices (mql)
self.cterm_global = None
self.cterm_transpose_global = None
# For FCT
self.extendedSourceTerm_hu=None
self.extendedSourceTerm_hv=None
self.new_SourceTerm_hu = None
#
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 != None:
self.timeIntegration.setFromOptions(options)
logEvent(memory("TimeIntegration","OneLevelTransport"),level=4)
logEvent("Calculating numerical quadrature formulas",2)
self.calculateQuadrature()
#lay out components/equations contiguously for now
self.offset = [0]
for ci in range(1,self.nc):
self.offset += [self.offset[ci-1]+self.nFreeDOF_global[ci-1]]
self.stride = [1 for ci in range(self.nc)]
#use contiguous layout of components for parallel, requires weak DBC's
comm = Comm.get()
self.comm=comm
if comm.size() > 1:
assert numericalFluxType != None and numericalFluxType.useWeakDirichletConditions,"You must use a numerical flux to apply weak boundary conditions for parallel runs"
self.offset = [0]
for ci in range(1,self.nc):
self.offset += [ci]
self.stride = [self.nc for ci in range(self.nc)]