sandbox/cselcuk/DLMFD_reverse_Uzawa.h
Fictitious-domain method with distributed Lagrangian multipliers
We are looking for a numerical approximation for the solution of the following equations \displaystyle \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 \displaystyle \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) \displaystyle \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) \displaystyle \mathbf{u}-\left(\mathbf{U}+\mathbf{\omega}\times \mathbf{r}\right)=0~\text{over}~P(t) \displaystyle \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. \displaystyle \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 \displaystyle {\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: \displaystyle \rho \partial_t\mathbf{u} = {\mathbf{\lambda}}~\text{over}~\Omega \displaystyle \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) \displaystyle \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) \displaystyle \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 <=> <lambda,v>_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 <=> - <lambda, V>_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 <=> - <lambda, xi^r_GM>_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 -<lambda,v>_P(t) */
/* Interior points qU = fU -(M_U^T)*lambda^0 */
/* qU = fU -(M_U^T)*lambda^0 = fU + <lambda,V>_P(t) */
/* Interior points qw = fw -(M_w^T)*lambda^0 */
/* qw = fw -(M_w^T)*lambda^0 = fw + <lambda,xi^r_GM>_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) */
/* <lambda,xi^r_GM>_P(t) = <r_GM, lambda^xi>_P(t) = <xi,r_GM^lambda>_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) */
/* <lambda,xi^r_GM>_P(t) = <r_GM, lambda^xi>_P(t) = <xi,r_GM^lambda>_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 *555/span>
Identity_Matrix */
/* Invert A_w w^0 = qw with A_w = ((1 - rho_f/rho_s)*Ip)/dt *558/span>
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 <=> <alpha, u>_P(t) */
/* Note that M_U*U^0 <=> - <alpha, U>_P(t) */
/* Note that M_w*w^0 <=> - <alpha, w^r_GM>_P(t) */
/* alpha is he test function for lambda the Lagrange multipliers */
/* So r^0 = G -<alpha, u>_P(t) + <alpha, U>_P(t) + <alpha, w^r_GM>_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 -<alpha, u>_P(t) + <alpha, U>_P(t) + <alpha, w^r_GM>_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 -<alpha, u>_P(t) + <alpha, U>_P(t) + <alpha, w^r_GM>_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 = <DLM_w, v>_P(t) */
/* Interior points qU = (M_U^T)DLM_w = - <DLM_w, V>_P(t) */
/* Interior points qw = (M_w^T)DLM_w = - <DLM_w, xi^r_GM>_P(t) */
/* - <DLM_w, xi^r_GM>_P(t) = - <r_GM, DLM_w^xi>_P(t) = - <xi, r_GM^DLM_w>_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 = <DLM_w, v>_P(t) */
/* Boundary points qU = (M_U^T)DLM_w = - <DLM_w, V>_P(t) */
/* Boundary points qw = (M_w^T)DLM_w = - <DLM_w, xi^r_GM>_P(t) */
/* - <DLM_w, xi^r_GM>_P(t) = - <r_GM, DLM_w^xi>_P(t) = - <xi, r_GM^DLM_w>_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 = <alpha, tu-(tU+tw^GM)>_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 = <alpha, tu>_P(t) - <alpha, tU>_P(t) - <alpha, tw^r_GM>_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 = <alpha, tu>_P(t) - <alpha, tU>_P(t) - <alpha, tw^r_GM>_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;
}