heun

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 = 1e-3                                       # time step
+dt = 5e-3                                       # time step
 tEnd = 100                                      # final time
 
 # Geometrical parameters
@@ -124,24 +124,26 @@ tic = time.time()
 # COMPUTATION LOOP
 #
 while t < tEnd :
-    s = np.copy(fluid.solution())
-    vp = np.copy(p.velocity())
+    sf0 = np.copy(fluid.solution())
+    vp0 = np.copy(p.velocity())
     ### predictor
     # Fluid solver
     fluid.implicit_euler(dt)
+    sf1 = np.copy(fluid.solution())
     # Computation of the new velocities
-    forces = fluid.compute_node_force(dt)
-    p.velocity()[:,:] = p.velocity()+forces*dt/p.mass()
+    f0 = np.copy(fluid.compute_node_force(dt))
+    p.velocity()[:,:] = p.velocity()+f0*dt/p.mass()
     fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
     
     ### corrector
-    fluid.solution()[:,:] = s[:,:]
+    fluid.solution()[:,:] = sf0
     # Fluid solver
     fluid.implicit_euler(dt)
+    fluid.solution()[:,:] = (fluid.solution()+sf1)/2
     # Computation of the new velocities
-    forces = fluid.compute_node_force(dt)
-    p.velocity()[:,:] = vp[:,:]
-    vn = p.velocity() + forces*dt/p.mass()
+    f1 = fluid.compute_node_force(dt)
+    p.velocity()[:,:] = vp0
+    vn = vp0 + (f0+f1)/2*dt/p.mass()
 
     vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
     # Number of contact sub time step
@@ -150,7 +152,7 @@ while t < tEnd :
     # NLGS iterations
     for i in range(nsub) :
         tol = 1e-8                                #numerical tolerance for particles intersection
-        p.iterate(dt/nsub, forces, tol)
+        p.iterate(dt/nsub, (f0+f1)/2, tol)
 
     # Localisation of the particles in the fluid
     fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())

testcases/drop-2d/SimpleDrop/drop.py

@@ -78,7 +78,7 @@ if not os.path.isdir(outputdir) :
 # Physical parameters
 g = -9.81                                       # gravity
 rhop = 2450                                     # particles density
-r = 154e-6                                      # particles radii
+r = 154e-6/2.5                                      # particles radii
 compacity = 0.20                                # solid volume fraction in the drop
 rho = 1030                                      # fluid density
 nu = 1.17/rho                                   # kinematic viscosity
@@ -88,7 +88,7 @@ print("RHOP = %g, r = %g, RHO = %g" % (rhop,r,rho))
 
 # Numerical parameters
 outf = 1                                        # number of iterations between output files
-dt = 5e-2                                       # time step
+dt = 2e-2                                       # time step
 tEnd = 100                                      # final time
 
 #
@@ -108,7 +108,7 @@ tEnd1 = 0.2
 #
 # FLUID PROBLEM
 #
-fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],petsc_solver_type="-pc_type lu -ksp_max_it 50")
+fluid = fluid.FluidProblem(2,g,[nu*rho],[rho],drag_in_stab=1)
 # Set the mesh geometry for the fluid computation
 fluid.load_msh("mesh.msh")
 fluid.set_wall_boundary("Top", pressure=0)
@@ -123,30 +123,48 @@ fluid.solution()[:,2] = (fluid.coordinates()[:,1]-0.6)*rho*g
 fluid.export_vtk(outputdir,0,0)
 tic = time.time()
 
+tic = time.time()
 #
 # COMPUTATION LOOP
 #
-while t < tEnd : 
-    # Fluid solver
-    fluid.implicit_euler(dt,newton_max_it=100)
+while t < tEnd :
 
-    # Adaptation of the mesh.
-    if (ii%15==0 and ii != 0):
-       fluid.adapt_mesh(5e-3,8e-4,50000)
+    sf0 = np.copy(fluid.solution())
+    vp0 = np.copy(p.velocity())
+    ### predictor
+    # Fluid solver
+    fluid.implicit_euler(dt)
+    sf1 = np.copy(fluid.solution())
+    # Computation of the new velocities
+    f0 = np.copy(fluid.compute_node_force(dt))
+    p.velocity()[:,:] = p.velocity()+f0*dt/p.mass()
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
     
+    ### corrector
+    fluid.solution()[:,:] = sf0
+    # Fluid solver
+    fluid.implicit_euler(dt)
+    fluid.solution()[:,:] = (fluid.solution()+sf1)/2
     # Computation of the new velocities
-    forces = fluid.compute_node_force(dt)
-    vn = p.velocity() + forces * dt / p.mass()
+    f1 = fluid.compute_node_force(dt)
+    p.velocity()[:,:] = vp0
+    vn = vp0 + (f0+f1)/2*dt/p.mass()
+
     vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
-    # Number of sub time step
-    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
+    # Number of contact sub time step
+    nsub = max(1, int(np.ceil((vmax*dt*4)/min(p.r()))))
     print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
     # NLGS iterations
     for i in range(nsub) :
-        p.iterate(dt/nsub, forces)  
+        tol = 1e-8                                #numerical tolerance for particles intersection
+        p.iterate(dt/nsub, (f0+f1)/2, tol)
 
-    t += dt
+    if (ii%15==0 and ii != 0):
+      fluid.adapt_mesh(2e-2,1e-3,10000)
+
+    # Localisation of the particles in the fluid
     fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
+    t += dt
 
     # Output files writting
     if ii %outf == 0 :