Commit 73a154c7 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

vector ok

parent 45a78fd8
......@@ -21,7 +21,7 @@
ALL: libmbfluid3.so libmbfluid2.so libscontact2.so libscontact3.so
CFLAGS=-Wno-unused-function -O3 -g -march=native -mtune=native -fPIC -std=gnu99 -Iscontact -Ihxt -Isrc
CFLAGS=-Wno-unused-function -g -march=native -mtune=native -fPIC -std=gnu99 -Iscontact -Ihxt -Isrc
LDFLAGS=-shared -lm
FLUID_C=src/fluid_problem.c src/mesh.c src/mesh_find.c scontact/quadtree.c hxt/hxt_linear_system.c hxt/hxt_linear_system_lu.c hxt/hxt_message.c src/mbxml.c src/yxml.c src/fluid_problem_io.c
......
......@@ -313,10 +313,9 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
const double detj = invDD(dxdxi, dxidx);
grad_shape_functions(dxidx, dphi);
double *local_vector = &all_local_vector[N_SF*n_fields*iel];
double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*iel];
double dq[D]={0}, du[D][D]={{0}}, dp[D]={0};
int P = (D+1)*n_fluids+1;
double *local_vector = &all_local_vector[local_size*iel];
double *local_matrix = &all_local_matrix[local_size*local_size*iel];
int P = (D+1)*n_fluids;
for (int iqp = 0; iqp< N_QP; ++iqp) {
double phi[N_SF];
shape_functions(QP[iqp],phi);
......@@ -326,10 +325,12 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
const double jw = QW[iqp]*detj;
for (int iphi = 0; iphi < N_SF; ++iphi){
#define LOCAL_VECTOR(a) local_vector[iphi+N_SF*(a)]
double phii = phi[iphi];
local_vector[N_SF*P+iphi] += (1-c)*phii*jw;
LOCAL_VECTOR(P) += (1-c)*phii*jw;
}
for (int nf = 0; nf < D; ++nf){
for (int nf = 0; nf < problem->n_fluids; ++nf){
double dq[D]={0}, du[D][D]={{0}}, dp[D]={0};
double rho = problem->rho[nf];
double mu = problem->mu[nf];
int Q = (D+1)*nf+D;
......@@ -436,7 +437,6 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
//filling
for (int iphi = 0; iphi < N_SF; ++iphi){
#define LOCAL_VECTOR(a) local_vector[iphi+N_SF*(a)]
double phii = phi[iphi];
const double *dphii = dphi[iphi];
for (int k = 0; k < D; ++k){
......@@ -509,9 +509,9 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
}
}
hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
}
hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
}
free(all_local_matrix);
free(all_local_vector);
......@@ -674,6 +674,8 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
problem->n_fluids = n_fluids;
int n_fields = problem_n_fields(problem);
problem->mu = (double*)malloc(sizeof(double)*n_fluids);
problem->rho = (double*)malloc(sizeof(double)*n_fluids);
for (int i = 0; i < n_fluids; ++i){
problem->rho[i] = rho[i];
problem->mu[i] = mu[i];
......@@ -702,7 +704,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
problem->solution[i] = 0.;
}
for (int i = 0; i < mesh->n_nodes; ++i){
problem->solution[i*n_fields+(D+1)] = 1.;
problem->solution[i*n_fields+D] = 1.;
}
problem->node_volume = NULL;
fluid_problem_compute_node_volume(problem);
......@@ -726,6 +728,8 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
void fluid_problem_free(FluidProblem *problem)
{
free(problem->mu);
free(problem->rho);
free(problem->solution);
free(problem->compacity);
free(problem->old_compacity);
......@@ -1061,11 +1065,6 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->old_compacity[i] = problem->compacity[i];
}
if (init) {
int Q = D+1;
for (int i = 0; i < problem->mesh->n_nodes; ++i)
problem->solution[i*n_fields+Q] = problem->compacity[i];
}
clock_get(&tic);
fluid_problem_compute_compacity(problem);
clock_get(&toc);
......
......@@ -59,8 +59,8 @@ outf = 10 # number of iterations between o
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("Left",1,0),("Right",1,0),("Top",0,0),("Bottom",0,0),("Right",2,0),("Left",0,lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]))]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries,1)
strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("Left",1,0),("Right",1,0),("Top",0,0),("Bottom",0,0),("Right",3,0),("Left",0,lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]))]
fluid = fluid.fluid_problem("mesh.msh",g,[nu*rho],[rho],epsilon,strong_boundaries,1,1)
ii = 0
t = 0
......@@ -79,4 +79,4 @@ while t < tEnd :
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
\ No newline at end of file
Markdown is supported
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