Commit f01f2717 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

implicit force velocity in force OK

parent d0bc9bac
......@@ -70,9 +70,9 @@ class ParticleProblem :
segmentTag = self._get_idx("SegmentTag")
segments[segmentTag == tag, 4:6] = v
def iterate(self, dt, forces) :
def iterate(self, dt, forces,tol=1e-8) :
self._get_matrix("ExternalForces",2)[:] = forces
lib.particleProblemIterate(self._p, c_double(numpy.min(self.r())), c_double(dt), c_double(1e-8), c_int(-1))
lib.particleProblemIterate(self._p, c_double(numpy.min(self.r())), c_double(dt), c_double(tol), c_int(-1))
def write_vtk(self, odir, i, t) :
lib.particleProblemWriteVtk(self._p, odir.encode(), i)
......
......@@ -134,7 +134,6 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
static void compute_node_force_implicit(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix) {
double *porosity = problem->porosity;
Mesh *mesh = problem->mesh;
const double *solution = problem->solution;
for (int i = 0; i < problem->n_particles; ++i) {
int iel = problem->particle_element_id[i];
if(iel < 0){
......@@ -145,7 +144,7 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
shape_functions(xi,phi);
uint32_t *el=problem->mesh->elements+iel*N_SF;
double *se = problem->solution_explicit;
double c=0, vf[DIMENSION]={0}, vfe[DIMENSION]={0}, p=0, gradp[DIMENSION]={0};
double c=0, vf[DIMENSION]={0}, vfe[DIMENSION]={0}, gradp[DIMENSION]={0};
double dxdxi[DIMENSION][DIMENSION],dxidx[DIMENSION][DIMENSION];
for (int k=0; k<DIMENSION; ++k)
for (int j=0; j<DIMENSION; ++j)
......@@ -153,14 +152,14 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
invDD(dxdxi,dxidx);
double dphi[N_SF][DIMENSION];
grad_shape_functions(dxidx, dphi);
double *si = problem->solution;
for (int iphi=0; iphi < N_SF; ++iphi) {
c += porosity[el[iphi]]*phi[iphi];
for (int j=0; j < DIMENSION; ++j) {
vf[j] += solution[el[iphi]*n_fields+j]*phi[iphi];
vf[j] += si[el[iphi]*n_fields+j]*phi[iphi];
vfe[j] += se[el[iphi]*n_fields+j]*phi[iphi];
gradp[j] += dphi[iphi][j]*solution[el[iphi]*n_fields+DIMENSION];
gradp[j] += dphi[iphi][j]*si[el[iphi]*n_fields+DIMENSION];
}
p += solution[el[iphi]*n_fields+DIMENSION]*phi[iphi];
}
double du[DIMENSION], due[DIMENSION];
for (int j = 0; j < DIMENSION; ++j) {
......@@ -183,8 +182,9 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
double *local_matrix = all_local_matrix+iel*n_fields*n_fields*N_SF*N_SF;
for (int iphi = 0; iphi < N_SF; ++iphi) {
for (int d = 0; d < DIMENSION; ++d)
local_vector[iphi+N_SF*d] += fcoeff*gamma*(du[d]/*-dt/mass*vol*gradp[d]*/)*phi[iphi];
local_vector[iphi+N_SF*1] += fcoeff*gamma*(-dt*g)*phi[iphi];
local_vector[iphi+N_SF*d] -= fcoeff*gamma*(du[d]-dt/mass*vol*gradp[d])*phi[iphi];
local_vector[iphi+N_SF*1] -= fcoeff*gamma*dt*g*phi[iphi];
}
for (int iphi = 0; iphi < N_SF; ++iphi) {
for (int jphi = 0; jphi < N_SF; ++jphi) {
......@@ -193,8 +193,8 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
double f = fcoeff*phi[iphi];
for (int d = 0; d < DIMENSION; ++d){
LOCAL_MATRIX(d,d) += -f/c*phi[jphi]*gamma;
LOCAL_MATRIX(d,DIMENSION) += -f*gamma*dt/mass*vol*dphi[jphi][d];
LOCAL_MATRIX(d,d) -= -f/c*phi[jphi]*gamma;
LOCAL_MATRIX(d,DIMENSION) -= -f*gamma*dt/mass*vol*dphi[jphi][d];
}
}
}
......@@ -242,12 +242,11 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int iqp = 0; iqp< N_QP; ++iqp) {
double phi[N_SF];
shape_functions(QP[iqp],phi);
double u[DIMENSION]={0}, dudt[DIMENSION]={0}, fu[DIMENSION]={0}, p=0, c=0, dcdt=0;
double u[DIMENSION]={0}, dudt[DIMENSION]={0}, p=0, c=0, dcdt=0;
for (int i = 0; i < N_SF; ++i) {
for (int j = 0; j < DIMENSION; ++j){
u[j] += solution[el[i]*n_fields+j]*phi[i];
dudt[j] += (solution[el[i]*n_fields+j]-solution_old[el[i]*n_fields+j])/dt*phi[i];
//fu[j] += node_force[el[i]*DIMENSION+j]*phi[i];
}
p += solution[el[i]*n_fields+DIMENSION]*phi[i];
c += porosity[el[i]]*phi[i];
......@@ -265,13 +264,15 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
double utau = 0;
double dphisigma=0;
for (int k = 0; k < DIMENSION; ++k) {
utau += u[k]*dphii[k]*u[j]/c;
//utau += u[k]*dphii[k]*u[j]/c;
utau += u[k]*tau[j][k];
dphisigma += 2*mu*dphii[k]*0.5*(tau[k][j]+tau[j][k]);
}
local_vector[iphi+N_SF*j] += jw*(
dphisigma
+rho*phii*dudt[j]-rho*utau
-p*dphii[j]-fu[j]*phii
//+rho*phii*dudt[j]-rho*utau
+rho*phii*(dudt[j]-u[j]/c*dcdt+utau/c)
-p*dphii[j]
);
}
double dphidp = 0, divu = 0;
......@@ -296,7 +297,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int j = 0; j < DIMENSION; ++j) {
dtau[j] = dphij[j]-phij*dc[j]/c;
dphiidphij += dphii[j]*dphij[j];
udtau += u[j]*phij*dphii[j]/c;
udtau += u[j]*dtau[j];
dphiidtau += dphii[j]*dtau[j];
}
for (int j = 0; j < DIMENSION; ++j) {
......@@ -305,9 +306,11 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
LOCAL_MATRIX(P,U) += jw*phii*dphij[j];
for (int k = 0; k < DIMENSION; ++k) {
int V = k;
LOCAL_MATRIX(U,V) += -jw*rho*u[j]*dphii[k]*phij/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
//LOCAL_MATRIX(U,V) += -jw*rho*u[j]*dphii[k]*phij/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
LOCAL_MATRIX(U,V) += jw*rho*phii*phij*tau[j][k]/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
}
LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*(phii*phij/dt-phij*dphii[j]*u[j]/c);
//LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*(phii*phij/dt-phij*dphii[j]*u[j]/c);
LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*phii*(phij/dt-phij/c*dcdt+udtau/c);
}
LOCAL_MATRIX(P,P) += jw*epsilon[el[iphi]]*dphiidphij;
}
......
L = 0.025;
H = 0.025;
L = 0.5;
H = 0.5;
y = 0;
lc = 0.0015;
lc = 0.03;
Point(1) = {-L, H, 0, lc};
Point(2) = {-L, -H, 0, lc};
......
......@@ -11,8 +11,8 @@ import random
def genInitialPosition(filename, r, ly, rhop) :
p = scontact2.ParticleProblem()
p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom","BottomOut","TopOut"])
x = np.arange(-0.025+r, 0.025-r, 2*r)
y = np.arange(-0.025+r, -0.025+r+ly, 2*r)
x = np.arange(-0.5+r, 0.5-r, 2*r)
y = np.arange(-ly/2+r, ly/2-r, 2*r)
x, y = np.meshgrid(x, y)
x = x.flat
y = y.flat
......@@ -29,8 +29,8 @@ if not os.path.isdir(outputdir) :
t = 0
ii = 0
r = 2e-4
ly = 5e-2
r = 5e-3
ly = 0.7
p = scontact2.ParticleProblem()
#R = np.random.uniform(45e-06, 90e-06, len(x))
......@@ -39,14 +39,12 @@ g = -9.81
rho = 1000
rhop = 1500
nu = 1e-6
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 100
#numerical parameters
lcmin = 0.001 # approx r*100 but should match the mesh size
dt = 1e-4*10
alpha = 2.5e-6
lcmin = 0.01 # approx r*100 but should match the mesh size
dt = 1e-4*100
alpha = 2.5e-6/50
epsilon = alpha*lcmin**2 /nu
print('epsilon',epsilon)
......@@ -61,21 +59,15 @@ p.read_vtk(outputdir,0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
outf = 5
outf = 10
outf1 = 100000
ii = 0
def outerBndV(x) :
bndLoop = 12000
if ii%bndLoop<11000 and ii%bndLoop>9500:
print(min(((ii-9500)%bndLoop)*5e-5,.05))
return min(((ii-9500)%bndLoop)*5e-5,.05)
else:
print(.05/((ii-11000)%bndLoop + 1))
return .05/((ii-11000)%bndLoop + 1)
strong_boundaries = [("Top",2,0.),("TopOut",1,outerBndV),("TopOut",0,0.),("Top",1,outerBndV),("BottomOut",1,outerBndV),("BottomOut",0,0.),("Bottom",1,outerBndV),("Lateral",0,0.)]
return 0.1 * max(np.sin(-t*np.pi*2./5-4*np.pi/5),0)
strong_boundaries = [("Top",2,0.),("TopOut",1,outerBndV),("TopOut",0,0.),("Top",1,outerBndV),("BottomOut",1,outerBndV),("BottomOut",0,0.),("Bottom",1,outerBndV),("Bottom",0,0),("Lateral",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())
......@@ -92,7 +84,7 @@ while t < tEnd :
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)
p.iterate(dt/nsub, forces,1e-6)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
t += dt
if ii %outf == 0 :
......
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