Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def getJacobian(self, jacobian):
cfemIntegrals.zeroJacobian_CSR(self.nNonzerosInJacobian, jacobian)
argsDict = cArgumentsDict.ArgumentsDict()
argsDict["mesh_trial_ref"] = self.u[0].femSpace.elementMaps.psi
argsDict["mesh_grad_trial_ref"] = self.u[0].femSpace.elementMaps.grad_psi
argsDict["mesh_dof"] = self.mesh.nodeArray
argsDict["mesh_l2g"] = self.mesh.elementNodesArray
argsDict["x_ref"] = self.elementQuadraturePoints
argsDict["dV_ref"] = self.elementQuadratureWeights[('u', 0)]
argsDict["u_trial_ref"] = self.u[0].femSpace.psi
argsDict["u_grad_trial_ref"] = self.u[0].femSpace.grad_psi
argsDict["u_test_ref"] = self.u[0].femSpace.psi
argsDict["u_grad_test_ref"] = self.u[0].femSpace.grad_psi
argsDict["mesh_trial_trace_ref"] = self.u[0].femSpace.elementMaps.psi_trace
argsDict["mesh_grad_trial_trace_ref"] = self.u[0].femSpace.elementMaps.grad_psi_trace
argsDict["dS_ref"] = self.elementBoundaryQuadratureWeights[('u', 0)]
argsDict["u_trial_trace_ref"] = self.u[0].femSpace.psi_trace
argsDict["u_grad_trial_trace_ref"] = self.u[0].femSpace.grad_psi_trace
argsDict["u_test_trace_ref"] = self.u[0].femSpace.psi_trace
def choose_dt(self):
maxCFL = 1.0e-6
# COMPUTE edge_based_cfl
rowptr_cMatrix, colind_cMatrix, Cx = self.transport.cterm_global[0].getCSRrepresentation()
rowptr_cMatrix, colind_cMatrix, Cy = self.transport.cterm_global[1].getCSRrepresentation()
rowptr_cMatrix, colind_cMatrix, CTx = self.transport.cterm_global_transpose[0].getCSRrepresentation()
rowptr_cMatrix, colind_cMatrix, CTy = self.transport.cterm_global_transpose[1].getCSRrepresentation()
numDOFsPerEqn = self.transport.u[0].dof.size
argsDict = cArgumentsDict.ArgumentsDict()
argsDict["g"] = self.transport.coefficients.g
argsDict["numDOFsPerEqn"] = numDOFsPerEqn
argsDict["lumped_mass_matrix"] = self.transport.ML
argsDict["h_dof_old"] = self.transport.u[0].dof
argsDict["hu_dof_old"] = self.transport.u[1].dof
argsDict["hv_dof_old"] = self.transport.u[2].dof
argsDict["heta_dof_old"] = self.transport.u[3].dof
argsDict["b_dof"] = self.transport.coefficients.b.dof
argsDict["csrRowIndeces_DofLoops"] = rowptr_cMatrix
argsDict["csrColumnOffsets_DofLoops"] = colind_cMatrix
argsDict["hEps"] = self.transport.hEps
argsDict["hReg"] = self.transport.hReg
argsDict["Cx"] = Cx
argsDict["Cy"] = Cy
argsDict["CTx"] = CTx
argsDict["CTy"] = CTy
def elementConstantSolve(self, u, r):
import pdb
import copy
"""
Calculate the element residuals and add in to the global residual
"""
r.fill(0.0)
# Load the unknowns into the finite element dof
self.setUnknowns(u)
# no flux boundary conditions
argsDict = cArgumentsDict.ArgumentsDict()
argsDict["mesh_trial_ref"] = self.u[0].femSpace.elementMaps.psi
argsDict["mesh_grad_trial_ref"] = self.u[0].femSpace.elementMaps.grad_psi
argsDict["mesh_dof"] = self.mesh.nodeArray
argsDict["mesh_l2g"] = self.mesh.elementNodesArray
argsDict["dV_ref"] = self.elementQuadratureWeights[('u', 0)]
argsDict["u_trial_ref"] = self.u[0].femSpace.psi
argsDict["u_grad_trial_ref"] = self.u[0].femSpace.grad_psi
argsDict["u_test_ref"] = self.u[0].femSpace.psi
argsDict["u_grad_test_ref"] = self.u[0].femSpace.grad_psi
argsDict["mesh_trial_trace_ref"] = self.u[0].femSpace.elementMaps.psi_trace
argsDict["mesh_grad_trial_trace_ref"] = self.u[0].femSpace.elementMaps.grad_psi_trace
argsDict["dS_ref"] = self.elementBoundaryQuadratureWeights[('u', 0)]
argsDict["u_trial_trace_ref"] = self.u[0].femSpace.psi_trace
argsDict["u_grad_trial_trace_ref"] = self.u[0].femSpace.grad_psi_trace
argsDict["u_test_trace_ref"] = self.u[0].femSpace.psi_trace
argsDict["u_grad_test_trace_ref"] = self.u[0].femSpace.grad_psi_trace
def FCTStep(self):
rowptr, colind, MassMatrix = self.MC_global.getCSRrepresentation()
limited_solution = np.zeros((len(rowptr) - 1),'d')
self.u_free_dof_stage_0_l = np.zeros((len(rowptr) - 1),'d')
self.timeIntegration.u_dof_stage[0][self.timeIntegration.lstage], # soln
fromFreeToGlobal=0
cfemIntegrals.copyBetweenFreeUnknownsAndGlobalUnknowns(fromFreeToGlobal,
self.offset[0],
self.stride[0],
self.dirichletConditions[0].global2freeGlobal_global_dofs,
self.dirichletConditions[0].global2freeGlobal_free_dofs,
self.u_free_dof_stage_0_l,
self.timeIntegration.u_dof_stage[0][self.timeIntegration.lstage])
argsDict = cArgumentsDict.ArgumentsDict()
argsDict["NNZ"] = self.nnz
argsDict["numDOFs"] = len(rowptr) - 1
argsDict["lumped_mass_matrix"] = self.ML
argsDict["soln"] = self.u_free_dof_stage_0_l
argsDict["solH"] = self.timeIntegration.u
argsDict["uLow"] = self.uLow
argsDict["limited_solution"] = limited_solution
argsDict["csrRowIndeces_DofLoops"] = rowptr
argsDict["csrColumnOffsets_DofLoops"] = colind
argsDict["MassMatrix"] = MassMatrix
argsDict["dt_times_dH_minus_dL"] = self.dt_times_dC_minus_dL
argsDict["min_u_bc"] = self.min_u_bc
argsDict["max_u_bc"] = self.max_u_bc
argsDict["LUMPED_MASS_MATRIX"] = self.coefficients.LUMPED_MASS_MATRIX
self.vos.FCTStep(argsDict)
self.timeIntegration.u[:] = limited_solution
def FCTStep(self):
rowptr, colind, MassMatrix = self.MassMatrix.getCSRrepresentation()
if (self.limited_L2p_vof_mass_correction is None):
self.limited_L2p_vof_mass_correction = np.zeros(self.LumpedMassMatrix.size, 'd')
argsDict = cArgumentsDict.ArgumentsDict()
argsDict["NNZ"] = self.nnz
argsDict["numDOFs"] = len(rowptr) - 1
argsDict["lumped_mass_matrix"] = self.LumpedMassMatrix
argsDict["solH"] = self.L2p_vof_mass_correction
argsDict["solL"] = self.lumped_L2p_vof_mass_correction
argsDict["limited_solution"] = self.limited_L2p_vof_mass_correction
argsDict["csrRowIndeces_DofLoops"] = rowptr
argsDict["csrColumnOffsets_DofLoops"] = colind
argsDict["matrix"] = MassMatrix
self.mcorr.FCTStep(argsDict)
self.ebqe[('advectiveFlux_bc', 0)][t[0], t[1]] = g(self.ebqe[('x')][t[0], t[1]], self.timeIntegration.t)
self.ebqe[('advectiveFlux_bc_flag', 0)][t[0], t[1]] = 1
for ck, diffusiveFluxBoundaryConditionsDict in list(self.fluxBoundaryConditionsObjectsDict[0].diffusiveFluxBoundaryConditionsDictDict.items()):
for t, g in list(diffusiveFluxBoundaryConditionsDict.items()):
self.ebqe[('diffusiveFlux_bc', ck, 0)][t[0], t[1]] = g(self.ebqe[('x')][t[0], t[1]], self.timeIntegration.t)
self.ebqe[('diffusiveFlux_bc_flag', ck, 0)][t[0], t[1]] = 1
# self.shockCapturing.lag=True
if self.forceStrongConditions:
for dofN, g in list(self.dirichletConditionsForceDOF.DOFBoundaryConditionsDict.items()):
self.u[0].dof[dofN] = g(self.dirichletConditionsForceDOF.DOFBoundaryPointDict[dofN], self.timeIntegration.t)
#
# mwf debug
#import pdb
# pdb.set_trace()
argsDict = cArgumentsDict.ArgumentsDict()
argsDict["mesh_trial_ref"] = self.u[0].femSpace.elementMaps.psi
argsDict["mesh_grad_trial_ref"] = self.u[0].femSpace.elementMaps.grad_psi
argsDict["mesh_dof"] = self.mesh.nodeArray
argsDict["mesh_velocity_dof"] = self.mesh.nodeVelocityArray
argsDict["MOVING_DOMAIN"] = self.MOVING_DOMAIN
argsDict["mesh_l2g"] = self.mesh.elementNodesArray
argsDict["dV_ref"] = self.elementQuadratureWeights[('u', 0)]
argsDict["u_trial_ref"] = self.u[0].femSpace.psi
argsDict["u_grad_trial_ref"] = self.u[0].femSpace.grad_psi
argsDict["u_test_ref"] = self.u[0].femSpace.psi
argsDict["u_grad_test_ref"] = self.u[0].femSpace.grad_psi
argsDict["mesh_trial_trace_ref"] = self.u[0].femSpace.elementMaps.psi_trace
argsDict["mesh_grad_trial_trace_ref"] = self.u[0].femSpace.elementMaps.grad_psi_trace
argsDict["dS_ref"] = self.elementBoundaryQuadratureWeights[('u', 0)]
argsDict["u_trial_trace_ref"] = self.u[0].femSpace.psi_trace
argsDict["u_grad_trial_trace_ref"] = self.u[0].femSpace.grad_psi_trace
def getMetricsAtEOS(self,u_exact):
import copy
"""
Calculate the element residuals and add in to the global residual
"""
degree_polynomial = 1
try:
degree_polynomial = self.u[0].femSpace.order
except:
pass
argsDict = cArgumentsDict.ArgumentsDict()
argsDict["mesh_trial_ref"] = self.u[0].femSpace.elementMaps.psi
argsDict["mesh_grad_trial_ref"] = self.u[0].femSpace.elementMaps.grad_psi
argsDict["mesh_dof"] = self.mesh.nodeArray
argsDict["mesh_l2g"] = self.mesh.elementNodesArray
argsDict["dV_ref"] = self.elementQuadratureWeights[('u',0)]
argsDict["u_trial_ref"] = self.u[0].femSpace.psi
argsDict["u_grad_trial_ref"] = self.u[0].femSpace.grad_psi
argsDict["u_test_ref"] = self.u[0].femSpace.psi
argsDict["nElements_global"] = self.mesh.nElements_global
argsDict["u_l2g"] = self.u[0].femSpace.dofMap.l2g
argsDict["elementDiameter"] = self.mesh.elementDiametersArray
argsDict["degree_polynomial"] = degree_polynomial
argsDict["epsFact_redist"] = self.coefficients.epsFact
argsDict["u_dof"] = self.u[0].dof
argsDict["u_exact"] = u_exact
argsDict["offset_u"] = self.offset[0]
#import np
import pdb
cfemIntegrals.zeroJacobian_CSR(self.nNonzerosInJacobian,jacobian)
# for now force time integration
useTimeIntegration = 1
if self.timeIntegration.__class__ == TimeIntegration.NoIntegration or not self.timeTerm:
useTimeIntegration = 0
if useTimeIntegration:
alpha_bdf = self.timeIntegration.alpha_bdf
beta_bdf = self.timeIntegration.beta_bdf
else:
alpha_bdf = self.timeIntegration.dt
beta_bdf = self.timeIntegration.m_last
argsDict = cArgumentsDict.ArgumentsDict()
argsDict["mesh_trial_ref"] = self.u[0].femSpace.elementMaps.psi
argsDict["mesh_grad_trial_ref"] = self.u[0].femSpace.elementMaps.grad_psi
argsDict["mesh_dof"] = self.mesh.nodeArray
argsDict["mesh_l2g"] = self.mesh.elementNodesArray
argsDict["x_ref"] = self.elementQuadraturePoints
argsDict["dV_ref"] = self.elementQuadratureWeights[('u', 0)]
argsDict["u_trial_ref"] = self.u[0].femSpace.psi
argsDict["u_grad_trial_ref"] = self.u[0].femSpace.grad_psi
argsDict["u_test_ref"] = self.u[0].femSpace.psi
argsDict["u_grad_test_ref"] = self.u[0].femSpace.grad_psi
argsDict["mesh_trial_trace_ref"] = self.u[0].femSpace.elementMaps.psi_trace
argsDict["mesh_grad_trial_trace_ref"] = self.u[0].femSpace.elementMaps.grad_psi_trace
argsDict["dS_ref"] = self.elementBoundaryQuadratureWeights[('u', 0)]
argsDict["u_trial_trace_ref"] = self.u[0].femSpace.psi_trace
argsDict["u_grad_trial_trace_ref"] = self.u[0].femSpace.grad_psi_trace
argsDict["u_test_trace_ref"] = self.u[0].femSpace.psi_trace
self.numericalFlux.setDirichletValues(self.ebqe)
# Flux boundary conditions
for ci, fbcObject in list(self.stressFluxBoundaryConditionsObjectsDict.items()):
for t, g in list(fbcObject.stressFluxBoundaryConditionsDict.items()):
self.ebqe[('stressFlux_bc', ci)][t[0], t[1]] = g(self.ebqe[('x')][t[0], t[1]], self.timeIntegration.t)
self.ebqe[('stressFlux_bc_flag', ci)][t[0], t[1]] = 1
r.fill(0.0)
self.elementResidual[0].fill(0.0)
self.elementResidual[1].fill(0.0)
if self.nSpace_global == 3:
self.elementResidual[2].fill(0.0)
if self.forceStrongConditions:
for cj in range(self.nc):
for dofN, g in list(self.dirichletConditionsForceDOF[cj].DOFBoundaryConditionsDict.items()):
self.u[cj].dof[dofN] = g(self.dirichletConditionsForceDOF[cj].DOFBoundaryPointDict[dofN], self.timeIntegration.t)
argsDict = cArgumentsDict.ArgumentsDict()
argsDict["mesh_trial_ref"] = self.u[0].femSpace.elementMaps.psi
argsDict["mesh_grad_trial_ref"] = self.u[0].femSpace.elementMaps.grad_psi
argsDict["mesh_dof"] = self.mesh.nodeArray
argsDict["mesh_l2g"] = self.mesh.elementNodesArray
argsDict["dV_ref"] = self.elementQuadratureWeights[('u', 0)]
argsDict["disp_trial_ref"] = self.u[0].femSpace.psi
argsDict["disp_grad_trial_ref"] = self.u[0].femSpace.grad_psi
argsDict["disp_test_ref"] = self.u[0].femSpace.psi
argsDict["disp_grad_test_ref"] = self.u[0].femSpace.grad_psi
argsDict["mesh_trial_trace_ref"] = self.u[0].femSpace.elementMaps.psi_trace
argsDict["mesh_grad_trial_trace_ref"] = self.u[0].femSpace.elementMaps.grad_psi_trace
argsDict["dS_ref"] = self.elementBoundaryQuadratureWeights[('u', 0)]
argsDict["disp_trial_trace_ref"] = self.u[0].femSpace.psi_trace
argsDict["disp_grad_trial_trace_ref"] = self.u[0].femSpace.grad_psi_trace
argsDict["disp_test_trace_ref"] = self.u[0].femSpace.psi_trace
argsDict["disp_grad_test_trace_ref"] = self.u[0].femSpace.grad_psi_trace
def getResidual(self, u, r):
import pdb
import copy
"""
Calculate the element residuals and add in to the global residual
"""
r.fill(0.0)
# Load the unknowns into the finite element dof
self.setUnknowns(u)
# no flux boundary conditions
argsDict = cArgumentsDict.ArgumentsDict()
argsDict["mesh_trial_ref"] = self.u[0].femSpace.elementMaps.psi
argsDict["mesh_grad_trial_ref"] = self.u[0].femSpace.elementMaps.grad_psi
argsDict["mesh_dof"] = self.mesh.nodeArray
argsDict["mesh_l2g"] = self.mesh.elementNodesArray
argsDict["x_ref"] = self.elementQuadraturePoints
argsDict["dV_ref"] = self.elementQuadratureWeights[('u', 0)]
argsDict["u_trial_ref"] = self.u[0].femSpace.psi
argsDict["u_grad_trial_ref"] = self.u[0].femSpace.grad_psi
argsDict["u_test_ref"] = self.u[0].femSpace.psi
argsDict["u_grad_test_ref"] = self.u[0].femSpace.grad_psi
argsDict["mesh_trial_trace_ref"] = self.u[0].femSpace.elementMaps.psi_trace
argsDict["mesh_grad_trial_trace_ref"] = self.u[0].femSpace.elementMaps.grad_psi_trace
argsDict["dS_ref"] = self.elementBoundaryQuadratureWeights[('u', 0)]
argsDict["u_trial_trace_ref"] = self.u[0].femSpace.psi_trace
argsDict["u_grad_trial_trace_ref"] = self.u[0].femSpace.grad_psi_trace
argsDict["u_test_trace_ref"] = self.u[0].femSpace.psi_trace