testcase for new adaptation

testcases/drop-2d/dropNewAdapt/adapt/mesh.geo

+L = 0.025;
+H = 0.5;
+y = 0;
+lc = 0.1;
+
+Point(1) = {-L, H, 0, lc};
+Point(2) = {-L, 0, 0, lc};
+Point(3) = {L, 0, 0, lc};
+Point(4) = {L, H, 0, lc};
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Line Loop(1) = {1:4};
+
+Plane Surface(1) = {1};
+Physical Line("Box") = {1,2,3};
+Physical Line("Top") = {4};
+Physical Surface("Domain") = {1};
+Physical Point("PtFix") = {1};
+Mesh.Algorithm = 5;
+
+Merge "lc.pos";
+Field[1] = PostView;
+Field[1].IView = 0;
+Background Field = 1;
+
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
+Mesh.CharacteristicLengthFromPoints = 0;

testcases/drop-2d/dropNewAdapt/drop-bell.py

+#!/usr/bin/env python
+from marblesbag import fluid as fluid
+from marblesbag import scontact2
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+
+
+def genInitialPosition(filename, r, rout, rhop) :
+    D = 3*r
+    p = scontact2.ParticleProblem()
+    p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
+    x = np.arange(rout, -rout, -D)
+    x, y = np.meshgrid(x, x)
+    R2 = x**2 + y**2
+    keep = R2 < (rout - r)**2
+    x = x[keep]
+    y = y[keep]
+    for i in range(x.shape[0]) :
+        if y[i]<0 :
+            p.add_particle((x[i]-D/2.1+2*random.uniform(0,D/5), y[i]-D/2.1+2*random.uniform(0,D/5)), r, r**2 * np.pi * rhop);
+    
+    x = np.arange(rout, -rout, -D)
+    y = np.arange(2*r, 4*rout, D)
+    x, y = np.meshgrid(x, y)
+    x = x.flat
+    y = y.flat
+    T = 1.5*rout
+    x1 = -4/6*rout
+    x2 = 4/6*rout
+    for i in range(len(x)) :
+        if y[i] < -T/x1*x[i]+T and y[i] < -T/x2*x[i]+T :
+            p.add_particle((x[i]-D/2.1+2*random.uniform(0,D/5), y[i]-D/2.1+2*random.uniform(0,D/5)), r, r**2 * np.pi * rhop);
+    p.position()[:, 1] += 0.45
+    print("N=%e\n"%len(p.position()[:, 1]))
+    p.write(filename)
+
+
+outputdir = "outputBell"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+t = 0
+ii = 0
+
+r=25e-6/(5)
+rout = 1e-3
+p = scontact2.ParticleProblem()
+#R = np.random.uniform(45e-06, 90e-06, len(x))
+
+#physical parameters
+g =  -9.81
+rho = 1200
+rhop = 2400
+nu = 1e-4
+V = 0.5 # todo : estimate V base on limit velocity
+print('V',V)
+tEnd = 10
+
+#numerical parameters
+lcmin = 0.0001 # approx r*100 but should match the mesh size
+dt = 5e-4
+alpha = 1e-2
+epsilon = alpha*lcmin**2 /nu#2e-2*c*h*2#2e-2*c*h # ?? not sure ??1e-5
+print('epsilon',epsilon)
+
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
+p.write(outputdir+"/part-00000")
+
+genInitialPosition(outputdir + "/part-00000", r, rout, rhop)
+
+p = scontact2.ParticleProblem(outputdir+"/part-00000")
+
+print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
+print("RHOP = %g" % rhop)
+outf = 50
+outf1 = 100000
+
+strong_boundaries = [("Top",2,0.),("Top",1,0.),("Top",0,0.),("Box",1,0.),("Box",0,0.)]
+fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
+
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+fluid.export(outputdir,0,0)
+
+tic = time.clock()
+
+mu = nu*rho
+porosity = len(p.velocity()[:,1])*r**2/(rout**2)
+p.velocity()[:,1] = 2./9. *porosity* rout**2*g*(rhop-rho)/mu
+print(p.velocity()[:,1])
+while ii<1000:
+    if (ii%100==0 and ii != 0):
+       fluid.adapt_mesh(0.01,50,20,500,1e-4,7500)
+    forces = fluid.compute_node_force(dt)
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    fluid.implicit_euler(dt)
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        p.write_vtk(outputdir, ioutput, t)
+        p.write(outputdir+"/part-%05d" % ioutput)
+        p.write_boundary_vtk(outputdir, ioutput, t)
+        fluid.export(outputdir, t, ioutput)
+    ii += 1
+
+
+forces = g*p.mass()
+while t < tEnd : 
+    if (ii%100==0 and ii != 0):
+       fluid.adapt_mesh(0.01,50,20,500,1e-4,7500)
+    forces = fluid.compute_node_force(dt)
+    vn = p.velocity() + forces * dt / p.mass()
+    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
+    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(10*p.r()))))
+    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
+    for i in range(nsub) :
+        p.iterate(dt/nsub, forces)  
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    fluid.implicit_euler(dt)
+    t += dt
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        p.write_vtk(outputdir, ioutput, t)
+        p.write(outputdir+"/part-%05d" % ioutput)
+        p.write_boundary_vtk(outputdir, ioutput, t)
+        fluid.export(outputdir, t, ioutput)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))

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

+#!/usr/bin/env python
+from marblesbag import fluid as fluid
+from marblesbag import scontact2
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+
+def genInitialPosition(filename, r, rout, rhop) :
+    p = scontact2.ParticleProblem()
+    p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
+    x = np.arange(rout, -rout, -2.5*r)
+    x, y = np.meshgrid(x, x)
+    R2 = x**2 + y**2
+    keep = R2 < (rout - r)**2
+    x = x[keep]
+    y = y[keep]
+    for i in range(x.shape[0]) :
+        p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
+    
+    p.position()[:, 1] += 0.45
+    p.write(filename)
+
+
+outputdir = "outputNew"
+outputdir1 = "outputNew/remesh"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+if not os.path.isdir(outputdir1) :
+    os.makedirs(outputdir1)
+
+t = 0
+ii = 0
+
+r=25e-6/2
+p = scontact2.ParticleProblem()
+#R = np.random.uniform(45e-06, 90e-06, len(x))
+
+#physical parameters
+g =  -9.81
+rho = 1200
+rhop = 2400
+nu = 1e-4
+V = 0.5 # todo : estimate V base on limit velocity
+print('V',V)
+tEnd = 100
+
+#numerical parameters
+lcmin = 0.0001 # approx r*100 but should match the mesh size
+dt = 5e-4
+alpha = 5e-3
+epsilon = alpha*lcmin**2 /nu
+print('epsilon',epsilon)
+
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
+p.write(outputdir+"/part-00000")
+rout = 1e-3
+mu = nu*rho
+
+genInitialPosition(outputdir + "/part-00000", r, rout, rhop)
+
+p = scontact2.ParticleProblem(outputdir+"/part-00000")
+
+print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
+print("RHOP = %g" % rhop)
+outf = 50
+outf1 = 100000
+
+strong_boundaries = [("Top",2,0.),("Top",1,0.),("Top",0,0.),("Box",1,0.),("Box",0,0.)]
+fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
+
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+fluid.export(outputdir,0,0)
+
+ii = 0
+jj = 0
+porosity = len(p.velocity()[:,1])*r**2/(rout**2)
+p.velocity()[:,1] = 2./9. *porosity* rout**2*g*(rhop-rho)/mu
+print(p.velocity()[:,1])
+while ii<1000:
+    if (ii%200==0 and ii != 0):
+       fluid.export(outputdir1, t, jj)
+       fluid.adapt_mesh(0.01,50,1e-4,7500)
+       forces = fluid.compute_node_force(dt)
+       fluid.export(outputdir1, t, jj+1)
+    forces = fluid.compute_node_force(dt)
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    fluid.implicit_euler(dt)
+    if (ii%200==0 and ii != 0):
+       fluid.export(outputdir1, t, jj+2)
+       jj = jj+3
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        p.write_vtk(outputdir, ioutput, t)
+        p.write(outputdir+"/part-%05d" % ioutput)
+        p.write_boundary_vtk(outputdir, ioutput, t)
+        fluid.export(outputdir, t, ioutput)
+    ii += 1
+
+
+tic = time.clock()
+forces = g*p.mass()
+while t < tEnd : 
+    if (ii%200==0 and ii != 0):
+       fluid.export(outputdir1, t, jj)
+       fluid.adapt_mesh(0.01,50,1e-4,7500)
+       forces = fluid.compute_node_force(dt)
+       fluid.export(outputdir1, t, jj+1)
+    forces = fluid.compute_node_force(dt)
+    vn = p.velocity() + forces * dt / p.mass()
+    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
+    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()))
+    for i in range(nsub) :
+        p.iterate(dt/nsub, forces)  
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    fluid.implicit_euler(dt)
+    if (ii%200==0 and ii != 0):
+       fluid.export(outputdir1, t, jj+2)
+       jj = jj+3
+    t += dt
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        p.write_vtk(outputdir, ioutput, t)
+        p.write(outputdir+"/part-%05d" % ioutput)
+        p.write_boundary_vtk(outputdir, ioutput, t)
+        fluid.export(outputdir, t, ioutput)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))

testcases/drop-2d/dropNewAdapt/dropE5.py

+#!/usr/bin/env python
+from marblesbag import fluid as fluid
+from marblesbag import scontact2
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+
+
+def genInitialPosition2(filename, r, rout, rhop) :
+    D = 3.8e-4/2
+    p = scontact2.ParticleProblem()
+    p.load_msh_boundaries("mesh.msh", ["Top", "Box"])
+    x = np.arange(rout, -rout, -D)
+    x, y = np.meshgrid(x, x)
+    R2 = x**2 + y**2
+    keep = R2 < (rout - r)**2
+    x = x[keep]
+    y = y[keep]
+    for i in range(x.shape[0]) :
+        if y[i]<0 :
+            p.add_particle((x[i]-D/2.1+2*random.uniform(0,D/4), y[i]-D/2.1+2*random.uniform(0,D/4)), r, r**2 * np.pi * rhop);
+    
+    x = np.arange(rout, -rout, -D)
+    y = np.arange(2*r, 4*rout, D)
+    x, y = np.meshgrid(x, y)
+    x = x.flat
+    y = y.flat
+    T = rout
+    x1 = -4/5*rout
+    x2 = 4/5*rout
+    for i in range(len(x)) :
+        if y[i] < -T/x1*x[i]+T and y[i] < -T/x2*x[i]+T :
+            p.add_particle((x[i]-D/2.1+2*random.uniform(0,D/4.1), y[i]-D/2.1+2*random.uniform(0,D/4.1)), r, r**2 * np.pi * rhop);
+    p.position()[:, 1] += 0.45
+    p.write(filename)
+
+outputdir = "outputE5"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+t = 0
+ii = 0
+
+r=25e-6
+p = scontact2.ParticleProblem()
+#R = np.random.uniform(45e-06, 90e-06, len(x))
+
+#physical parameters
+g =  -9.81
+rho = 2103
+rhop = 2400
+nu = 1e-4
+V = 0.5 # todo : estimate V base on limit velocity
+print('V',V)
+tEnd = 100
+
+#numerical parameters
+lcmin = 0.0001 # approx r*100 but should match the mesh size
+dt = 5e-4
+alpha = 1e-2
+epsilon = alpha*lcmin**2 /nu#2e-2*c*h*2#2e-2*c*h # ?? not sure ??1e-5
+print('epsilon',epsilon)
+
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
+p.write(outputdir+"/part-00000")
+rout = 3.5e-3
+genInitialPosition2(outputdir + "/part-00000", r, rout, rhop)
+
+p = scontact2.ParticleProblem(outputdir+"/part-00000")
+
+print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
+print("RHOP = %g" % rhop)
+outf = 50
+outf1 = 100000
+
+strong_boundaries = [("Top",2,0.),("Top",1,0.),("Top",0,0.),("Box",1,0.),("Box",0,0.)]
+fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
+
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+fluid.export(outputdir,0,0)
+
+mu = nu*rho
+porosity = len(p.velocity()[:,1])*r**2/(rout**2)
+p.velocity()[:,1] = 2./9. *porosity* rout**2*g*(rhop-rho)/mu
+while ii<300:
+    if (ii%100==0 and ii != 0) or ii==50:
+       fluid.adapt_mesh(0.01,50,0.01,50,1e-4,8000)
+    forces = fluid.compute_node_force(dt)
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    fluid.implicit_euler(dt)
+    ii += 1
+    
+ii = 0
+tic = time.clock()
+while t < tEnd : 
+    if (ii%100==0 and ii != 0) or ii==50:
+       fluid.adapt_mesh(0.01,50,0.01,50,3e-4,8000)
+    forces = fluid.compute_node_force(dt)
+    vn = p.velocity() + forces * dt / p.mass()
+    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
+    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()))
+    for i in range(nsub) :
+        p.iterate(dt/nsub, forces)  
+    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    fluid.implicit_euler(dt)
+    t += dt
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        p.write_vtk(outputdir, ioutput, t)
+        p.write(outputdir+"/part-%05d" % ioutput)
+        p.write_boundary_vtk(outputdir, ioutput, t)
+        fluid.export(outputdir, t, ioutput)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))

testcases/drop-2d/dropNewAdapt/mesh.geo

+L = 0.025;
+H = 0.5;
+y = 0;
+lc = 0.001;
+
+Point(1) = {-L, H, 0};
+Point(2) = {-L, 0, 0};
+Point(3) = {L, 0, 0};
+Point(4) = {L, H, 0};
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Point(5) = {0,0.4,0};
+Point(6) = {0,H,0};
+Line(5) = {5,6};
+
+Line Loop(1) = {1:4};
+
+Plane Surface(1) = {1};
+Physical Line("Box") = {1,2,3};
+Physical Line("Top") = {4};
+Physical Surface("Domain") = {1};
+Physical Point("PtFix") = {1};
+
+Field[1] = Attractor;
+Field[1].EdgesList = {5};
+Field[1].NNodesByEdge = 200;
+
+Field[2] = Threshold;
+Field[2].DistMax = 0.08;
+Field[2].DistMin = 0.008;
+Field[2].LcMax = 0.05;
+Field[2].LcMin = 0.001;
+Field[2].IField = 1;
+
+Background Field = 2;
+
+Mesh.CharacteristicLengthFromPoints = 0;
+Mesh.CharacteristicLengthExtendFromBoundary = 0;