predictor-corrector for drag in stab

testcases/depot-2d/depotPredCorr.py

@@ -75,7 +75,7 @@ nu = 1e-6                                       # kinematic viscosity
 
 # Numerical parameters
 outf = 1                                        # number of iterations between output files
-dt = 0.25e-2                                       # time step
+dt = 1e-3                                       # time step
 tEnd = 100                                      # final time
 
 # Geometrical parameters
@@ -116,11 +116,10 @@ fluid.set_strong_boundary("Bottom",0,0)
 fluid.set_strong_boundary("Bottom",1,0)
 fluid.solution()[:,2] = (fluid.coordinates()[:,1]-0.6)*rho*g
 # Set location of the particles in the mesh and compute the porosity in each computation cell
-fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces(), init=True)
 fluid.export_vtk(outputdir,0,0)
 
 tic = time.time()
-fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
 #
 # COMPUTATION LOOP
 #
@@ -132,8 +131,9 @@ while t < tEnd :
     fluid.implicit_euler(dt)
     # Computation of the new velocities
     forces = fluid.compute_node_force(dt)
-    p.velocity()[:,:] = p.velocity() + forces*dt/p.mass()
-
+    p.velocity()[:,:] = p.velocity()+forces*dt/p.mass()
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
+    
     ### corrector
     fluid.solution()[:,:] = s[:,:]
     # Fluid solver