/** # Fictitious-domain method with distributed Lagrangian multipliers We are looking for a numerical approximation for the solution of the following equations $$ \rho\left(\partial_t{\mathbf{u}} +\mathbf{u}\cdot \mathbf{\nabla}\mathbf{u}\right) = -\mathbf{\nabla} p + \mathbf{\nabla} \cdot\left(2\mu\mathbf{D}\right) + \rho\mathbf{a} + \mathbf{\lambda}~\text{over}~\Omega $$ $$ \left(1-\frac{\rho}{\rho_s}\right)\left(M\left(\partial_t{\mathbf{U}}-\mathbf{g}\right)\right)=-\int_{P(t)} {\mathbf{\lambda}} dv~\text{over}~P(t) $$ $$ \left(1-\frac{\rho}{\rho_s}\right)\left(\mathbf{I}\partial_t{\mathbf{\omega}} + \mathbf{\omega} \times\mathbf{I}\mathbf{\omega}\right) = -\int_{P(t)}\mathbf{r}\times{\mathbf{\lambda}} dv~\text{over}~P(t) $$ $$ \mathbf{u}-\left(\mathbf{U}+\mathbf{\omega}\times \mathbf{r}\right)=0~\text{over}~P(t) $$ $$ \mathbf{\nabla}\cdot\mathbf{u} = \mathbf{0} ~\text{over}~\Omega $$ From a numerical point of view, this is a complicated set of equations to solve direclty. One has to deal with three main difficulties: * an advection-diffusion equation * the imcompressibility constraint with the pressure $p$ as unknown * the particle's solid body's constraint with the Lagrange multipliers as unknow. The classic strategy is to split the problem into subproblems and solve them successively. We chose here a two-steps time-spliting. */ # define BGHOSTS 2 # define BSIZE 128 # ifndef DLM_Moving_particle # define DLM_Moving_particle 0 # define TRANSLATION 0 # define ROTATION 0 # endif # if DLM_Moving_particle # ifndef TRANSLATION # define TRANSLATION 1 # endif # ifndef ROTATION # define ROTATION 1 # endif # endif # ifndef NPARTICLES # define NPARTICLES 1 # endif # ifndef DLM_alpha_coupling # define DLM_alpha_coupling 0 # endif #ifndef debugBD # define debugBD 0 // set 1 to desactivate boundary points #endif #ifndef debugInterior # define debugInterior 0 // set 1 to desactivate interior points #endif /** Functions and structures needed for the implementation of fictitious domain method. */ # include "save_data.h" # include "dlmfd_functions.h" /** Basilisk scalars and vectors needed for the implementation of the fictitious domain method. */ vector DLM_lambda[]; scalar flagfield[]; scalar flagfield_mailleur[]; vector index_lambda[]; vector DLM_periodic_shift[]; static FILE * pdata[NPARTICLES]; static FILE * fdata[NPARTICLES]; particle particles[NPARTICLES] = {0}; #if DLM_alpha_coupling vector DLM_explicit[]; #endif /** ## First sub-problem: Navier Stokes problem We are using the centred solver for this problem. $$ \rho\left(\partial_t{\mathbf{u}}+\mathbf{u}\cdot {\mathbf{\nabla}}\mathbf{u}\right) = -{\mathbf{\nabla}}p + {\mathbf{\nabla}}\cdot\left(2\mu\mathbf{D}\right) + \rho\mathbf{a}~\text{over}~\Omega $$ $$ {\mathbf{\nabla}}\cdot\mathbf{u} = \mathbf{0} ~\text{over}~\Omega. $$ */ # include "navier-stokes/centered.h" /* Adding this small macro because dv() breaks the compability with embed.h */ #if dimension == 1 # define dlmfd_dv() (Delta) #elif dimension == 2 # define dlmfd_dv() (sq(Delta)) #else // dimension == 3 # define dlmfd_dv() (cube(Delta)) #endif /** ## Second sub-problem: Fictitious-domain problem We are solving a fictitious problem given by: $$ \rho \partial_t\mathbf{u} = {\mathbf{\lambda}}~\text{over}~\Omega $$ $$ \left(1-\frac{\rho}{\rho_s}\right)\left(M\left(\partial_t\mathbf{U}-\mathbf{g}\right)\right)=-\int_{P(t)}{\mathbf{\lambda}} dv~\text{over}~P(t) $$ $$ \left(1-\frac{\rho}{\rho_s}\right)\left(\mathbf{I}\partial_t {\mathbf{\omega}} + {\mathbf{\omega}} \times\mathbf{I}{\mathbf{\omega}}\right) = -\int_{P(t)}\mathbf{r}\times{\mathbf{\lambda}} dv~\text{over}~P(t) $$ $$ \mathbf{u}-\left(\mathbf{U}+{\mathbf{\omega}}\times \mathbf{r}\right)=0~\text{over}~P(t), $$ with unknowns $\mathbf{u}, \mathbf{U}, {\mathbf{\omega}}$ and ${\mathbf{\lambda}}$ being respectively the fluid velocity field, the particle's translational velocity, the particle's rotational velocity and the Lagrange multipliers. The particle occupies the domain $P(t)$ with density $\rho_s$, Inertia tensor $\mathbf{I}$ and mass $M$. The vector $\mathbf{g}$ is the gravity acceleration here. This leads to a saddle-point problem which is solved with an iterative solver (Uzawa, conjugate gradient algorithm). It is implemented in the function below. */ /* Timers and Timings */ timing dlmfd_globaltiming = {0.}; # include "dlmfd_perf.h" void DLMFD_subproblem (particle * p, const int i, const double rho_f) { /* Timers and Timings */ timer dlmfd_timer = timer_start (); double mpitimings[npe()]; double DLM_alpha = 0., DLM_beta = 0.; double DLM_tol = 1.e-5, DLM_nr2 = 0., DLM_nr2_km1 = 0., DLM_wv = 0.; int ki = 0; int DLM_maxiter = 100; vector qu[], tu[]; coord ppshift = {0, 0, 0}; static FILE * converge = fopen ("converge-uzawa", "w"); static FILE * cellvstime = fopen ("dlmfd-cells", "w"); /* Allocate and initialize particles */ allocate_and_init_particles (p, NPARTICLES, index_lambda, flagfield, flagfield_mailleur, DLM_periodic_shift); #if debugInterior == 0 Cache * Interior[NPARTICLES]; #endif #if debugBD == 0 Cache * Boundary[NPARTICLES]; SolidBodyBoundary * sbm[NPARTICLES]; #endif GeomParameter * gci[NPARTICLES]; #if DLM_Moving_particle #if TRANSLATION coord * qU[NPARTICLES]; coord * U[NPARTICLES]; coord * tU[NPARTICLES]; #endif #if ROTATION coord * qw[NPARTICLES]; coord * w[NPARTICLES]; coord * tw[NPARTICLES]; #endif #endif for (int k = 0; k < NPARTICLES; k++) { gci[k] = &(p[k].g); #if debugInterior == 0 Interior[k] = &(p[k].Interior); #endif #if debugBD == 0 Boundary[k] = &(p[k].reduced_domain); sbm[k] = &(p[k].s); #endif #if DLM_Moving_particle #if TRANSLATION qU[k] = &(p[k].qU); U[k] = &(p[k].U); tU[k] = &(p[k].tU); #endif #if ROTATION qw[k] = &(p[k].qw); w[k] = &(p[k].w); tw[k] = &(p[k].tw); #endif #endif } if (NPARTICLES > 1) remove_too_close_multipliers (p, index_lambda); #if debugBD == 0 #if DLM_Moving_particle reverse_fill_flagfield (p, flagfield_mailleur, index_lambda, 0); #endif /* Consider fictitious domain's boundary points only if they are far enough of the domain's solid boundaries (i.e a distance greater than 2*Delta, with Delta being the smallest cell size) otherwise the fictitious domain is pixalated. This does not apply when the box is periodic */ double twodelta; # if dimension == 2 if (!Period.x && !Period.y) { #elif dimension == 3 if (!Period.x && !Period.y && !Period.z) { #endif foreach_level (depth()) { twodelta = 2.*Delta; #if dimension == 2 if (( x > (L0 - twodelta + X0)) || (x < (twodelta + X0)) || ( y > (L0 - twodelta + Y0)) || (y < (twodelta + Y0))) #elif dimension == 3 if ((x > (L0 - twodelta + X0)) || (x < (twodelta + X0)) || (y > (L0 - twodelta + Y0)) || (y < (twodelta + Y0)) || (z > (L0 - twodelta +Z0)) || (z < (twodelta + Z0))) #endif { if (index_lambda.x[] > -1.) { index_lambda.x[] = -1.; } } } } boundary((scalar*) {index_lambda}); /* Tagg the cells belonging to the stencil of the boundary-multipliers points */ reverse_fill_flagfield (p, flagfield, index_lambda, 1); #endif /* Alllocate arrays for CG/Uzawa. */ vector DLM_r[]; vector DLM_w[]; vector DLM_v[]; /* Create fictitious-domain's cache for the interior domain. */ #if debugInterior == 0 for (int k = 0; k < NPARTICLES; k++) { if (!(p[k].iswall) && !(p[k].iscube)) { create_FD_Interior_Particle (&p[k], index_lambda, DLM_periodic_shift, flagfield); } if (p[k].iswall) { create_FD_Interior_Wall(&p[k], index_lambda, flagfield); } if (p[k].iscube) { create_FD_Interior_Cube_v2 (&p[k], index_lambda, DLM_periodic_shift); } } #endif boundary ((scalar*) {DLM_periodic_shift}); int allpts = 0; int lm = 0; int tcells = 0; foreach() { tcells ++; DLM_lambda.x[] = 0.; DLM_r.x[] = 0.; DLM_w.x[] = 0.; DLM_v.x[] = 0.; qu.x[] = 0.; tu.x[] = 0.; DLM_lambda.y[] = 0.; DLM_r.y[] = 0.; DLM_w.y[] = 0.; DLM_v.y[] = 0.; qu.y[] = 0.; tu.y[] = 0.; DLM_lambda.z[] = 0.; DLM_r.z[] = 0.; DLM_w.z[] = 0.; DLM_v.z[] = 0.; qu.z[] = 0.; tu.z[] = 0.; } #if _MPI mpi_all_reduce (tcells, MPI_INTEGER, MPI_SUM); #endif lm = total_dlmfd_multipliers (p, NPARTICLES); /* Get track of the number of multipliers for statistics */ for (int k = 0; k < NPARTICLES; k++) p[k].tmultipliers += lm; allpts = total_dlmfd_cells (p, NPARTICLES); /* Get track of the number of cells involved in the dlmfd solver for statistics */ for (int k = 0; k < NPARTICLES; k++) p[k].tcells += allpts; fprintf (stderr,"# Total Lagrange multipliers: %d, constraint cells: %d, Total cells: %d \n", lm, allpts, tcells); if (i == 0) fprintf (cellvstime,"# Time Iteration \t Lagrange multipliers \t Constraint cells \t Total cells\n"); fprintf (cellvstime,"%d \t %d \t %d \t %d \n", i, lm, allpts, tcells); /* fflush (cellvstime); */ #if DLM_Moving_particle #if TRANSLATION for (int k = 0; k < NPARTICLES; k++) { (*qU[k]).x = 0.; (*tU[k]).x = 0.; (*qU[k]).y = 0.; (*tU[k]).y = 0.; (*qU[k]).z = 0.; (*tU[k]).z = 0.; } #endif #if ROTATION for (int k = 0; k < NPARTICLES; k++) { (*qw[k]).x = 0.; (*tw[k]).x = 0.; (*qw[k]).y = 0.; (*tw[k]).y = 0.; (*qw[k]).z = 0.; (*tw[k]).z = 0.; } #endif #endif /** # Iterative Uzawa/conjugate-gradient solver Implementation of the iterative Uzawa/conjugate gradient solver for the fictitious domain problem. The code below follows closely the steps given in the annex of the paper "Wachs et al J.Eng. Math, 2011". (minus all the sign errors)*/ /* Iterative CG/Uzawa */ /* Initialization */ /* ------------- */ /* Compute qu = fu - (M_u^T)*lambda^0 with fu = rho*dV/dt*u^n+1/2 + alpha*lambda^n. */ /* Note that (M_u^T)*lambda <=> _P(t) (v is the test function for u)*/ /* For moving particles add: */ /* qU = fU -(M_U^T)*lambda^0 with fU = (1 - rho_f / rho_s)*M*(U^{n+1/2}/dt + g) with M the particle's mass and g the gravity. */ /* Note that (M_U^T)*lambda <=> - _P(t) (V is the test function for U)*/ /* qw = fw -(M_w^T)*lambda^0 with fw=(1-rho_f/rho_s)/dt*Ip*w^n+1/2 * with Ip the Inertia tensor.*/ /* Note that (M_w^T)*lambda <=> - _P(t) (xi is the test * function for w)*/ foreach() { foreach_dimension() { qu.x[] = rho_f*dlmfd_dv()*u.x[]/dt; #if DLM_alpha_coupling qu.x[] += DLM_explicit.x[]; DLM_explicit.x[] = 0.; #endif } } boundary((scalar*){qu}); /* Interior points qu = fu -(M_u^T)*lambda^0 */ /* qu = fu -(M_u^T)*lambda^0 = fu -_P(t) */ /* Interior points qU = fU -(M_U^T)*lambda^0 */ /* qU = fU -(M_U^T)*lambda^0 = fU + _P(t) */ /* Interior points qw = fw -(M_w^T)*lambda^0 */ /* qw = fw -(M_w^T)*lambda^0 = fw + _P(t) */ /* For moving particles the fU and fw parts are added after the scalar product (for mpi purpose) */ #if debugInterior == 0 for (int k = 0; k < NPARTICLES; k++) { foreach_cache((*Interior[k])) { if ((flagfield[] < 1) && ((int)index_lambda.y[] == k)) { foreach_dimension() { qu.x[] -= DLM_lambda.x[]; #if DLM_Moving_particle #if TRANSLATION (*qU[k]).x += DLM_lambda.x[]; #endif #endif } #if DLM_Moving_particle #if ROTATION /* Modify temporerly the particle center position for periodic boundary condition */ foreach_dimension() (*gci[k]).center.x += DLM_periodic_shift.x[]; /* as a.(b^c) = c.(a^b) = b.(c^a) */ /* _P(t) = _P(t) = _P(t) */ (*qw[k]).x += (DLM_lambda.z[]*(y - (*gci[k]).center.y) - DLM_lambda.y[]*(z - (*gci[k]).center.z)); // lambda_z*r_y - lambda_y*r_z (*qw[k]).y += (DLM_lambda.x[]*(z - (*gci[k]).center.z) - DLM_lambda.z[]*(x - (*gci[k]).center.x)); // lambda_x*r_z - lambda_z*r_x (*qw[k]).z += (DLM_lambda.y[]*(x - (*gci[k]).center.x) - DLM_lambda.x[]*(y - (*gci[k]).center.y)); // lambda_y*r_x - lambda_x*r_y foreach_dimension() (*gci[k]).center.x -= DLM_periodic_shift.x[]; #endif #endif } } } #endif /* Boundary points qu = fu -(M_u^T)*lambda^0 */ /* Boundary points qU = fU -(M_U^T)*lambda^0 */ /* Boundary points qw = fw -(M_w^T)*lambda^0 */ # if debugBD == 0 double weight = 0.; coord weightcellpos = {0, 0, 0}; coord lambdacellpos = {0, 0, 0}; coord lambdapos = {0, 0, 0}; coord sum = {0, 0, 0}; boundary ((scalar*) {DLM_lambda}); for (int k = 0; k < NPARTICLES; k++) { particle * pp = &p[k]; foreach_cache ((*Boundary[k])) { sum.x = 0; sum.y = 0; sum.z = 0; weightcellpos.x = x; weightcellpos.y = y; weightcellpos.z = 0.; #if dimension == 3 weightcellpos.z = z; #endif foreach_neighbor() { if (((int)index_lambda.x[] > -1) && (level == depth()) && is_leaf(cell) && ((int)index_lambda.y[] == pp->pnum)) { lambdacellpos.x = x; lambdacellpos.y = y; lambdacellpos.z = 0.; lambdapos.x = (*sbm[k]).x[(int)index_lambda.x[]]; lambdapos.y = (*sbm[k]).y[(int)index_lambda.x[]]; lambdapos.z = 0.; #if dimension == 3 lambdacellpos.z = z; lambdapos.z = (*sbm[k]).z[(int)index_lambda.x[]]; #endif foreach_dimension() { ppshift.x = 0.; if (lambdacellpos.x < 0.) ppshift.x = -L0; if (lambdacellpos.x > L0) ppshift.x = +L0; } weight = reversed_weight (pp, weightcellpos, lambdacellpos, lambdapos, Delta, ppshift,1,1); foreach_dimension() sum.x += weight*DLM_lambda.x[]; } } qu.x[] -= sum.x; qu.y[] -= sum.y; qu.z[] -= sum.z; #if DLM_Moving_particle if (index_lambda.x[] > -1 && ((int)index_lambda.y[] == pp->pnum)) { lambdapos.x = (*sbm[k]).x[(int)index_lambda.x[]]; lambdapos.y = (*sbm[k]).y[(int)index_lambda.x[]]; lambdapos.z = 0.; #if dimension == 3 lambdapos.z = (*sbm[k]).z[(int)index_lambda.x[]]; #endif #if TRANSLATION foreach_dimension() (*qU[k]).x += DLM_lambda.x[]; #endif #if ROTATION /* Modify temporerly the particle center position for periodic boundary condition */ foreach_dimension() (*gci[k]).center.x += DLM_periodic_shift.x[]; /* as a.(b^c) = c.(a^b) = b.(c^a) */ /* _P(t) = _P(t) = _P(t) */ (*qw[k]).x += (DLM_lambda.z[]*(lambdapos.y - (*gci[k]).center.y) - DLM_lambda.y[]*(lambdapos.z - (*gci[k]).center.z)); //lambda_z*r_y - lambda_y*r_z (*qw[k]).y += (DLM_lambda.x[]*(lambdapos.z - (*gci[k]).center.z) - DLM_lambda.z[]*(lambdapos.x - (*gci[k]).center.x)); //lambda_x*r_z - lambda_z*r_x (*qw[k]).z += (DLM_lambda.y[]*(lambdapos.x - (*gci[k]).center.x) - DLM_lambda.x[]*(lambdapos.y - (*gci[k]).center.y)); //lambda_y*r_x - lambda_x*r_y foreach_dimension() (*gci[k]).center.x -= DLM_periodic_shift.x[]; #endif } #endif } } #endif /* Broadcast through mpi */ #if DLM_Moving_particle #if _MPI for (int k = 0; k < NPARTICLES; k++) { #if TRANSLATION mpi_all_reduce ((*qU[k]).x, MPI_DOUBLE, MPI_SUM); mpi_all_reduce ((*qU[k]).y, MPI_DOUBLE, MPI_SUM); mpi_all_reduce ((*qU[k]).z, MPI_DOUBLE, MPI_SUM); #endif #if ROTATION mpi_all_reduce ((*qw[k]).x, MPI_DOUBLE, MPI_SUM); mpi_all_reduce ((*qw[k]).y, MPI_DOUBLE, MPI_SUM); mpi_all_reduce ((*qw[k]).z, MPI_DOUBLE, MPI_SUM); #endif } #endif /* Adding here the fU and fw parts to qU and qw */ for (int k = 0; k < NPARTICLES; k++) { #if TRANSLATION /* (*qU[k]).x += (1. - (rho_f/p[k].rho_s))*(p[k].M)*(p[k].gravity.x + p[k].adforce.x + (*U[k]).x/dt) ; */ /* (*qU[k]).y += (1. - (rho_f/p[k].rho_s))*(p[k].M)*(p[k].gravity.y + p[k].adforce.y + (*U[k]).y/dt) ; */ /* (*qU[k]).z += (1. - (rho_f/p[k].rho_s))*(p[k].M)*(p[k].gravity.z + p[k].adforce.z + (*U[k]).z/dt) ; */ /* Modification for the explicit added term */ (*qU[k]).x += (1. - (rho_f/p[k].rho_s))*(p[k].M)*(p[k].gravity.x + p[k].adforce.x ) + p[k].DLMFD_couplingfactor * p[k].M * (*U[k]).x / dt ; (*qU[k]).y += (1. - (rho_f/p[k].rho_s))*(p[k].M)*(p[k].gravity.y + p[k].adforce.y ) + p[k].DLMFD_couplingfactor * p[k].M * (*U[k]).y / dt ; (*qU[k]).z += (1. - (rho_f/p[k].rho_s))*(p[k].M)*(p[k].gravity.z + p[k].adforce.z ) + p[k].DLMFD_couplingfactor * p[k].M * (*U[k]).z / dt ; #endif #if ROTATION /* The inertia tensor is */ /* Ixx -Ixy -Ixz */ /* -Iyx Iyy -Iyz */ /* -Izx -Izy Izz */ /* with */ /* Ip[0] = Ixx */ /* Ip[1] = Iyy */ /* Ip[2] = Izz */ /* Ip[3] = Ixy */ /* Ip[4] = Ixz */ /* Ip[5] = Iyz */ /* (*qw[k]).x += (1. - (rho_f/p[k].rho_s))*( (p[k].Ip[0])*(*w[k]).x - (p[k].Ip[3])*(*w[k]).y - (p[k].Ip[4])*(*w[k]).z)/dt; */ /* (*qw[k]).y += (1. - (rho_f/p[k].rho_s))*(-(p[k].Ip[3])*(*w[k]).x + (p[k].Ip[1])*(*w[k]).y - (p[k].Ip[5])*(*w[k]).z)/dt; */ /* (*qw[k]).z += (1. - (rho_f/p[k].rho_s))*(-(p[k].Ip[4])*(*w[k]).x - (p[k].Ip[5])*(*w[k]).y + (p[k].Ip[2])*(*w[k]).z)/dt; */ /* Modification for the explicit added term */ (*qw[k]).x += p[k].DLMFD_couplingfactor * ( (p[k].Ip[0])*(*w[k]).x - (p[k].Ip[3])*(*w[k]).y - (p[k].Ip[4])*(*w[k]).z) / dt; (*qw[k]).y += p[k].DLMFD_couplingfactor * (-(p[k].Ip[3])*(*w[k]).x + (p[k].Ip[1])*(*w[k]).y - (p[k].Ip[5])*(*w[k]).z) / dt; (*qw[k]).z += p[k].DLMFD_couplingfactor * (-(p[k].Ip[4])*(*w[k]).x - (p[k].Ip[5])*(*w[k]).y + (p[k].Ip[2])*(*w[k]).z) / dt; #endif } #endif /* Invert L*u^0 = qu with L=dv*rho/dt*Identity_Matrix */ foreach() { foreach_dimension() { u.x[] = qu.x[]*dt/(rho_f*dlmfd_dv()); } } /* Invert A_U U^0 = qU with A_U = ((1 - rho_f/rho_s)*M)/dt * Identity_Matrix */ /* Invert A_w w^0 = qw with A_w = ((1 - rho_f/rho_s)*Ip)/dt * Identity_Matrix */ #if DLM_Moving_particle for (int k = 0; k < NPARTICLES; k++) { #if TRANSLATION /* (*U[k]).x = (dt*(*qU[k]).x)/((1. - rho_f/p[k].rho_s)*p[k].M); */ /* (*U[k]).y = (dt*(*qU[k]).y)/((1. - rho_f/p[k].rho_s)*p[k].M); */ /* (*U[k]).z = (dt*(*qU[k]).z)/((1. - rho_f/p[k].rho_s)*p[k].M); */ /* Modification for the explicit added term */ (*U[k]).x = (dt*(*qU[k]).x)/(p[k].DLMFD_couplingfactor*p[k].M); (*U[k]).y = (dt*(*qU[k]).y)/(p[k].DLMFD_couplingfactor*p[k].M); (*U[k]).z = (dt*(*qU[k]).z)/(p[k].DLMFD_couplingfactor*p[k].M); #endif #if ROTATION /* The inertia tensor is */ /* Ixx -Ixy -Ixz */ /* -Iyx Iyy -Iyz */ /* -Izx -Izy Izz */ /* with */ /* Ip[0] = Ixx */ /* Ip[1] = Iyy */ /* Ip[2] = Izz */ /* Ip[3] = Ixy */ /* Ip[4] = Ixz */ /* Ip[5] = Iyz */ /* Allocate temporly a matrix for Ip */ double ** Imat = malloc(3 * sizeof(double*)); for (int ii = 0; ii < 3; ii++) { Imat[ii] = malloc(3 * sizeof(double)); } Imat[0][0] = p[k].Ip[0]; Imat[1][1] = p[k].Ip[1]; Imat[2][2] = p[k].Ip[2]; Imat[0][1] = -p[k].Ip[3]; Imat[0][2] = -p[k].Ip[4]; Imat[1][2] = -p[k].Ip[5]; Imat[1][0] = Imat[0][1]; Imat[2][0] = Imat[0][2]; Imat[2][1] = Imat[1][2]; /* A_w = ((1 - rho_f/rho_s)*Ip)/dt */ double Aw[3][3]; for (int ll = 0; ll < 3; ll++) { for (int nn = 0; nn < 3; nn++) { /* Aw[ll][nn] = (1 - rho_f/p[k].rho_s)*Imat[ll][nn]/dt; */ /* Modification for the explicit added term */ Aw[ll][nn] = p[k].DLMFD_couplingfactor * Imat[ll][nn] / dt; } } /* Takes in Aw and return Aw^-1 as Imat */ inverse3by3matrix(Aw, Imat); (*w[k]).x = Imat[0][0]*(*qw[k]).x + Imat[0][1]*(*qw[k]).y + Imat[0][2]*(*qw[k]).z; (*w[k]).y = Imat[1][0]*(*qw[k]).x + Imat[1][1]*(*qw[k]).y + Imat[1][2]*(*qw[k]).z; (*w[k]).z = Imat[2][0]*(*qw[k]).x + Imat[2][1]*(*qw[k]).y + Imat[2][2]*(*qw[k]).z; for (int ii = 0; ii < 3; ii++) free(Imat[ii]); free(Imat); #endif } #endif /* Compute residual r^0 = G - M_u*u^0 - M_U*U^0 - M_w*w^0 */ /* Note that M_u*u^0 <=> _P(t) */ /* Note that M_U*U^0 <=> - _P(t) */ /* Note that M_w*w^0 <=> - _P(t) */ /* alpha is he test function for lambda the Lagrange multipliers */ /* So r^0 = G -_P(t) + _P(t) + _P(t) */ /* In the 2D case when rotation is enabled, the z-component of the * residual becomes involved in the algorithm, this is why the * foreach_dimension() iterators are being removed from this * point. Otherwise the z components would be discarded */ /* Interior points: r^0 = G - M_u*u^0 - M_U*U^0 - M_w*w^0 */ /* So r^0 = G -_P(t) + _P(t) + _P(t) */ #if debugInterior == 0 for (int k = 0; k < NPARTICLES; k++) { #if !DLM_Moving_particle coord imposedU = (p[k]).imposedU; coord imposedw = (p[k]).imposedw; #endif foreach_cache((*Interior[k])) { if ((flagfield[] < 1) && ((int)index_lambda.y[] == k)) { foreach_dimension() { DLM_r.x[] = - u.x[]; #if DLM_Moving_particle #if TRANSLATION DLM_r.x[] += (*U[k]).x; #endif #else DLM_r.x[] += imposedU.x; #endif } #if DLM_Moving_particle #if ROTATION /* Modify temporerly the particle center position for periodic boundary condition */ foreach_dimension() (*gci[k]).center.x += DLM_periodic_shift.x[]; /* if (DLM_periodic_shift.x[] < 0) { */ /* printf ("position (%f,%f,%f), center shifted of %f\n", x, y, z, DLM_periodic_shift.x[]); */ /* printf ("new sphere position (%f,%f,%f)\n", (*gci[k]).center.x, (*gci[k]).center.y, (*gci[k]).center.z); */ /* } */ DLM_r.x[] += ((*w[k]).y*(z - (*gci[k]).center.z) - (*w[k]).z*(y - (*gci[k]).center.y)); // w_y*r_z - w_z*r_y DLM_r.y[] += ((*w[k]).z*(x - (*gci[k]).center.x) - (*w[k]).x*(z - (*gci[k]).center.z)); // w_z*r_x - w_x*r_z DLM_r.z[] += ((*w[k]).x*(y - (*gci[k]).center.y) - (*w[k]).y*(x - (*gci[k]).center.x)); // w_x*r_y - w_y*r_x foreach_dimension() (*gci[k]).center.x -= DLM_periodic_shift.x[]; #endif #else /* Modify temporerly the particle center position for periodic boundary condition */ foreach_dimension() (*gci[k]).center.x += DLM_periodic_shift.x[]; DLM_r.x[] += (imposedw.y*(z - (*gci[k]).center.z) - imposedw.z*(y - (*gci[k]).center.y)); DLM_r.y[] += (imposedw.z*(x - (*gci[k]).center.x) - imposedw.x*(z - (*gci[k]).center.z)); DLM_r.z[] += (imposedw.x*(y - (*gci[k]).center.y) - imposedw.y*(x - (*gci[k]).center.x)); foreach_dimension() (*gci[k]).center.x -= DLM_periodic_shift.x[]; #endif } } } #endif /* Boundary points: r^0 = G - M_u*u^0 - M_U*U^0 - M_w*w^0 */ /* So r^0 = G -_P(t) + _P(t) + _P(t) */ #if debugBD == 0 double testweight = 0.; boundary((scalar*){u}); for (int k = 0; k < NPARTICLES; k++) { particle * pp = &p[k]; #if !DLM_Moving_particle coord imposedU = (p[k]).imposedU; coord imposedw = (p[k]).imposedw; #endif foreach_cache((*Boundary[k])) { if (index_lambda.x[] > -1 && ((int)index_lambda.y[] == pp->pnum)) { lambdacellpos.x = x; lambdacellpos.y = y; lambdacellpos.z = 0.; lambdapos.x = (*sbm[k]).x[(int)index_lambda.x[]]; lambdapos.y = (*sbm[k]).y[(int)index_lambda.x[]]; lambdapos.z = 0.; #if dimension == 3 lambdacellpos.z = z; lambdapos.z = (*sbm[k]).z[(int)index_lambda.x[]]; #endif sum.x = 0; sum.y = 0; sum.z = 0; testweight = 0.; foreach_dimension() ppshift.x = DLM_periodic_shift.x[]; foreach_neighbor() { if ((level == depth())) { weightcellpos.x = x; weightcellpos.y = y; weightcellpos.z = 0.; #if dimension == 3 weightcellpos.z = z; #endif weight = reversed_weight (pp, weightcellpos, lambdacellpos, lambdapos, Delta, ppshift, 1., 0.); testweight += weight; foreach_dimension() sum.x += weight*u.x[]; } } if ((testweight < 1. - 0.000001) || (testweight > 1. + 0.000001)) printf("testweight = %f\n",testweight); DLM_r.x[] = -sum.x; DLM_r.y[] = -sum.y; DLM_r.z[] = -sum.z; #if DLM_Moving_particle #if TRANSLATION foreach_dimension() DLM_r.x[] += (*U[k]).x; #endif #if ROTATION /* Modify temporerly the particle center position for periodic boundary condition */ foreach_dimension() (*gci[k]).center.x += DLM_periodic_shift.x[]; DLM_r.x[] += ((*w[k]).y*(lambdapos.z - (*gci[k]).center.z) - (*w[k]).z*(lambdapos.y - (*gci[k]).center.y)); // w_y*r_z - w_z*r_y DLM_r.y[] += ((*w[k]).z*(lambdapos.x - (*gci[k]).center.x) - (*w[k]).x*(lambdapos.z - (*gci[k]).center.z)); // w_z*r_x - w_x*r_z DLM_r.z[] += ((*w[k]).x*(lambdapos.y - (*gci[k]).center.y) - (*w[k]).y*(lambdapos.x - (*gci[k]).center.x)); // w_x*r_y - w_y*r_x foreach_dimension() (*gci[k]).center.x -= DLM_periodic_shift.x[]; #endif #else foreach_dimension() { DLM_r.x[] += imposedU.x; } /* Modify temporerly the particle center position for periodic boundary condition */ foreach_dimension() (*gci[k]).center.x += DLM_periodic_shift.x[]; DLM_r.x[] += (imposedw.y*(lambdapos.z - (*gci[k]).center.z) - imposedw.z*(lambdapos.y - (*gci[k]).center.y)); DLM_r.y[] += (imposedw.z*(lambdapos.x - (*gci[k]).center.x) - imposedw.x*(lambdapos.z - (*gci[k]).center.z)); DLM_r.z[] += (imposedw.x*(lambdapos.y - (*gci[k]).center.y) - imposedw.y*(lambdapos.x - (*gci[k]).center.x)); foreach_dimension() (*gci[k]).center.x -= DLM_periodic_shift.x[]; #endif } } } #endif boundary((scalar*){DLM_r}); DLM_nr2 = 0.; /* set DLM_w = DLM_r */ foreach() { DLM_w.x[] = DLM_r.x[]; DLM_w.y[] = DLM_r.y[]; DLM_w.z[] = DLM_r.z[]; DLM_nr2 += sq(DLM_r.x[]) + sq(DLM_r.y[]) + sq(DLM_r.z[]); } #if _MPI mpi_all_reduce(DLM_nr2, MPI_DOUBLE, MPI_SUM); #endif /* Iterative loop */ /* -------------- */ for (ki = 1; ki < DLM_maxiter && (sqrt(DLM_nr2) > DLM_tol); ki++) { /* (1) Initialize and compute qu = (M^T)*DLM_w */ foreach() { qu.x[] = 0.; qu.y[] = 0.; qu.z[] = 0.; tu.x[] = 0.; tu.y[] = 0.; tu.z[] = 0.; } #if DLM_Moving_particle for (int k = 0; k < NPARTICLES; k++) { #if TRANSLATION (*qU[k]).x = 0.; (*tU[k]).x = 0.; (*qU[k]).y = 0.; (*tU[k]).y = 0.; (*qU[k]).z = 0.; (*tU[k]).z = 0.; #endif #if ROTATION (*qw[k]).x = 0.; (*tw[k]).x = 0.; (*qw[k]).y = 0.; (*tw[k]).y = 0.; (*qw[k]).z = 0.; (*tw[k]).z = 0.; #endif } #endif /* Interior points qu = (M_u^T)DLM_w = _P(t) */ /* Interior points qU = (M_U^T)DLM_w = - _P(t) */ /* Interior points qw = (M_w^T)DLM_w = - _P(t) */ /* - _P(t) = - _P(t) = - _P(t) */ #if debugInterior == 0 for (int k = 0; k < NPARTICLES; k++) { foreach_cache((*Interior[k])) { if ((flagfield[] < 1) && ((int)index_lambda.y[] == k)) { qu.x[] = DLM_w.x[]; qu.y[] = DLM_w.y[]; qu.z[] = DLM_w.z[]; #if DLM_Moving_particle #if TRANSLATION (*qU[k]).x += -DLM_w.x[]; (*qU[k]).y += -DLM_w.y[]; (*qU[k]).z += -DLM_w.z[]; #endif #endif #if DLM_Moving_particle #if ROTATION /* Modify temporerly the particle center position for periodic boundary condition */ foreach_dimension() (*gci[k]).center.x += DLM_periodic_shift.x[]; (*qw[k]).x += -(DLM_w.z[]*(y - (*gci[k]).center.y) - DLM_w.y[]*(z - (*gci[k]).center.z)); // -(w_z*r_y - w_y*r_z) (*qw[k]).y += -(DLM_w.x[]*(z - (*gci[k]).center.z) - DLM_w.z[]*(x - (*gci[k]).center.x)); // -(w_x*r_z - w_z*r_x) (*qw[k]).z += -(DLM_w.y[]*(x - (*gci[k]).center.x) - DLM_w.x[]*(y - (*gci[k]).center.y)); // -(w_y*r_x - w_x*r_y) foreach_dimension() (*gci[k]).center.x -= DLM_periodic_shift.x[]; #endif #endif } } } #endif /* Boundary points qu = (M_u^T)DLM_w = _P(t) */ /* Boundary points qU = (M_U^T)DLM_w = - _P(t) */ /* Boundary points qw = (M_w^T)DLM_w = - _P(t) */ /* - _P(t) = - _P(t) = - _P(t) */ # if debugBD == 0 boundary ((scalar*){DLM_w, qu}); for (int k = 0; k < NPARTICLES; k++) { particle * pp = &p[k]; foreach_cache((*Boundary[k])) { sum.x = 0; sum.y = 0; sum.z = 0; weightcellpos.x = x; weightcellpos.y = y; weightcellpos.z = 0.; #if dimension == 3 weightcellpos.z = z; #endif foreach_neighbor() { if(((int)index_lambda.x[] > -1) && (level == depth()) && (is_leaf(cell)) && ((int)index_lambda.y[] == pp->pnum)) { lambdacellpos.x = x; lambdacellpos.y = y; lambdacellpos.z = 0.; lambdapos.x = (*sbm[k]).x[(int)index_lambda.x[]]; lambdapos.y = (*sbm[k]).y[(int)index_lambda.x[]]; lambdapos.z = 0.; #if dimension == 3 lambdacellpos.z = z; lambdapos.z = (*sbm[k]).z[(int)index_lambda.x[]]; #endif foreach_dimension() { ppshift.x = 0.; if (lambdacellpos.x < 0.) ppshift.x = -L0; if (lambdacellpos.x > L0) ppshift.x = +L0; } weight = reversed_weight (pp, weightcellpos, lambdacellpos, lambdapos, Delta, ppshift, 1., 1.); sum.x += weight*DLM_w.x[]; sum.y += weight*DLM_w.y[]; sum.z += weight*DLM_w.z[]; } } qu.x[] += sum.x; qu.y[] += sum.y; qu.z[] += sum.z; //+= here as one fluid cell can be affected by multiples particle's boundary multipliers #if DLM_Moving_particle if (index_lambda.x[] > -1 && ((int)index_lambda.y[] == pp->pnum)) { lambdapos.x = (*sbm[k]).x[(int)index_lambda.x[]]; lambdapos.y = (*sbm[k]).y[(int)index_lambda.x[]]; lambdapos.z = 0.; #if dimension == 3 lambdapos.z = (*sbm[k]).z[(int)index_lambda.x[]]; #endif #if TRANSLATION (*qU[k]).x += -DLM_w.x[]; (*qU[k]).y += -DLM_w.y[]; (*qU[k]).z += -DLM_w.z[]; #endif #if ROTATION /* Modify temporerly the particle center position for periodic boundary condition */ foreach_dimension() (*gci[k]).center.x += DLM_periodic_shift.x[]; (*qw[k]).x += -(DLM_w.z[]*(lambdapos.y - (*gci[k]).center.y) - DLM_w.y[]*(lambdapos.z - (*gci[k]).center.z)); // -(w_z*r_y - w_y*r_z) (*qw[k]).y += -(DLM_w.x[]*(lambdapos.z - (*gci[k]).center.z) - DLM_w.z[]*(lambdapos.x - (*gci[k]).center.x)); // -(w_x*r_z - w_z*r_x) (*qw[k]).z += -(DLM_w.y[]*(lambdapos.x - (*gci[k]).center.x) - DLM_w.x[]*(lambdapos.y - (*gci[k]).center.y)); // -(w_y*r_x - w_x*r_y) foreach_dimension() (*gci[k]).center.x -= DLM_periodic_shift.x[]; #endif } #endif } } #endif #if DLM_Moving_particle #if _MPI for (int k = 0; k < NPARTICLES; k++) { #if TRANSLATION mpi_all_reduce ((*qU[k]).x, MPI_DOUBLE, MPI_SUM); mpi_all_reduce ((*qU[k]).y, MPI_DOUBLE, MPI_SUM); mpi_all_reduce ((*qU[k]).z, MPI_DOUBLE, MPI_SUM); #endif #if ROTATION mpi_all_reduce ((*qw[k]).x, MPI_DOUBLE, MPI_SUM); mpi_all_reduce ((*qw[k]).y, MPI_DOUBLE, MPI_SUM); mpi_all_reduce ((*qw[k]).z, MPI_DOUBLE, MPI_SUM); #endif } #endif #endif /* (2) Invert L*t = qu with L=rho*dV/dt*I */ foreach() { tu.x[] = qu.x[]*dt/(rho_f*dlmfd_dv()); tu.y[] = qu.y[]*dt/(rho_f*dlmfd_dv()); tu.z[] = qu.z[]*dt/(rho_f*dlmfd_dv()); } #if DLM_Moving_particle for (int k = 0; k < NPARTICLES; k++) { #if TRANSLATION /* (*tU[k]).x = ((*qU[k]).x*dt)/((1. - (rho_f/p[k].rho_s))*p[k].M); */ /* (*tU[k]).y = ((*qU[k]).y*dt)/((1. - (rho_f/p[k].rho_s))*p[k].M); */ /* (*tU[k]).z = ((*qU[k]).z*dt)/((1. - (rho_f/p[k].rho_s))*p[k].M); */ /* Modification for the explicit added term */ (*tU[k]).x = ((*qU[k]).x*dt)/( p[k].DLMFD_couplingfactor * p[k].M ); (*tU[k]).y = ((*qU[k]).y*dt)/( p[k].DLMFD_couplingfactor * p[k].M ); (*tU[k]).z = ((*qU[k]).z*dt)/( p[k].DLMFD_couplingfactor * p[k].M ); #endif #if ROTATION /* The inertia tensor is */ /* Ixx -Ixy -Ixz */ /* -Iyx Iyy -Iyz */ /* -Izx -Izy Izz */ /* with */ /* Ip[0] = Ixx */ /* Ip[1] = Iyy */ /* Ip[2] = Izz */ /* Ip[3] = Ixy */ /* Ip[4] = Ixz */ /* Ip[5] = Iyz */ /* Allocate temporly a matrix for Ip */ double ** Imat = malloc(3 * sizeof(double*)); for (int ii = 0; ii < 3; ii++) { Imat[ii] = malloc(3 * sizeof(double)); } Imat[0][0] = p[k].Ip[0]; Imat[1][1] = p[k].Ip[1]; Imat[2][2] = p[k].Ip[2]; Imat[0][1] = -p[k].Ip[3]; Imat[0][2] = -p[k].Ip[4]; Imat[1][2] = -p[k].Ip[5]; Imat[1][0] = Imat[0][1]; Imat[2][0] = Imat[0][2]; Imat[2][1] = Imat[1][2]; double Aw[3][3]; for (int ll = 0; ll < 3; ll++) { for (int nn = 0; nn < 3; nn++) { /* Aw[ll][nn] = (1. - rho_f/p[k].rho_s)*Imat[ll][nn]/dt; */ /* Modification for the explicit added term */ Aw[ll][nn] = p[k].DLMFD_couplingfactor * Imat[ll][nn] / dt; } } /* Takes in Aw and return Aw^-1 as Imat */ inverse3by3matrix(Aw, Imat); (*tw[k]).x = Imat[0][0]*(*qw[k]).x + Imat[0][1]*(*qw[k]).y + Imat[0][2]*(*qw[k]).z; (*tw[k]).y = Imat[1][0]*(*qw[k]).x + Imat[1][1]*(*qw[k]).y + Imat[1][2]*(*qw[k]).z; (*tw[k]).z = Imat[2][0]*(*qw[k]).x + Imat[2][1]*(*qw[k]).y + Imat[2][2]*(*qw[k]).z; for (int ii = 0; ii < 3; ii++) free(Imat[ii]); free(Imat); /* (*tw[k]).x = ((*qw[k]).x*dt)/((1. - (rho_f/p[k].rho_s))*p[k].Ip[0]); */ /* (*tw[k]).y = ((*qw[k]).y*dt)/((1. - (rho_f/p[k].rho_s))*p[k].Ip[0]); */ /* (*tw[k]).z = ((*qw[k]).z*dt)/((1. - (rho_f/p[k].rho_s))*p[k].Ip[0]); */ #endif } #endif /* (3) Compute residual y = M*t, the residual vector y = _P(t) */ /* Sign error in eq (A.9) in J.Eng. Math, 2011 */ /* Written -M*t, should be +M*t */ #if debugInterior == 0 /* Interior points: y = M*t */ /* Interior points: y = M_u*tu + M_U*tU + M_w*tw */ /* So y = _P(t) - _P(t) - _P(t) */ for (int k = 0; k < NPARTICLES; k++) { foreach_cache((*Interior[k])) { if ((flagfield[] < 1) && ((int)index_lambda.y[] == k)) { DLM_v.x[] = tu.x[]; DLM_v.y[] = tu.y[]; DLM_v.z[] = tu.z[]; #if DLM_Moving_particle #if TRANSLATION DLM_v.x[] += -(*tU[k]).x; DLM_v.y[] += -(*tU[k]).y; DLM_v.z[] += -(*tU[k]).z; #endif #if ROTATION /* Modify temporerly the particle center position for periodic boundary condition */ foreach_dimension() (*gci[k]).center.x += DLM_periodic_shift.x[]; DLM_v.x[] += -((*tw[k]).y*(z - (*gci[k]).center.z) - (*tw[k]).z*(y - (*gci[k]).center.y)); // -(t_y*r_z - t_z*r_y) DLM_v.y[] += -((*tw[k]).z*(x - (*gci[k]).center.x) - (*tw[k]).x*(z - (*gci[k]).center.z)); // -(t_z*r_x - t_x*r_z) DLM_v.z[] += -((*tw[k]).x*(y - (*gci[k]).center.y) - (*tw[k]).y*(x - (*gci[k]).center.x)); // -(t_x*r_y - t_y*r_x) foreach_dimension() (*gci[k]).center.x -= DLM_periodic_shift.x[]; #endif #endif } } } #endif #if debugBD == 0 /* Boundary points: y = M*t */ /* Boundary points: y = M_u*tu + M_U*tU + M_w*tw */ /* So y = _P(t) - _P(t) - _P(t) */ boundary ((scalar*){tu}); for (int k = 0; k < NPARTICLES; k++) { particle * pp = &p[k]; foreach_cache((*Boundary[k])) { if (index_lambda.x[] > -1 && ((int)index_lambda.y[] == pp->pnum)) { lambdacellpos.x = x; lambdacellpos.y = y;lambdacellpos.z = 0.; lambdapos.x = (*sbm[k]).x[(int)index_lambda.x[]]; lambdapos.y = (*sbm[k]).y[(int)index_lambda.x[]]; lambdapos.z = 0.; #if dimension == 3 lambdacellpos.z = z; lambdapos.z = (*sbm[k]).z[(int)index_lambda.x[]]; #endif sum.x = 0; sum.y = 0; sum.z = 0; double testweight = 0.; foreach_dimension() ppshift.x = DLM_periodic_shift.x[]; foreach_neighbor() { if (level == depth()) { weightcellpos.x = x; weightcellpos.y = y; weightcellpos.z = 0.; #if dimension == 3 weightcellpos.z = z; #endif weight = reversed_weight (pp, weightcellpos, lambdacellpos, lambdapos, Delta, ppshift, 1., 0.); testweight += weight; sum.x += weight*tu.x[]; sum.y += weight*tu.y[]; sum.z += weight*tu.z[]; } } if ((testweight < 1. - 0.000001) || (testweight > 1. + 0.000001)) printf("testweight = %f\n",testweight); DLM_v.x[] = sum.x; DLM_v.y[] = sum.y; DLM_v.z[] = sum.z; #if DLM_Moving_particle #if TRANSLATION DLM_v.x[] += -(*tU[k]).x; DLM_v.y[] += -(*tU[k]).y; DLM_v.z[] += -(*tU[k]).z; #endif #if ROTATION /* Modify temporerly the particle center position for periodic boundary condition */ foreach_dimension() (*gci[k]).center.x += DLM_periodic_shift.x[]; DLM_v.x[] += -((*tw[k]).y*(lambdapos.z - (*gci[k]).center.z) - (*tw[k]).z*(lambdapos.y - (*gci[k]).center.y)); // -(t_y*r_z - t_z*r_y) DLM_v.y[] += -((*tw[k]).z*(lambdapos.x - (*gci[k]).center.x) - (*tw[k]).x*(lambdapos.z - (*gci[k]).center.z)); // -(t_z*r_x - t_x*r_z) DLM_v.z[] += -((*tw[k]).x*(lambdapos.y - (*gci[k]).center.y) - (*tw[k]).y*(lambdapos.x - (*gci[k]).center.x)); // -(t_x*r_y - t_y*r_x) foreach_dimension() (*gci[k]).center.x -= DLM_periodic_shift.x[]; #endif #endif } } } #endif /* (4) Compute alpha = r^(k-1)*r^(k-1) / DLM_w.y */ DLM_nr2_km1 = DLM_nr2; DLM_wv = 0.; foreach() { DLM_wv += DLM_w.x[]*DLM_v.x[] + DLM_w.y[]*DLM_v.y[] + DLM_w.z[]*DLM_v.z[]; } #if _MPI mpi_all_reduce (DLM_wv, MPI_DOUBLE, MPI_SUM); #endif DLM_alpha = DLM_nr2_km1/DLM_wv; /* (5) Update lambda = lambda - alpha*DLM_w and r = r - alpha*y */ foreach() { DLM_lambda.x[] -= DLM_alpha*DLM_w.x[]; DLM_lambda.y[] -= DLM_alpha*DLM_w.y[]; DLM_lambda.z[] -= DLM_alpha*DLM_w.z[]; DLM_r.x[] -= DLM_alpha*DLM_v.x[]; DLM_r.y[] -= DLM_alpha*DLM_v.y[]; DLM_r.z[] -= DLM_alpha*DLM_v.z[]; /* (7) Update u = u + alpha*t */ u.x[] += DLM_alpha*tu.x[]; u.y[] += DLM_alpha*tu.y[]; u.z[] += DLM_alpha*tu.z[]; } #if DLM_Moving_particle for (int k = 0; k < NPARTICLES; k++) { #if TRANSLATION (*U[k]).x += DLM_alpha*(*tU[k]).x; (*U[k]).y += DLM_alpha*(*tU[k]).y; (*U[k]).z += DLM_alpha*(*tU[k]).z; #endif #if ROTATION (*w[k]).x += DLM_alpha*(*tw[k]).x; (*w[k]).y += DLM_alpha*(*tw[k]).y; (*w[k]).z += DLM_alpha*(*tw[k]).z; #endif } #endif /* (8) Compute beta = nr2^k / nr2^(k-1) */ DLM_nr2 = 0.; /* (6) compute norme ||r^k||^2 = DLM_nr2 */ foreach() { DLM_nr2 += sq(DLM_r.x[]) + sq(DLM_r.y[]) + sq(DLM_r.z[]); } #if _MPI mpi_all_reduce (DLM_nr2, MPI_DOUBLE, MPI_SUM); #endif DLM_beta = DLM_nr2/DLM_nr2_km1; /* (9) Update DLM_w = r + beta*DLM_w */ /* Sign error in eq (A.15) of J.Eng. Math, 2011 */ /* Written -beta*DLM_w, should be +beta*DLM_w */ foreach() { DLM_w.x[] = DLM_r.x[] + DLM_beta*DLM_w.x[]; DLM_w.y[] = DLM_r.y[] + DLM_beta*DLM_w.y[]; DLM_w.z[] = DLM_r.z[] + DLM_beta*DLM_w.z[]; } /* // Output nr */ /* fprintf (stderr,"Iteration %d ||u-u_imposed|| = %12.10e\n",ki,sqrt(DLM_nr2)); */ } if (i == 0) fprintf (converge, "# Time Iteration \t Uzawa Iterations \t ||u-u_imposed||\n"); fprintf (converge,"%d \t %d \t %12.10e\n", i, ki, sqrt(DLM_nr2)); fflush(converge); #if DLM_Moving_particle #if TRANSLATION fprintf (stderr,"DLMFD: particle-0's velocity on thread %d is (%20.18f, %20.18f, %20.18f)\n",pid(), (*U[0]).x, (*U[0]).y, (*U[0]).z); coord cc = (*gci[0]).center; fprintf (stderr,"DLMFD: particle-0's position on thread %d is (%20.18f, %20.18f, %20.18f)\n",pid(), cc.x, cc.y, cc.z); #endif #if ROTATION fprintf (stderr,"DLMFD: particle-0's rotation on thread %d is (%20.18f, %20.18f, %20.18f)\n",pid(), (*w[0]).x, (*w[0]).y, (*w[0]).z); #endif #endif /* Compute the explicit term here */ #if DLM_alpha_coupling boundary ((scalar*) {DLM_lambda}); coord nn = {0., 0., 0.}; foreach_dimension() { nn.x = 0.; norm n = normf (DLM_lambda.x); nn.x = n.max; } fprintf (stderr, "# Lambda explicit term: max of each component (%g, %g, %g)\n", nn.x, nn.y, nn.z); // log for (int k = 0; k < NPARTICLES; k++) { # if debugBD == 0 particle * pp = &p[k]; foreach_cache ((*Boundary[k])) { sum.x = 0.; sum.y = 0.; sum.z = 0.; weightcellpos.x = x; weightcellpos.y = y; #if dimension == 3 weightcellpos.z = z; #endif foreach_neighbor() { if (((int)index_lambda.x[] > -1) && (level == depth()) && is_leaf(cell) && ((int)index_lambda.y[]) == k) { lambdacellpos.x = x; lambdacellpos.y = y; lambdapos.x = (*sbm[k]).x[(int)index_lambda.x[]]; lambdapos.y = (*sbm[k]).y[(int)index_lambda.x[]]; #if dimension == 3 lambdacellpos.z = z; lambdapos.z = (*sbm[k]).z[(int)index_lambda.x[]]; #endif foreach_dimension() ppshift.x = DLM_periodic_shift.x[]; weight = reversed_weight (pp, weightcellpos, lambdacellpos, lambdapos, Delta, ppshift,1.,1.); foreach_dimension() sum.x += weight*DLM_lambda.x[]; } } foreach_dimension() { DLM_explicit.x[] += sum.x; } } #endif #if debugInterior == 0 #if POISEUILLE /* For the poiseuille flow, the explicit forcing term is too sharp and needs to be "smoothed" */ #define smooth 0.5 #else #define smooth 1. #endif foreach_cache ((*Interior[k])) { if ((flagfield[] < 1) && ((int)index_lambda.y[] == k)) foreach_dimension() { DLM_explicit.x[] = smooth*DLM_lambda.x[]; } } #endif } boundary((scalar *) {DLM_explicit}); #endif /* Added by Arthur Ghigo, October 10, 2019 */ boundary ((scalar *) {u}); /* Timers and Timings */ /* give 1 as number of cell so that timer_timing does not compute the total number of cells*/ /* the total number of cell is used to compute the speed of the solver which does not make sens for a single step here, we will compute it globally at the end of the run*/ timing dlmfd_timing = timer_timing (dlmfd_timer, 1, 1, mpitimings); dlmfd_globaltiming.cpu += dlmfd_timing.cpu; dlmfd_globaltiming.real += dlmfd_timing.real; dlmfd_globaltiming.min += dlmfd_timing.min; dlmfd_globaltiming.avg += dlmfd_timing.avg; dlmfd_globaltiming.max += dlmfd_timing.max; }