add check_residual

python/petsclsys.py

@@ -51,6 +51,7 @@ class LinearSystem :
         self.mat.zeroEntries()
         np.add.at(rhs.reshape([-1]),self.idx,localv)
         vmap = np.tile(self.elements*3,[1,3]) + np.repeat(np.arange(3,dtype=np.int32),3)[None,:]
+        # print(vmap)
         for e,m in zip(vmap,localm.reshape([-1,self.localsize**2])) :
             self.mat.setValues(e,e,m,PETSc.InsertMode.ADD)
         #for e,m in zip(self.elements,localm.reshape([-1,self.localsize**2])) :

python/time_integration.py

@@ -72,7 +72,7 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,max_nsub=None,afte
                 after_sub_iter(n_divide)
 
       
-def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, alpha=0.5,after_sub_iter=None,max_nsub=None) :  
+def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, alpha=0.5,after_sub_iter=None,max_nsub=None,check_residual_norm=-1) :  
     """Predictor-corrector scheme to solve fluid and grains.
 
     Keyword arguments:
@@ -101,7 +101,7 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
     # Predictor
     #
     # Compute the fluid solution and keep a copy of this solution and the forces applied by the fluid on the grains
-    fluid.implicit_euler(dt)
+    fluid.implicit_euler(dt,check_residual_norm)
     sf1 = np.copy(fluid.solution())
     f0 = np.copy(fluid.compute_node_force(dt))+external_particles_forces
     _advance_particles(particles,f0,dt,min_nsub,contact_tol,max_nsub=max_nsub)
@@ -113,7 +113,7 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
     #
     # Compute the fluid solution from the previous time-step one using the predicted solid velocities in the fluid-grain interaction force
     fluid.solution()[:,:] = sf0
-    fluid.implicit_euler(dt)
+    fluid.implicit_euler(dt,check_residual_norm)
     # Fluid solution is a weighted average of the predicted one and the corrected one.
     # The alpha parametre is the weight
     fluid.solution()[:,:] = (alpha*fluid.solution()+(1-alpha)*sf1)
@@ -129,7 +129,7 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
     fluid.move_particles(particles.position(), particles.velocity(), particles.contact_forces())
 
 
-def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False, after_sub_iter=None,max_nsub=None) :  
+def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False, after_sub_iter=None,max_nsub=None,check_residual_norm=-1) :  
     """Migflow solver: the solver type depends on the fluid and particles given as arguments.
 
     Keyword arguments:
@@ -153,7 +153,7 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
         if external_particles_forces.shape!=(particles.n_particles(),particles.dim()):
             raise ValueError("external_particles_forces must have shape (number of particles,dimension!")
     if (particles is None or particles.r() is None or particles.r().size == 0) and fluid is not None:
-        fluid.implicit_euler(dt)
+        fluid.implicit_euler(dt,check_residual_norm)
         # For two fluids flows, advance the concentration field with the fluid velocity.
         # The number of sub-time steps for the advection is automatically computed.
         if fluid.n_fluids() == 2 :
@@ -162,10 +162,10 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
         _advance_particles(particles,external_particles_forces,dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_nsub=max_nsub)
     if fluid is not None and particles is not None and particles.r() is not None and particles.r().size !=0:
         if fixed_grains:
-            fluid.implicit_euler(dt)
+            fluid.implicit_euler(dt,check_residual_norm)
             fluid.move_particles(particles.position(), particles.velocity(),-fluid.compute_node_force(dt))
         else:
-            predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces,after_sub_iter=after_sub_iter,max_nsub=max_nsub)
+            predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces,after_sub_iter=after_sub_iter,max_nsub=max_nsub,check_residual_norm=check_residual_norm)
 
 
 class FreeSurface():