Commit 4caffdf3 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

forgot validation again

parent 06ddaea0
Pipeline #8251 passed with stages
in 3 minutes and 3 seconds
......@@ -32,12 +32,12 @@ import time
import shutil
import random
import sys
listeVit = [14.39,36.32,64.06,95.57,130.75,169.32,211.05,255.63,303.22,354.05,408.80,467.48,530.59]
listeVit = [14,36,64,95,130,169,211,255,303,354,408,467]
inter = int(sys.argv[1])
index = int(inter/5)
numDep = inter%5
# Physical parameters
vit = -0.001*listeVit[index] #stream velocity
vit = -0.001*listeVit[index] #stream velocity
muwall = 0.3 #friction coefficient between the walls and the particles
mupart = 0.15 #friction coefficient between particles
rho = 1000 #fluid density
......@@ -46,9 +46,9 @@ g = np.array([0,-9.81,0]) #gravity
# Numerical parameters
dt = 1e-3 #time step
t = 0 #initial time
tEnd = 10 #final time
tEnd = 30 #final time
ii = 0 #iteration number
outf = 100 #iterations between data frames
outf = 2500 #iterations between data frames
# Define output directory
outputdir = "drag" + str(int(listeVit[index])) + "/" + str(numDep)
if not os.path.isdir(outputdir) :
......@@ -64,12 +64,12 @@ if os.path.exists(filename):
#
#Injected particles
p = scontact.ParticleProblem(3,True,True)
p.read_vtk("outputhahalol",5)
p.read_vtk("depot/"+ str(numDep),2)
p.remove_particles_flag((p.position()[:,1]-p.r()[:,0]>0.5)*(p.position()[:,1]+p.r()[:,0] <0.52))
p.velocity()[:,1] = vit
#Flow particles
p2 = scontact.ParticleProblem(3,True,True)
p2.read_vtk("outputhahalol",5)
p.read_vtk("depot/"+ str(numDep),2)
p2.remove_particles_flag(p2.position()[:,1]+p2.r()[:,0] <0.5)
p2.clear_boundaries()
p2.load_msh_boundaries("mesh.msh", ["Inner","Left", "Right","Front","Rear"],material = "Steel")
......@@ -79,17 +79,17 @@ p2.write_vtk(outputdir, 0, 0)
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(3,g,[nu*rho],[rho],drag_in_stab=0)
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Inner")
fluid.set_wall_boundary("Left", velocity = [0,0,0])
fluid.set_wall_boundary("Right", velocity = [0,0,0])
fluid.set_open_boundary("Bottom", velocity = [0,vit,0])
fluid.set_wall_boundary("Front", velocity = [0,0,0])
fluid.set_wall_boundary("Rear", velocity = [0,0,0])
fluid.set_open_boundary("Top",pressure = 0)
fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces(),init=True)
fluid.export_vtk(outputdir, 0.,0)
#fluid = fluid.FluidProblem(3,g,[nu*rho],[rho],drag_in_stab=0)
#fluid.load_msh("mesh.msh")
#fluid.set_wall_boundary("Inner")
#fluid.set_wall_boundary("Left", velocity = [0,0,0])
#fluid.set_wall_boundary("Right", velocity = [0,0,0])
#fluid.set_open_boundary("Bottom", velocity = [0,vit,0])
#fluid.set_wall_boundary("Front", velocity = [0,0,0])
#fluid.set_wall_boundary("Rear", velocity = [0,0,0])
#fluid.set_open_boundary("Top",pressure = 0)
#fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces(),init=True)
#fluid.export_vtk(outputdir, 0.,0)
#
# COMPUTATION LOOP
#
......@@ -103,9 +103,10 @@ while t < tEnd :
p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
p2.remove_particles_flag( (p2.position()[:,1] + p2.r()[:,0] >-0.52))
F = np.zeros(3)
time_integration.particle_changed(fluid,p2)
time_integration.iterate(fluid,p2,dt,1,contact_tol=1e-6,external_particles_forces=g*p2.mass(),after_sub_iter=lambda subdt : accumulate(F))
#time_integration.particle_changed(fluid,p2)
time_integration.iterate(None,p2,dt,1,contact_tol=1e-6,external_particles_forces=g*p2.mass(),after_sub_iter=lambda subdt : accumulate(F))
with open(filename,"a") as file1:
F /= dt
F += fluid.boundary_force()[0,:]
meanVel = abs((np.mean(np.linalg.norm(p2.velocity()[(p2.position()[:,1] > -0.45)*(p2.position()[:,1]<-0.3),:2],axis=1)) + np.mean(np.linalg.norm(p2.velocity()[(p2.position()[:,1] < 0.45)*(p2.position()[:,1]>0.3),:2],axis=1)))/2)
file1.write( str(F[0])+";"+str(F[1])+";"+str(F[2])+";"+str(t)+";"+str(meanVel)+"\n")
......@@ -114,6 +115,6 @@ while t < tEnd :
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p2.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
#fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g" % (ii, t, tEnd))
......@@ -124,7 +124,7 @@ class Weight(unittest.TestCase) :
ioutput = int(ii/outf) + 1
#p.write_vtk(outputdir, ioutput, t)
print("%i : %.2g/%.2g" % (ii, t, tEnd))
error = np.sqrt(sum(sum(forces) - bnd_forces)**2)
error = np.sqrt(sum(sum(forces) - bnd_forces/dt)**2)
error /= np.sqrt(sum(sum(forces)**2))
print(error*100,"%")
self.assertLess(error,.5/100., "error is too large in weight")
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment