Commit dbd81b51 by Matthieu Constant

### Merge branch 'a_eq' of git.immc.ucl.ac.be:fluidparticles/migflow into a_eq

parents d2f2242c 110c3aa3
Pipeline #4887 passed with stage
in 26 seconds
 ... ... @@ -33,6 +33,8 @@ #define P D #define U 0 #define LOCAL_MATRIX(i,j,a,b) local_matrix[((i)*n_fields*N_SF+N_SF*(a)+(j))*n_fields+(b)] int fluid_problem_n_fields(const FluidProblem *p) {return D+1;} #if D==2 ... ... @@ -207,8 +209,6 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old, particle_force_f(problem,h0+jd*n_fields,s,ds,sold,c,dc,a,dt,iel,ip,1,reduced_gravity); ds[jfield*D+jd] = store; } int N2 = n_fields*N_SF; #define LOCAL_MATRIX2(a,b) local_matrix[IDX+N2*(a)+(b)] for (int ifield = 0; ifield < n_fields; ++ifield){ for (int iphi = 0; iphi < N_SF; ++iphi){ for (int jphi = 0; jphi < N_SF; ++jphi){ ... ... @@ -218,8 +218,7 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old, double df01 = (h0[jd*n_fields+ifield]-f0[ifield])/deps; m += phi[iphi]*dphi[jphi][jd]*df01; } int IDX = (iphi*N2)*n_fields + jphi*n_fields; LOCAL_MATRIX2(ifield,jfield) += m; LOCAL_MATRIX(iphi,jphi,ifield,jfield) += m; } } } ... ... @@ -267,7 +266,6 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n for (int id = 0; id < D; ++id) { f[U+id] = (pid<0?p:data[pid])*n[id] + (un> 0 ? rho*u[id]*un/c : 0) + (vid<0?0:(1+D)/(D*h)*mu*(1+problem->tauc[eid])*N_N*((u[id]-data[vid+id])+s_c*n[id])) - mu/c*(dun[id]-u[id]/c*dcn+dunt[id]-un/c*dc[id]); } //printf("%s:velocity=[%g,%g]\n",wbnd->tag,data[vid],data[vid+1]); } else { double p = pid<0 ? s[P] : data[pid]; ... ... @@ -421,7 +419,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d double *ds = malloc(sizeof(double)*n_fields*D); double *f0 = malloc(sizeof(double)*n_fields); double *g0 = malloc(sizeof(double)*n_fields); double *h0 = malloc(sizeof(double)*n_fields*D); double *h0 = malloc(sizeof(double)*n_fields*N_LSF); for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){ WeakBoundary *wbnd = problem->weak_boundaries + ibnd; int bndid= -1; ... ... @@ -442,11 +440,12 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4]; const int iel = interface[0]; const int icl = interface[1]; int nodes[D]; double *x[D]; const int *cl = elbnd[icl]; int nodes[N_LSF]; double *x[N_LSF]; const int *el = &mesh->elements[iel*N_N]; for (int i = 0; i < D; ++i){ nodes[i] = el[elbnd[icl][i]]; nodes[i] = el[cl[i]]; x[i] = &mesh->x[nodes[i]*3]; } double dxdxi[D][D], dxidx[D][D], dphi[N_N][D]; ... ... @@ -457,59 +456,41 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d double n[D],detbnd; get_normal_and_det(x,n,&detbnd); grad_shape_functions(dxidx, dphi); double dc[D] = {0}; for (int j = 0; j < n_fields*D; ++j) { ds[j] = 0; } for (int i = 0; i < N_SF; ++i) { for (int j = 0; j < n_fields; ++j) { double dof = solution[el[i]*n_fields+j]; for (int k = 0; k < D; ++k) ds[j*D+k] += dphi[i][k]*dof; } for (int k=0; kporosity[el[i]]*dphi[i][k]; } } for (int iqp = 0; iqp < N_LQP; ++iqp) { double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)]; double xi[D]; { #if DIMENSION == 2 double a = LQP[iqp][0]; switch(icl) { case 0 : xi[0] = a; xi[1] = 0; break; case 1 : xi[0] = 1-a; xi[1] = a; break; case 2 : xi[0] = 0; xi[1] = 1-a; break; default:; } #else double a = LQP[iqp][0]; double b = LQP[iqp][1]; switch(icl) { case 0 : xi[0] = a; xi[1] = b; xi[2] = 0; break; case 1 : xi[0] = 0; xi[1] = a; xi[2] = b; break; case 2 : xi[0] = b; xi[1] = 0; xi[2] = a; break; case 3 : xi[0] = 1-a-b; xi[1] = b; xi[2] = a; break; default:; } #endif } double dp[D]={0}; double *local_vector = &all_local_vector[local_size*iel]; double *local_matrix = &all_local_matrix[local_size*local_size*iel]; double phi[N_SF]; shape_functions(xi,phi); double phi[DIMENSION]; l_shape_functions(LQP[iqp],phi); for (int i = 0; i < n_fields; ++i) { s[i] = 0; sold[i] = 0; for (int j = 0; j < D; ++j) { ds[i*D+j] = 0; } } double c = 0; double dc[D] = {0}; double a = 0; for (int i = 0; i < N_SF; ++i) { c += problem->porosity[el[i]]*phi[i]; for (int i = 0; i < DIMENSION; ++i) { c += problem->porosity[nodes[i]]*phi[i]; if (problem->n_fluids == 2) { a += problem->concentration[iel*N_N+i]*phi[i]; a += problem->concentration[iel*N_N+cl[i]]*phi[i]; } for (int j = 0; j < n_fields; ++j) { double dof = solution[el[i]*n_fields+j]; double dof = solution[nodes[i]*n_fields+j]; s[j] += phi[i]*dof; for (int k = 0; k < D; ++k) ds[j*D+k] += dphi[i][k]*dof; } for (int k=0; kporosity[el[i]]*dphi[i][k]; } } double rho = problem->rho[0], mu = problem->mu[0]; ... ... @@ -520,8 +501,8 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d const double jw = LQW[iqp]*detbnd; f_boundary(wbnd,problem,n,f0,s,ds,c,dc,rho,mu,dt,iel,dataqp); for (int ifield = 0; ifield < n_fields; ++ifield) { for (int iphi = 0; iphi < N_SF; ++iphi){ local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield]*jw; for (int iphi = 0; iphi < D; ++iphi){ local_vector[cl[iphi]+N_SF*ifield] += phi[iphi]*f0[ifield]*jw; } } for (int jfield = 0; jfield < n_fields; ++jfield) { ... ... @@ -535,19 +516,19 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d f_boundary(wbnd,problem,n,h0+jd*n_fields,s,ds,c,dc,rho,mu,dt,iel,dataqp); ds[jfield*D+jd] = store; } int N2 = n_fields*N_SF; #define LOCAL_MATRIX2(a,b) local_matrix[IDX+N2*(a)+(b)] for (int ifield = 0; ifield < n_fields; ++ifield){ for (int iphi = 0; iphi < N_SF; ++iphi){ for (int jphi = 0; jphi < N_SF; ++jphi){ for (int iphi = 0; iphi < N_LSF; ++iphi){ for (int jphi = 0; jphi < N_LSF; ++jphi){ double df00 = (g0[ifield]-f0[ifield])/deps; double m = jw*phi[iphi]*phi[jphi]*df00; for (int jd = 0; jd < D; ++jd) { LOCAL_MATRIX(cl[iphi],cl[jphi],ifield,jfield) += jw*phi[iphi]*phi[jphi]*df00; } for (int jphi = 0; jphi < N_SF; ++jphi){ double m = 0; for (int jd = 0; jd < N_LSF; ++jd) { double df01 = (h0[jd*n_fields+ifield]-f0[ifield])/deps; m += jw*phi[iphi]*dphi[jphi][jd]*df01; } int IDX = (iphi*N2)*n_fields + jphi*n_fields; LOCAL_MATRIX2(ifield,jfield) += m; LOCAL_MATRIX(cl[iphi],jphi,ifield,jfield) += m; } } } ... ... @@ -592,19 +573,36 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o grad_shape_functions(dxidx, dphi); double *local_vector = &all_local_vector[local_size*iel]; double *local_matrix = &all_local_matrix[local_size*local_size*iel]; for (int i = 0; i < n_fields; ++i) { for (int j = 0; j < D; ++j) { ds[i*D+j] = 0; } } double dc[D] = {0}, da[D]={0}; for (int i = 0; i < N_SF; ++i) { for (int j = 0; j < n_fields; ++j) { double dof = solution[el[i]*n_fields+j]; for (int k = 0; k < D; ++k) ds[j*D+k] += dphi[i][k]*dof; } for (int k=0; kporosity[el[i]]*dphi[i][k]; } if (problem->n_fluids==2) { for (int k=0; kconcentration[iel*N_N+i]*dphi[i][k]; } } } for (int iqp = 0; iqp< N_QP; ++iqp) { double phi[N_SF]; shape_functions(QP[iqp],phi); for (int i = 0; i < n_fields; ++i) { s[i] = 0; sold[i] = 0; for (int j = 0; j < D; ++j) { ds[i*D+j] = 0; } } double c = 0, a = 0; double cold = 0, aold = 0; double dc[D] = {0}, da[D]={0}; for (int i = 0; i < N_SF; ++i) { c += problem->porosity[el[i]]*phi[i]; cold += problem->oldporosity[el[i]]*phi[i]; ... ... @@ -613,17 +611,9 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o double dofold = solution_old[el[i]*n_fields+j]; s[j] += phi[i]*dof; sold[j] += phi[i]*dofold; for (int k = 0; k < D; ++k) ds[j*D+k] += dphi[i][k]*dof; } for (int k=0; kporosity[el[i]]*dphi[i][k]; } if (problem->n_fluids==2) { a += problem->concentration[iel*N_N+i]*phi[i]; for (int k=0; kconcentration[iel*N_N+i]*dphi[i][k]; } } } double mu = problem->mu[0]; ... ... @@ -657,8 +647,6 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o fluid_problem_f(problem,h0+jd*n_fields,h1+jd*D*n_fields,s,ds,sold,c,dc,cold,rho,drho,mu,dt,iel,reduced_gravity,stab_param); ds[jfield*D+jd] = store; } int N2 = n_fields*N_SF; #define LOCAL_MATRIX2(a,b) local_matrix[IDX+N2*(a)+(b)] for (int ifield = 0; ifield < n_fields; ++ifield){ for (int iphi = 0; iphi < N_SF; ++iphi){ for (int jphi = 0; jphi < N_SF; ++jphi){ ... ... @@ -676,8 +664,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o double df01 = (h0[jd*n_fields+ifield]-f0[ifield])/deps; m += jw*phi[iphi]*dphi[jphi][jd]*df01; } int IDX = (iphi*N2)*n_fields + jphi*n_fields; LOCAL_MATRIX2(ifield,jfield) += m; LOCAL_MATRIX(iphi,jphi,ifield,jfield) += m; } } } ... ... @@ -1311,8 +1298,10 @@ static void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) { free(problem->solution); int n_fields = fluid_problem_n_fields(problem); problem->solution = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double)); free(problem->concentration); problem->concentration = (double*)malloc(mesh->n_elements*N_N*sizeof(double)); if (problem->n_fluids == 2) { free(problem->concentration); problem->concentration = (double*)malloc(mesh->n_elements*N_N*sizeof(double)); } for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){ problem->solution[i] = 0.; } ... ...
 ... ... @@ -73,10 +73,10 @@ void fluid_problem_advance_concentration(FluidProblem *problem, double dt) { if (iel1 < 0) continue; const int *cl0 = elbnd[interface[1]]; const int *cl1 = elbnd[interface[3]]; double *x[D]; double *x[N_LSF]; const int *el0 = &mesh->elements[iel0*N_N]; const int *el1 = &mesh->elements[iel1*N_N]; for (int i = 0; i < D; ++i){ for (int i = 0; i < N_LSF; ++i){ x[i] = &mesh->x[el0[cl0[i]]*3]; } double detbnd, n[D]; ... ... @@ -85,7 +85,7 @@ void fluid_problem_advance_concentration(FluidProblem *problem, double dt) { double a0 = 0, a1 = 0, un = 0; double phil[N_LQP]; l_shape_functions(LQP[iqp],phil); for (int i = 0; i < D; ++i) { for (int i = 0; i < N_LSF; ++i) { a0 += problem->concentration[iel0*N_N+cl0[i]]*phil[i]; a1 += problem->concentration[iel1*N_N+cl1[i]]*phil[i]; for (int id = 0; id < D; ++id) { ... ... @@ -96,7 +96,7 @@ void fluid_problem_advance_concentration(FluidProblem *problem, double dt) { const double r0 = -un*(a-a0); const double r1 = un*(a-a1); const double jw = LQW[iqp]*detbnd; for (int i = 0; i < D; ++i){ for (int i = 0; i < N_LSF; ++i){ rhs[iel0*N_N+cl0[i]] += phil[i]*r0*jw; rhs[iel1*N_N+cl1[i]] += phil[i]*r1*jw; } ... ... @@ -119,7 +119,7 @@ void fluid_problem_advance_concentration(FluidProblem *problem, double dt) { const int iel = interface[0]; const int *cl = elbnd[interface[1]]; const int *el = &mesh->elements[iel*N_N]; double *x[D]; double *x[N_LSF]; for (int i = 0; i < D; ++i){ x[i] = &mesh->x[el[cl[i]]*3]; } ... ... @@ -130,7 +130,7 @@ void fluid_problem_advance_concentration(FluidProblem *problem, double dt) { double phil[N_LQP]; l_shape_functions(LQP[iqp],phil); double a = 0; for (int i = 0; i < D; ++i) { for (int i = 0; i < N_LSF; ++i) { a += problem->concentration[iel*N_N+cl[i]]*phil[i]; } double un = 0; ... ... @@ -139,14 +139,14 @@ void fluid_problem_advance_concentration(FluidProblem *problem, double dt) { un += n[id]*dataqp[wbnd->vid+id]; } else { for (int i = 0; i < D; ++i) { for (int i = 0; i < N_LSF; ++i) { un += phil[i]*problem->solution[el[cl[i]]*n_fields+id]*n[id]; } } } const double jw = LQW[iqp]*detbnd; double f0 = -un * (un > 0 ? a : dataqp[wbnd->aid]); for (int iphi = 0; iphi < D; ++iphi){ for (int iphi = 0; iphi < N_LSF; ++iphi){ rhs[iel*N_N+cl[iphi]] += phil[iphi]*f0*jw*dt; } } ... ...
 ... ... @@ -8,6 +8,7 @@ #define N_N 3 #define N_QP 3 #define N_LQP 2 #define N_LSF 2 static double LQP[][1] = {{0.21132486540518713}, {0.7886751345948129}}; static double LQW[] = {0.5,0.5}; ... ... @@ -40,7 +41,6 @@ static double invDD(const double m[2][2], double inv[2][2]) inv[1][1] = m[0][0]/d; return d; } static int elbnd[][2] = {{0,1},{1,2},{2,0},{1,0},{2,1},{0,2}}; static void l_shape_functions(double *uvw, double *sf) { sf[0] = 1-uvw[0]; sf[1] = uvw[0]; ... ... @@ -63,6 +63,7 @@ static void get_normal_and_det(double *x[2], double n[2], double *det){ #define N_N 4 #define N_QP 5 #define N_LQP 3 #define N_LSF 3 static double QP[][3] = { {0.25, 0.25, 0.25}, {0.5, 1./6, 1./6}, ... ... @@ -109,7 +110,6 @@ static double invDD(const double m[3][3], double inv[3][3]) { return d; }; static int elbnd[4][3] = {{0,1,2},{0,2,3},{0,3,1},{1,3,2}}; static void l_shape_functions(double *uvw, double *sf) { sf[0] = 1-uvw[0]-uvw[1]; sf[1] = uvw[0]; ... ... @@ -133,5 +133,6 @@ static void get_normal_and_det(double *x[3], double n[3], double *det) { n[2] = -nn[2]/detbnd; *det = detbnd; } #endif
 ... ... @@ -27,22 +27,27 @@ #include #include "tools.h" #if DIMENSION == 2 static int elbnd[][2] = {{0,1},{1,2},{2,0},{1,0},{2,1},{0,2}}; #else static int elbnd[4][3] = {{0,1,2},{0,2,3},{0,3,1},{1,3,2}}; #endif typedef struct { int tag; int node[DIMENSION]; } BoundaryElement; int edgecmp(const void *pe0, const void *pe1) { int *e0 = (int*)pe0; int *e1 = (int*)pe1; typedef struct { int n[DIMENSION]; int element; int pos; // or tag if element == -1 }HalfEdge; int hedgecmp(const void *pe0, const void *pe1) { int *e0 = ((HalfEdge*)pe0)->n; int *e1 = ((HalfEdge*)pe1)->n; if (e0[0] != e1[0]) return e0[0]-e1[0]; #if DIMENSION==2 return e0[1]-e1[1]; #else if (e0[1] != e1[1]) return e0[1]-e1[1]; return e0[2]-e1[2]; #endif } static void swapint(int *a, int *b) { ... ... @@ -51,82 +56,86 @@ static void swapint(int *a, int *b) { *b = c; } static void sort_edge_nodes(int *n) { if(n[0]>n[1]) swapint(n,n+1); #if DIMENSION == 3 if(n[0]>n[2]) swapint(n,n+2); if(n[1]>n[2]) swapint(n+1,n+2); #endif } static void mesh_gen_edges(Mesh *m, int n_bnd_elements, BoundaryElement *bnd_elements) { #if DIMENSION == 2 int n_half_edges = m->n_elements*3+n_bnd_elements; int *halfedges = malloc(sizeof(int)*n_half_edges*4); //n0 n1 {(el,closure),(-1,tag)} int edgemap[3][2] = {{0,1},{1,2},{2,0}}; int n_half_edges = m->n_elements*(DIMENSION+1)+n_bnd_elements; HalfEdge *halfedges = malloc(sizeof(HalfEdge)*n_half_edges); HalfEdge *he = halfedges; // current int id = 0; for (int i=0; in_elements; ++i) { for (int j = 0; j < 3; ++j) { int n0 = m->elements[i*3+edgemap[j][0]]; int n1 = m->elements[i*3+edgemap[j][1]]; if (n0 < n1) { halfedges[id*4+0] = n0; halfedges[id*4+1] = n1; halfedges[id*4+3] = j; for (int j = 0; j < DIMENSION+1; ++j) { for (int k = 0; k < DIMENSION; ++k) { he->n[k] = m->elements[i*(DIMENSION+1)+elbnd[j][k]]; } else { halfedges[id*4+0] = n1; halfedges[id*4+1] = n0; halfedges[id*4+3] = j+3; } halfedges[id*4+2] = i; id++; sort_edge_nodes(he->n); he->element = i; he->pos = j; he++; } } for (int i=0; inode[0] < e->node[1]) { halfedges[id*4+0] = e->node[0]; halfedges[id*4+1] = e->node[1]; } else { halfedges[id*4+0] = e->node[1]; halfedges[id*4+1] = e->node[0]; for (int k = 0; k < DIMENSION; ++k){ he->n[k] = e->node[k]; } halfedges[id*4+2] = -1; halfedges[id*4+3] = e->tag; id ++; sort_edge_nodes(he->n); he->element = -1; he->pos = e->tag; he ++; } qsort(halfedges,n_half_edges,4*sizeof(int),edgecmp); qsort(halfedges,n_half_edges,sizeof(HalfEdge),hedgecmp); int n_edges = 0; for (int i = 0; i < n_half_edges; ++i) { if (i == 0 || edgecmp(halfedges+4*(i-1),halfedges+4*i)) n_edges++; if (i == 0 || hedgecmp(halfedges+i-1,halfedges+i)) n_edges++; } int *edges = malloc(n_edges*4*sizeof(int)); n_edges = 0; for (int i = 0; i < n_half_edges; ++i) { if (i == 0 || edgecmp(halfedges+4*(i-1),halfedges+4*i)){ if (i == 0 || hedgecmp(halfedges+i-1,halfedges+i)){ n_edges += 1; edges[(n_edges-1)*4+0] = halfedges[4*i+2]; edges[(n_edges-1)*4+1] = halfedges[4*i+3]; edges[(n_edges-1)*4+0] = halfedges[i].element; edges[(n_edges-1)*4+1] = halfedges[i].pos; edges[(n_edges-1)*4+2] = -1; edges[(n_edges-1)*4+3] = -1; } else { edges[(n_edges-1)*4+2] = halfedges[4*i+2]; edges[(n_edges-1)*4+3] = halfedges[4*i+3]; edges[(n_edges-1)*4+2] = halfedges[i].element; edges[(n_edges-1)*4+3] = halfedges[i].pos; } } m->n_interfaces = n_edges; m->interfaces = edges; free(halfedges); for (int i = 0; i < n_edges; ++i) { int *e = edges + i*4; // ensure boundary tag is on the second part of the closure if (edges[i*4+0] < 0) { swapint(edges+i*4,edges+i*4+2); swapint(edges+i*4+1,edges+i*4+3); if (e[0] < 0) { swapint(&e[0],&e[2]); swapint(&e[1],&e[3]); } // ensure edge is alligned with first element if (edges[i*4+1] >= 3) edges[i*4+1] -= 3; if (edges[i*4+2] >= 0 && edges[i*4+3] < 3) edges[i*4+3] += 3; } // compute the closure id of the second element if (e[2] >= 0) { #if DIMENSION==2 e[3] += 3; #else // todo Jon int k = 0; const int *e0 = &m->elements[e[0]*4]; const int *e1 = &m->elements[e[2]*4]; e[3] += 12; while (e0[elbnd[e[1]][0]] != e1[elbnd[e[3]][k]]) e[3] += 4; #endif } } } static int check_word_in_file(FILE *f, const char *w) { ... ...
 ... ... @@ -58,5 +58,21 @@ typedef struct { MeshBoundary *mesh_boundaries_create(const Mesh *m); void mesh_boundaries_free(MeshBoundary *mbnd, int n); #if DIMENSION == 2 static int elbnd[6][2] = { {0,1},{1,2},{2,0}, {1,0},{2,1},{0,2}}; #else static int elbnd[24][3] = { {0,1,2},{0,2,3},{0,3,1},{1,3,2}, {2,0,1},{3,0,2},{1,0,3},{2,1,3}, {1,2,0},{2,3,0},{3,1,0},{3,2,1}, {0,2,1},{0,3,2},{0,1,3},{1,2,3}, {2,1,0},{3,2,0},{1,3,0},{2,3,1}, {1,0,2},{2,0,3},{3,0,1},{3,1,2} }; #endif #endif
 ... ... @@ -50,11 +50,13 @@ dt = .01 # time step outf = 1 # number of iterations between output files def outerBndV(x) : a = np.where(np.abs(x[:,0]-0.111) < 0.01, 0.4*max(np.sin(t*np.pi*2./1),0),0) l2 = x[:,0]**2+x[:,2]**2 a = np.where(l2 < 0.03**2, 0.4*max(np.sin(t*np.pi*2./1),0),0) return a fluid = fluid.FluidProblem(3,g,[nu0*rho0,nu1*rho1],[rho0,rho1]) #fluid = fluid.FluidProblem(3,g,[nu1*rho1],[rho1]) fluid.load_msh("mesh.msh") fluid.set_wall_boundary("XVel") fluid.set_wall_boundary("ZVel") ... ... @@ -65,6 +67,7 @@ fluid.set_open_boundary("Bottom",velocity=[0,outerBndV,0],concentration=1)