diff --git a/proteus/richards/Richards.h b/proteus/richards/Richards.h index 1adaaefeef..497ba6c03b 100644 --- a/proteus/richards/Richards.h +++ b/proteus/richards/Richards.h @@ -32,12 +32,15 @@ namespace proteus } namespace proteus + { class Richards_base { //The base class defining the interface public: - virtual ~Richards_base(){} + virtual ~Richards_base(){ + double anb_seepage_flux =1e-16; + } virtual void calculateResidual(arguments_dict& args)=0; virtual void calculateJacobian(arguments_dict& args)=0; virtual void invert(arguments_dict& args)=0; @@ -113,6 +116,7 @@ namespace proteus psiC = -u; m_vg = 1.0 - 1.0 / n_vg; thetaS = thetaR + thetaSR; + //std::cout<< "Thetas"< 0.0) { pcBar = alpha * psiC; @@ -172,7 +176,7 @@ namespace proteus } } } - + inline void evaluateInverseCoefficients(const int rowptr[nSpace], const int colind[nnz], @@ -183,7 +187,7 @@ namespace proteus const double n_vg, const double thetaR, const double thetaSR, - const double KWs[nnz], + const double KWs[nnz], double& u, const double& m, const double& dm, @@ -191,7 +195,7 @@ namespace proteus const double df[nSpace], const double a[nnz], const double da[nnz]) - { + { double psiC, pcBar,pcBar_n, sBar, @@ -222,7 +226,6 @@ namespace proteus <<"psiC "<& csrRowIndeces_DofLoops, + const xt::pyarray& csrColumnOffsets_DofLoops, + const xt::pyarray& uLow, + const xt::pyarray& dt_times_fH_minus_fL, + double dt, + double mi +) { + double Pposi = 0.0, Pnegi = 0.0; + double mini = 0.0, maxi = 1.0; + int ij = csrRowIndeces_DofLoops(i); // assuming using xtensor array + + for (int offset = csrRowIndeces_DofLoops(i); offset < csrRowIndeces_DofLoops(i + 1); ++offset) { + int j = csrColumnOffsets_DofLoops(offset); + // Compute P vectors + Pposi += dt_times_fH_minus_fL(offset) * (dt_times_fH_minus_fL(offset) > 0 ? 1.0 : 0.0); + Pnegi += dt_times_fH_minus_fL(offset) * (dt_times_fH_minus_fL(offset) < 0 ? 1.0 : 0.0); + } + + double Qposi = mi * (maxi - uLow(i)); + double Qnegi = mi * (mini - uLow(i)); + double Rposi = (Pposi == 0.0) ? 1.0 : std::fmin(1.0, Qposi / Pposi); + double Rnegi = (Pnegi == 0.0) ? 1.0 : std::fmin(1.0, Qnegi / Pnegi); + + double ith_limited_flux_correction = 0.0; + for (int offset = csrRowIndeces_DofLoops(i); offset < csrRowIndeces_DofLoops(i + 1); ++offset) { + int j = csrColumnOffsets_DofLoops(offset); + double Lij = dt_times_fH_minus_fL(offset) > 0 ? std::fmin(Rposi, uLow(j)) : std::fmin(Rnegi, uLow(j)); + ith_limited_flux_correction += Lij * dt_times_fH_minus_fL(offset); + } + + if (std::isnan(ith_limited_flux_correction)) { + ith_limited_flux_correction = 0.0; + } + + return ith_limited_flux_correction; +} + + + void calculateResidual(arguments_dict& args) { xt::pyarray& mesh_trial_ref = args.array("mesh_trial_ref"); @@ -501,11 +562,26 @@ namespace proteus xt::pyarray& quantDOFs = args.array("quantDOFs"); xt::pyarray& sLow = args.array("sLow"); xt::pyarray& sn = args.array("sn"); - + assert(a_rowptr.data()[nSpace] == nnz); assert(a_rowptr.data()[nSpace] == nSpace); //cek should this be read in? double Ct_sge = 4.0; + //For flux Calculation + //double anb_seepage_flux=0.0; + //anb_seepage_flux = args["anb_seepage_flux"]; + //double anb_seepage_flux = args.scalar("anb_seepage_flux"); + //anb_seepage_flux=0.0; + + xt::pyarray& anb_seepage_flux_n = args.array("anb_seepage_flux_n"); + + + //double anb_seepage_flux=0.0; + double & anb_seepage_flux (args.scalar("anb_seepage_flux")) ; + //double anb_seepage_flux = args.scalar("anb_seepage_flux_n"); + + //double anb_seepage_flux=0.0; + anb_seepage_flux=0.0; //loop over elements to compute volume integrals and load them into element and global residual // @@ -648,8 +724,6 @@ namespace proteus /* // */ /* //calculate shock capturing diffusion */ /* // */ - - /* ck.calculateNumericalDiffusion(shockCapturingDiffusion,elementDiameter[eN],pdeResidual_u,grad_u,numDiff0); */ /* //ck.calculateNumericalDiffusion(shockCapturingDiffusion,G,pdeResidual_u,grad_u_old,numDiff1); */ /* ck.calculateNumericalDiffusion(shockCapturingDiffusion,sc_uref, sc_alpha,G,G_dd_G,pdeResidual_u,grad_u,numDiff1); */ @@ -718,6 +792,7 @@ namespace proteus da_ext[nnz], as_ext[nnz], flux_ext=0.0, + //anb_seepage_flux=0.0, // for flux calculation bc_u_ext=0.0, bc_grad_u_ext[nSpace], bc_m_ext=0.0, @@ -852,17 +927,32 @@ namespace proteus ebqe_penalty_ext.data()[ebNE_kb],// penalty, flux_ext); ebqe_flux.data()[ebNE_kb] = flux_ext; - ebqe_u.data()[ebNE_kb] = u_ext; + + anb_seepage_flux= seepagefluxcalculator(anb_seepage_flux, + isSeepageFace.data()[ebNE], + dS, + flux_ext); + //std::cout<<"The seepage flux is "<0) + { + std::cout<<"The seepage flux is "<0 ? Rposi : Rpos[j]); + Lij = (Fluxij>0 ? Rposi : Rpos[j]); // compute limited flux ith_Limiter_times_FluxCorrectionMatrix += Lij*Fluxij; - + // update limited flux limitedFlux.data()[ij] = Lij*Fluxij; - + //update FluxMatrix - FluxMatrix.data()[ij] = Fluxij; - + FluxMatrix.data()[ij] = Fluxij; + //update ij ij+=1; } @@ -1575,7 +1669,7 @@ namespace proteus double Qnegi = mi*(mini-solLim.data()[i]); // compute R vectors // Rpos[i] = ((Pposi==0) ? 1. : fmin(1.0,Qposi/Pposi)); - Rneg[i] = ((Pnegi==0) ? 1. : fmin(1.0,Qnegi/Pnegi)); + Rneg[i] = ((Pnegi==0) ? 1. : fmin(1.0,Qnegi/Pnegi)); } // COMPUTE LIMITERS // @@ -1623,6 +1717,7 @@ namespace proteus xt::pyarray& mesh_grad_trial_trace_ref = args.array("mesh_grad_trial_trace_ref"); xt::pyarray& dS_ref = args.array("dS_ref"); xt::pyarray& u_trial_trace_ref = args.array("u_trial_trace_ref"); + xt::pyarray& u_grad_trial_trace_ref = args.array("u_grad_trial_trace_ref"); xt::pyarray& u_test_trace_ref = args.array("u_test_trace_ref"); xt::pyarray& u_grad_test_trace_ref = args.array("u_grad_test_trace_ref"); @@ -1717,6 +1812,18 @@ namespace proteus xt::pyarray& quantDOFs = args.array("quantDOFs"); xt::pyarray& sLow = args.array("sLow"); xt::pyarray& sn = args.array("sn"); + + xt::pyarray& anb_seepage_flux_n = args.array("anb_seepage_flux_n"); + + + //double anb_seepage_flux=0.0; + double & anb_seepage_flux (args.scalar("anb_seepage_flux")) ; + //double anb_seepage_flux = args.scalar("anb_seepage_flux_n"); + + //double anb_seepage_flux=0.0; + anb_seepage_flux=0.0; + + double Rpos[numDOFs], Rneg[numDOFs]; //double FluxCorrectionMatrix[NNZ]; // NOTE: This function follows a different (but equivalent) implementation of the smoothness based indicator than NCLS.h @@ -1727,6 +1834,8 @@ namespace proteus std::valarray u_free_dof(numDOFs); std::valarray u_free_dof_old(numDOFs); std::valarray ML2(numDOFs); + + for(int eN=0;eN= 0 && isFluxBoundary_u[ebNE_kb] != 1 ) //outflow. This is handled via the transport matrices. Then flux_ext=0 and dflux_ext!=0 */ - /* { */ - /* dflux_ext = flow; */ - /* flux_ext = 0; */ - /* // save external u */ - /* ebqe_u[ebNE_kb] = u_ext; */ - /* } */ - /* else // inflow. This is handled via the boundary integral. Then flux_ext!=0 and dflux_ext=0 */ - /* { */ - /* dflux_ext = 0; */ - /* // save external u */ - /* ebqe_u[ebNE_kb] = isDOFBoundary_u[ebNE_kb]*ebqe_bc_u_ext[ebNE_kb]+(1-isDOFBoundary_u[ebNE_kb])*u_ext; */ - /* if (isDOFBoundary_u[ebNE_kb] == 1) */ - /* flux_ext = ebqe_bc_u_ext[ebNE_kb]*flow; */ - /* else if (isFluxBoundary_u[ebNE_kb] == 1) */ - /* flux_ext = ebqe_bc_flux_u_ext[ebNE_kb]; */ - /* else */ - /* { */ - /* std::cout<<"warning: VOF open boundary with no external trace, setting to zero for inflow"<0) + // { + // std::cout<<"The seepage flux is "<("anb_seepage_flux_n"); + // //anb_seepage_flux_n=0.0; + + // //anb_seepage_flux_n= anb_seepage_flux; + + + + // //if (isSeepageFace) + // //{ + // // anb_seepage_flux+= dS* flux_ext; + // // std::cout<<"The seepage flux is "<0){ + // // std::cout<<"The seepage flux is "<= 0 && isFluxBoundary_u[ebNE_kb] != 1 ) //outflow. This is handled via the transport matrices. Then flux_ext=0 and dflux_ext!=0 + // // { + // // dflux_ext = flow; + // // flux_ext = 0; + // // // save external u + // // ebqe_u[ebNE_kb] = u_ext; + // // } + // //else // inflow. This is handled via the boundary integral. Then flux_ext!=0 and dflux_ext=0 + // //{ + // //dflux_ext = 0; + // // save external u + // //ebqe_u[ebNE_kb] = isDOFBoundary_u[ebNE_kb]*ebqe_bc_u_ext[ebNE_kb]+(1-isDOFBoundary_u[ebNE_kb])*u_ext; + // //if (isDOFBoundary_u[ebNE_kb] == 1) + // // flux_ext = ebqe_bc_u_ext[ebNE_kb]*flow; + // //else if (isFluxBoundary_u[ebNE_kb] == 1) + // //flux_ext = ebqe_bc_flux_u_ext[ebNE_kb]; + // // flux_ext = ebqe_bc_u_ext[ebNE_kb]; + // //else + // // { + // // std::cout<<"warning: VOF open boundary with no external trace, setting to zero for inflow"< mMax) - { - std::cout<<"mass out of bounds "< mMax) + { + std::cout<<"mass out of bounds "< 1.0 || Rposi < 0.0) + //if (Rposi > 1.0 || Rposi < 0.0) //std::cout << "Rposi: " << Rposi << std::endl; - //if (Rnegi > 1.0 || Rnegi < 0.0) - //std::cout << "Rnegi: " << Rnegi << std::endl; - + //if (Rnegi > 1.0 || Rnegi < 0.0) + //std::cout << "Rnegi: " << Rnegi << std::endl; + // LOOP OVER THE SPARSITY PATTERN (j-LOOP)// // for (int offset=csrRowIndeces_DofLoops.data()[i]; // offset& a_rowptr = args.array("a_rowptr"); @@ -2926,6 +3202,8 @@ namespace proteus }//computeMassMatrix };//Richards + + inline Richards_base* newRichards(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, @@ -2936,11 +3214,11 @@ namespace proteus { if (nSpaceIn == 1) return proteus::chooseAndAllocateDiscretization1D(nSpaceIn, - nQuadraturePoints_elementIn, - nDOF_mesh_trial_elementIn, - nDOF_trial_elementIn, - nDOF_test_elementIn, - nQuadraturePoints_elementBoundaryIn, + nQuadraturePoints_elementIn, + nDOF_mesh_trial_elementIn, + nDOF_trial_elementIn, + nDOF_test_elementIn, + nQuadraturePoints_elementBoundaryIn, CompKernelFlag); else if (nSpaceIn == 2) return proteus::chooseAndAllocateDiscretization2D(nSpaceIn, @@ -2962,5 +3240,5 @@ namespace proteus CompKernelFlag); } } -}//proteus -#endif +};//proteus +#endif \ No newline at end of file diff --git a/proteus/richards/Richards.py b/proteus/richards/Richards.py index 841fcd0884..0cdd506784 100644 --- a/proteus/richards/Richards.py +++ b/proteus/richards/Richards.py @@ -12,6 +12,8 @@ from proteus.LinearAlgebraTools import SparseMat from proteus import TimeIntegration from proteus.NonlinearSolvers import ExplicitLumpedMassMatrixForRichards +from proteus.NonlinearSolvers import Newton + class ThetaScheme(TimeIntegration.BackwardEuler): def __init__(self,transport,integrateInterpolationPoints=False): @@ -51,6 +53,8 @@ def __init__(self, transport, timeOrder=1, runCFL=0.1, integrateInterpolationPoi self.u_dof_stage[ci] = [] for k in range(self.nStages + 1): self.u_dof_stage[ci].append(transport.u[ci].dof.copy()) + #print() + # def set_dt(self, DTSET): # self.dt = DTSET # don't update t @@ -63,6 +67,8 @@ def choose_dt(self): self.dtLast = self.dt self.t = self.tLast + self.dt self.substeps = [self.t for i in range(self.nStages)] # Manuel is ignoring different time step levels for now + + def initialize_dt(self, t0, tOut, q): """ @@ -134,6 +140,7 @@ def updateStage(self): for ci in range(self.nc): self.u_dof_stage[ci][self.lstage][:] = self.transport.u[ci].dof[:] self.transport.u_dof_old[:] = self.transport.u[ci].dof + def initializeTimeHistory(self, resetFromDOF=True): """ @@ -199,6 +206,8 @@ def setFromOptions(self, nOptions): setattr(self, flag, val) if flag == 'timeOrder': self.resetOrder(self.timeOrder) + + class Coefficients(proteus.TransportCoefficients.TC_base): """ @@ -223,8 +232,7 @@ def __init__(self, ENTROPY_TYPE=2, # logarithmic LUMPED_MASS_MATRIX=False, MONOLITHIC=True, - FCT=False, - forceStrongBoundaryConditions=False, + FCT=True, num_fct_iter=1, # FOR ENTROPY VISCOSITY cE=1.0, @@ -233,7 +241,9 @@ def __init__(self, # FOR ARTIFICIAL COMPRESSION cK=1.0, # OUTPUT quantDOFs - outputQuantDOFs=False): + outputQuantDOFs=False, ): + self.anb_seepage_flux= 0.00 + #self.anb_seepage_flux_n =0.0 variableNames=['pressure_head'] nc=1 mass={0:{0:'nonlinear'}} @@ -294,7 +304,7 @@ def __init__(self, self.uL = uL self.uR = uR self.cK = cK - self.forceStrongConditions = forceStrongBoundaryConditions + self.forceStrongConditions = True self.cE = cE self.outputQuantDOFs = outputQuantDOFs TC_base.__init__(self, @@ -320,12 +330,15 @@ def initializeMesh(self,mesh): #mwf missing ebNE-->ebN? ebN = mesh.exteriorElementBoundariesArray[ebNE] #print "eb flag",mesh.elementBoundaryMaterialTypes[ebN] + #print self.getSeepageFace(mesh.elementBoundaryMaterialTypes[ebN]) self.isSeepageFace[ebNE] = self.getSeepageFace(mesh.elementBoundaryMaterialTypes[ebN]) - #print self.isSeepageFace + #print (self.isSeepageFace) def initializeElementQuadrature(self,t,cq): self.materialTypes_q = self.elementMaterialTypes self.q_shape = cq[('u',0)].shape + #self.anb_seepage_flux= anb_seepage_flux + #print("The seepage is ", anb_seepage_flux) # cq['Ks'] = np.zeros(self.q_shape,'d') # for k in range(self.q_shape[1]): # cq['Ks'][:,k] = self.Ksw_types[self.elementMaterialTypes,0] @@ -342,6 +355,8 @@ def initializeGlobalExteriorElementBoundaryQuadrature(self,t,cebqe): self.ebqe_shape = cebqe[('u',0)].shape self.ebqe[('vol_frac',0)] = np.zeros(self.ebqe_shape,'d') # + + def evaluate(self,t,c): if c[('u',0)].shape == self.q_shape: materialTypes = self.materialTypes_q @@ -402,11 +417,13 @@ def evaluate(self,t,c): np.isnan(c[('dm',0,0)]).any()): import pdb pdb.set_trace() - -# #mwf debug -# if c[('u',0)].shape == self.q_shape: -# c[('visPerm',0)]=c[('a',0,0)][:,:,0,0] - + + def postStep(self, t, firstStep=False): + # #anb_seepage_flux_n[:]= self.anb_seepage_flux + with open('seepage_stab_0', "a") as f: + # f.write("\n Time"+ ",\t" +"Seepage\n") + f.write(repr(t)+ ",\t")# +repr(np.sum(self.LevelModel.anb_seepage_flux_n))) + class LevelModel(proteus.Transport.OneLevelTransport): nCalls=0 def __init__(self, @@ -538,6 +555,7 @@ def __init__(self, for I in self.coefficients.elementIntegralKeys: if I in elementQuadrature: elementQuadratureDict[I] = elementQuadrature[I] + else: elementQuadratureDict[I] = elementQuadrature['default'] else: @@ -557,6 +575,7 @@ def __init__(self, if elemQuadIsDict: if ('numDiff',ci,ci) in elementQuadrature: elementQuadratureDict[('numDiff',ci,ci)] = elementQuadrature[('numDiff',ci,ci)] + else: elementQuadratureDict[('numDiff',ci,ci)] = elementQuadrature['default'] else: @@ -716,6 +735,17 @@ def __init__(self, assert cond, "Use lumped mass matrix just with: STABILIZATION_TYPE=2 (smoothness based stab.)" cond = 'levelNonlinearSolver' in dir(options) and options.levelNonlinearSolver == ExplicitLumpedMassMatrixForRichards assert cond, "Use levelNonlinearSolver=ExplicitLumpedMassMatrixForRichards when the mass matrix is lumped" + + ################################################################# + ####################ARNOB_FCT_EDIT############################### + ################################################################# + if not self.coefficients.LUMPED_MASS_MATRIX and self.coefficients.STABILIZATION_TYPE == 2: + cond = 'levelNonlinearSolver' in dir(options) and options.levelNonlinearSolver == Newton + + + + + if self.coefficients.FCT == True: cond = self.coefficients.STABILIZATION_TYPE > 0, "Use FCT just with STABILIZATION_TYPE>0; i.e., edge based stabilization" # END OF ASSERTS @@ -739,6 +769,7 @@ def __init__(self, self.sLow = np.zeros(self.u[0].dof.shape, 'd') self.sHigh = np.zeros(self.u[0].dof.shape, 'd') self.sn = np.zeros(self.u[0].dof.shape, 'd') + self.anb_seepage_flux_n = np.zeros(self.u[0].dof.shape, 'd') comm = Comm.get() self.comm=comm if comm.size() > 1: @@ -749,6 +780,8 @@ def __init__(self, self.stride = [self.nc for ci in range(self.nc)] # logEvent(memory("stride+offset","OneLevelTransport"),level=4) + + if numericalFluxType != None: if options is None or options.periodicDirichletConditions is None: self.numericalFlux = numericalFluxType(self, @@ -821,7 +854,7 @@ def __init__(self, self.MOVING_DOMAIN=0.0 if self.mesh.nodeVelocityArray is None: self.mesh.nodeVelocityArray = np.zeros(self.mesh.nodeArray.shape,'d') - self.forceStrongConditions=self.coefficients.forceStrongConditions + self.forceStrongConditions=True self.dirichletConditionsForceDOF = {} if self.forceStrongConditions: for cj in range(self.nc): @@ -860,6 +893,7 @@ def FCTStep(self): argsDict["max_s_bc"] = np.ones_like(self.min_s_bc)*self.coefficients.rho*(self.coefficients.thetaR_types[0] + self.coefficients.thetaSR_types[0]) #cek hack argsDict["LUMPED_MASS_MATRIX"] = self.coefficients.LUMPED_MASS_MATRIX argsDict["MONOLITHIC"] =0#cek hack self.coefficients.MONOLITHIC + argsDict["anb_seepage_flux_n"]= self.anb_seepage_flux_n #pdb.set_trace() self.richards.FCTStep(argsDict) @@ -1177,8 +1211,7 @@ def getResidual(self,u,r): pass argsDict = cArgumentsDict.ArgumentsDict() - if self.forceStrongConditions: - argsDict["bc_mask"] = self.bc_mask + argsDict["bc_mask"] = self.bc_mask argsDict["dt"] = self.timeIntegration.dt argsDict["Theta"] = 1.0 argsDict["mesh_trial_ref"] = self.u[0].femSpace.elementMaps.psi @@ -1292,12 +1325,36 @@ def getResidual(self,u,r): argsDict["quantDOFs"] = self.quantDOFs argsDict["sLow"] = self.sLow argsDict["sn"] = self.sn + argsDict["anb_seepage_flux_n"]= self.anb_seepage_flux_n +###################################################################################### + argsDict["pn"] = self.u[0].dof + argsDict["solH"] = self.sHigh + + rowptr, colind, MassMatrix = self.MC_global.getCSRrepresentation() + argsDict["MassMatrix"] = MassMatrix + argsDict["solH"] = self.sHigh + +###################################################################################### + argsDict["anb_seepage_flux"] = self.coefficients.anb_seepage_flux + #print(anb_seepage_flux) + #argsDict["anb_seepage_flux_n"] = self.coefficients.anb_seepage_flux_n + #if np.sum(anb_seepage_flux_n)>0: + + #logEvent("Hi, this is Arnob", self.anb_seepage_flux_n[0]) + #print("Seepage Flux from Python file", np.sum(self.anb_seepage_flux_n)) + seepage_text_variable= np.sum(self.anb_seepage_flux_n) + + with open('seepage_stab_0',"a" ) as f: + #f.write("\n Time"+ ",\t" +"Seepage\n") + #f.write(repr(self.coefficients.t)+ ",\t" +repr(seepage_text_variable), "\n") + f.write(repr(seepage_text_variable)+ "\n") if (self.coefficients.STABILIZATION_TYPE == 0): # SUPG self.calculateResidual = self.richards.calculateResidual self.calculateJacobian = self.richards.calculateJacobian else: self.calculateResidual = self.richards.calculateResidual_entropy_viscosity self.calculateJacobian = self.richards.calculateMassMatrix + if self.delta_x_ij is None: self.delta_x_ij = -np.ones((self.nNonzerosInJacobian*3,),'d') self.calculateResidual(argsDict) @@ -1316,6 +1373,8 @@ def getResidual(self,u,r): if self.globalResidualDummy is None: self.globalResidualDummy = np.zeros(r.shape,'d') + + def invert(self,u,r): #u=s #r=p @@ -1468,6 +1527,10 @@ def invert(self,u,r): argsDict["quantDOFs"] = self.quantDOFs argsDict["sLow"] = self.sLow argsDict["sn"] = self.sn + #Arnob trying to print flux + argsDict["anb_seepage_flux"] = self.coefficients.anb_seepage_flux + + self.richards.invert(argsDict) # self.timeIntegration.dt, # self.u[0].femSpace.elementMaps.psi, @@ -1591,6 +1654,10 @@ def invert(self,u,r): #self.nonlinear_function_evaluations += 1 #if self.globalResidualDummy is None: # self.globalResidualDummy = numpy.zeros(r.shape,'d') + def postStep(self, t, firstStep=False): + with open('seepage_flux_nnnn', "a") as f: + f.write("\n Time"+ ",\t" +"Seepage\n") + f.write(repr(t)+ ",\t" +repr(self.coefficients.anb_seepage_flux)) def getJacobian(self,jacobian): if (self.coefficients.STABILIZATION_TYPE == 0): # SUPG cfemIntegrals.zeroJacobian_CSR(self.nNonzerosInJacobian, @@ -1666,6 +1733,8 @@ def getJacobian(self,jacobian): argsDict["ebqe_bc_flux_ext"] = self.ebqe[('advectiveFlux_bc',0)] argsDict["csrColumnOffsets_eb_u_u"] = self.csrColumnOffsets_eb[(0,0)] argsDict["LUMPED_MASS_MATRIX"] = self.coefficients.LUMPED_MASS_MATRIX + argsDict["anb_seepage_flux"] = self.coefficients.anb_seepage_flux + self.calculateJacobian(argsDict) if self.forceStrongConditions: for dofN in list(self.dirichletConditionsForceDOF[0].DOFBoundaryConditionsDict.keys()): @@ -1737,9 +1806,27 @@ def calculateExteriorElementBoundaryQuadrature(self): getDiffusiveFluxBoundaryConditions=self.diffusiveFluxBoundaryConditionsSetterDictDict[cj])) for cj in list(self.advectiveFluxBoundaryConditionsSetterDict.keys())]) self.coefficients.initializeGlobalExteriorElementBoundaryQuadrature(self.timeIntegration.t,self.ebqe) + #argsDict = cArgumentsDict.ArgumentsDict() + #argsDict["anb_seepage_flux"] = self.coefficients.anb_seepage_flux + #print("Hi", self.coefficients.anb_seepage_flux) + + #print("The seepage is ", anb_seepage_flux) def estimate_mt(self): pass def calculateSolutionAtQuadrature(self): pass def calculateAuxiliaryQuantitiesAfterStep(self): pass + def postStep(self, t, firstStep=False): + with open('seepage_flux_nnnnk', "a") as f: + f.write("\n Time"+ ",\t" +"Seepage\n") + f.write(repr(t)+ ",\t" +repr(self.coefficients.anb_seepage_flux)) + + + +#argsDict["anb_seepage_flux"] = self.coefficients.anb_seepage_flux +#anb_seepage_flux= self.coefficients.anb_seepage_flux +#print("Hi",anb_seepage_flux) + + +#print("Hello from the python file", self.coefficients.anb_seepage_flux)