Commit 9574aeda authored by Matthieu Constant's avatar Matthieu Constant
Browse files

correction of minor errors + source term in velocity

parent 056803bb
......@@ -241,7 +241,7 @@ static void f_u_0(
for (int k = 0; k < D; ++k) {
utau += u[k]*(du[j][k]-u[j]*dq[k]/q);
}
f[j] = rho*((u[j]-uold[j])/dt + u[j]/q*divu + utau/q);
f[j] = rho*((u[j]-uold[j])/dt + u[j]/q*divu + utau/q) + (j==1 ? .1 : 0);
}
}
......@@ -304,7 +304,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
all_local_vector[i] = 0;
for (size_t i = 0; i < N_E*local_size*local_size; ++i)
all_local_matrix[i] = 0;
compute_node_force_implicit(problem, dt, all_local_vector, all_local_matrix);
//compute_node_force_implicit(problem, dt, all_local_vector, all_local_matrix);
for (int iel=0; iel < mesh->n_elements; ++iel) {
const unsigned int *el = &mesh->elements[iel*N_N];
double dxdxi[D][D], dxidx[D][D], dphi[N_N][D];
......@@ -376,6 +376,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
}
buf = dq[i];
dq[i] += deps;
f_u_0(fu0e,u,du,q,dq,uold,dt,rho);
f_u_1(fu1e,u,du,q,dq,p,mu);
dq[i] = buf;
......@@ -431,13 +432,14 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
local_vector[iphi+N_SF*j] += jw*f;
}
double f = fp0*phii;
double ff = fp0*phii;
for (int j = 0; j < D; ++j){
f += fp1[j]*dphii[j];
ff += fp1[j]*dphii[j];
}
local_vector[iphi+N_SF*P] += jw*f;
local_vector[iphi+N_SF*P] += jw*ff;
local_vector[iphi+N_SF*Q] += jw*fq0*phii;
//printf("iphi=%d,iel=%d: local_vectorU=%g,local_vectorV=%g,local_vectorP=%g,local_vectorQ=%g\n",iphi,iel,local_vector[iphi+N_SF*0],local_vector[iphi+N_SF*1],local_vector[iphi+N_SF*P],local_vector[iphi+N_SF*Q]);
for (int jphi = 0; jphi < N_SF; ++jphi){
double phij = phi[jphi];
......@@ -465,7 +467,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
mu_u += fu0_du[j][k][i]*dphij[i]*phii+fu1_u[j][k][i]*phij*dphii[i];
mu_q += fu1_dq[j][k][i]*dphij[k]*dphii[i];
for (int l = 0; l < D; ++l){
mu_u+= fu1_du[j][k][i][l]*dphij[i]*dphii[l];
mu_u += fu1_du[j][k][i][l]*dphij[i]*dphii[l];
}
}
LOCAL_MATRIX(U,V) = jw*mu_u;
......
......@@ -58,7 +58,7 @@ def genInitialPosition(filename, r, rout, rhop, compacity) :
p.write_vtk(filename,0,0)
outputdir = "output"
outputdir = "outputTest"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
......@@ -169,5 +169,5 @@ while t < tEnd :
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
if ii==2:
if ii==1:
break
......@@ -28,8 +28,8 @@ Field[1].NodesList = {5};
Field[2] = Threshold;
Field[2].DistMax = 0.08;
Field[2].DistMin = 0.02;
Field[2].LcMax = 0.02;
Field[2].LcMin = 0.001;
Field[2].LcMax = 0.8;
Field[2].LcMin = 0.04;
Field[2].IField = 1;
Background Field = 2;
......
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