1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14
15
16 /* ----------------------------------------------------------------------
17 Contributing authors: Frances Mackay, Santtu Ollila, Colin Denniston (UWO)
18 ------------------------------------------------------------------------- */
19
20 #include "fix_lb_fluid.h"
21
22 #include "atom.h"
23 #include "comm.h"
24 #include "domain.h"
25 #include "error.h"
26 #include "force.h"
27 #include "group.h"
28 #include "memory.h"
29 #include "modify.h"
30 #include "random_mars.h"
31 #include "update.h"
32
33 #include <cmath>
34 #include <cstring>
35 #include <algorithm>
36 #include <utility>
37
38 using namespace LAMMPS_NS;
39 using namespace FixConst;
40
41 static const double kappa_lb=0.0;
42
FixLbFluid(LAMMPS * lmp,int narg,char ** arg)43 FixLbFluid::FixLbFluid(LAMMPS *lmp, int narg, char **arg) :
44 Fix(lmp, narg, arg)
45 {
46 //=====================================================================================================
47 // Sample inputfile call:
48 // fix # group lb/fluid nevery typeLB viscosity densityinit_real
49 //
50 // where: nevery: call this fix every nevery timesteps.
51 // (keep this set to 1 for now).
52 // typeLB: there are two different integrators
53 // in the code labeled "1" and "2".
54 // viscosity: the viscosity of the fluid.
55 // densityinit_real: the density of the fluid.
56 //
57 // optional arguments:
58 // "setArea" type node_area: set the surface area per node associated with a
59 // given atom type. By default the surface area
60 // is set at 1.0*dx_lb^2.
61 // "setGamma" gamma: specify a user-defined value for the force
62 // coupling constant, instead of using the default
63 // value.
64 // "scaleGamma" type scale_factor: scale the user provided force coupling constant
65 // by the factor, scale_factor, for the given atom
66 // type.
67 // "dx" dx_lb: the lattice-Boltzmann grid spacing.
68 // "dm" dm_lb: the lattice-Boltzmann mass unit.
69 // "a0" a_0_real: the square of the sound speed in the fluid.
70 // "noise" Temperature seed: include noise in the system.
71 // Temperature is the temperature for the fluid.
72 // seed is the seed for the random number generator.
73 // "calcforce" N group: print the force acting on a given group every
74 // N timesteps.
75 // "trilinear": use the trilinear interpolation stencil.
76 // "read_restart" restart_file: restart a fluid run from restart_file.
77 // "write_restart" N: write a fluid restart file every N timesteps.
78 // "zwall_velocity" velocity_bottom velocity_top: assign velocities to the z-walls
79 // in the system.
80 // "bodyforce" bodyforcex bodyforcey bodyforcez: add a constant body force to the
81 // fluid.
82 // "printfluid" N: print the fluid density and velocity at each
83 // grid point every N timesteps.
84 // "D3Q19": use the 19 velocity D3Q19 model. By default,
85 // the 15 velocity D3Q15 model is used.
86 //=====================================================================================================
87
88 if (narg <7) error->all(FLERR,"Illegal fix lb/fluid command");
89
90 if (comm->style != 0)
91 error->universe_all(FLERR,"Fix lb/fluid can only currently be used with "
92 "comm_style brick");
93
94 MPI_Comm_rank(world,&me);
95 MPI_Comm_size(world,&nprocs);
96
97 nevery = atoi(arg[3]);
98 typeLB = atoi(arg[4]);
99 viscosity = atof(arg[5]);
100 densityinit_real = atof(arg[6]);
101
102 // Default values for optional arguments:
103 force_diagnostic=0;
104 noisestress = 0;
105 trilinear_stencil = 0;
106 readrestart = 0;
107 printrestart = 0;
108 bodyforcex = bodyforcey = bodyforcez = 0.0;
109 vwtp = vwbt = 0.0;
110 printfluid = 0;
111 T = 300.0;
112 dm_lb = 1.0;
113 fixviscouslb = 0;
114 setdx = 1;
115 seta0 = 1;
116 setGamma = 0;
117 setArea = 0;
118 numvel = 15;
119
120 Gamma = nullptr;
121 NodeArea = nullptr;
122
123 int iarg = 7;
124 while (iarg < narg) {
125 if (strcmp(arg[iarg],"setArea")==0) {
126 if (setGamma == 1)
127 error->all(FLERR,"Illegal fix lb/fluid command: cannot use a combination of default and user-specified gamma values");
128 setArea = 1;
129 int itype = atoi(arg[iarg+1]);
130 double areafactor = atof(arg[iarg+2]);
131 if (itype <= 0 || itype > atom->ntypes || areafactor < 0.0)
132 error->all(FLERR,"Illegal fix lb/fluid command: setArea");
133 if (NodeArea == nullptr) {
134 NodeArea = new double[atom->ntypes+1];
135 for (int i=0; i<=atom->ntypes; i++) NodeArea[i] = -1.0;
136 }
137 NodeArea[itype] = areafactor;
138 iarg += 3;
139 }
140 else if (strcmp(arg[iarg],"setGamma")==0) {
141 if (setArea == 1)
142 error->all(FLERR,"Illegal fix lb/fluid command: cannot use a combination of default and user-specified gamma values");
143 setGamma = 1;
144 double Gammaone;
145 Gammaone = atof(arg[iarg+1]);
146 if (Gamma == nullptr)
147 Gamma = new double[atom->ntypes+1];
148 for (int i=0; i<=atom->ntypes; i++) Gamma[i] = Gammaone;
149 iarg += 2;
150 }
151 else if (strcmp(arg[iarg],"scaleGamma")==0) {
152 if (setGamma == 0)
153 error->all(FLERR,"Illegal fix lb/fluid command: must set a value for Gamma before scaling it");
154 int itype = atoi(arg[iarg+1]);
155 double scalefactor = atof(arg[iarg+2]);
156 if (itype <= 0 || itype > atom->ntypes || scalefactor < 0.0)
157 error->all(FLERR,"Illegal fix lb/fluid command: scaleGamma");
158 Gamma[itype] *= scalefactor;
159 iarg += 3;
160 }
161 else if (strcmp(arg[iarg],"dx")==0) {
162 dx_lb = atof(arg[iarg+1]);
163 iarg += 2;
164 setdx = 0;
165 }
166 else if (strcmp(arg[iarg],"dm")==0) {
167 dm_lb = atof(arg[iarg+1]);
168 iarg += 2;
169 }
170 else if (strcmp(arg[iarg],"a0")==0) {
171 a_0_real = atof(arg[iarg+1]);
172 iarg += 2;
173 seta0 = 0;
174 }
175 else if (strcmp(arg[iarg],"noise")== 0) {
176 noisestress = 1;
177 T = atof(arg[iarg+1]);
178 seed = atoi(arg[iarg+2]);
179 iarg += 3;
180 }
181 else if (strcmp(arg[iarg],"calcforce")==0) {
182 force_diagnostic = atoi(arg[iarg+1]);
183 igroupforce=group->find(arg[iarg+2]);
184 iarg += 3;
185 }
186 else if (strcmp(arg[iarg],"trilinear")==0) {
187 trilinear_stencil = 1;
188 iarg += 1;
189 }
190 else if (strcmp(arg[iarg],"read_restart")==0) {
191 readrestart = 1;
192 int nlength = strlen(arg[iarg+1]) + 16;
193 char *filename = new char[nlength];
194 strcpy(filename,arg[iarg+1]);
195 MPI_File_open(world,filename,MPI_MODE_RDONLY,MPI_INFO_NULL,&pFileRead);
196 delete [] filename;
197 iarg += 2;
198 }
199 else if (strcmp(arg[iarg],"write_restart")==0) {
200 printrestart = atoi(arg[iarg+1]);
201 iarg += 2;
202 }
203 else if (strcmp(arg[iarg],"zwall_velocity")==0) {
204 if (domain->periodicity[2]!=0) error->all(FLERR,"fix lb/fluid error: setting \
205 a z wall velocity without implementing fixed BCs in z");
206 vwbt = atof(arg[iarg+1]);
207 vwtp = atof(arg[iarg+2]);
208 iarg += 3;
209 }
210 else if (strcmp(arg[iarg],"bodyforce")==0) {
211 bodyforcex = atof(arg[iarg+1]);
212 bodyforcey = atof(arg[iarg+2]);
213 bodyforcez = atof(arg[iarg+3]);
214 iarg += 4;
215 }
216 else if (strcmp(arg[iarg],"printfluid")==0) {
217 printfluid = atoi(arg[iarg+1]);
218 iarg += 2;
219 }
220 else if (strcmp(arg[iarg],"D3Q19")==0) {
221 numvel = 19;
222 iarg += 1;
223 }
224 else error->all(FLERR,"Illegal fix lb/fluid command");
225 }
226
227 //--------------------------------------------------------------------------
228 //Choose between D3Q15 and D3Q19 functions:
229 //--------------------------------------------------------------------------
230 if (numvel == 15) {
231 initializeLB = &FixLbFluid::initializeLB15;
232 equilibriumdist = &FixLbFluid::equilibriumdist15;
233 update_full = &FixLbFluid::update_full15;
234 } else {
235 initializeLB = &FixLbFluid::initializeLB19;
236 equilibriumdist = &FixLbFluid::equilibriumdist19;
237 update_full = &FixLbFluid::update_full19;
238 }
239
240 //--------------------------------------------------------------------------
241 // perform initial allocation of atom-based array register
242 // with Atom class
243 //--------------------------------------------------------------------------
244 hydroF = nullptr;
245 grow_arrays(atom->nmax);
246 atom->add_callback(Atom::GROW);
247
248 for (int i=0; i<atom->nmax; i++)
249 for (int j=0; j<3; j++)
250 hydroF[i][j] = 0.0;
251
252 Ng_lb = nullptr;
253 w_lb = nullptr;
254 mg_lb = nullptr;
255 e = nullptr;
256 feq = nullptr;
257 feqold = nullptr;
258 feqn = nullptr;
259 feqoldn = nullptr;
260 f_lb = nullptr;
261 fnew = nullptr;
262 density_lb = nullptr;
263 u_lb = nullptr;
264 altogether = nullptr;
265 buf = nullptr;
266 Ff = nullptr;
267 Fftempx = nullptr;
268 Fftempy = nullptr;
269 Fftempz = nullptr;
270
271 //--------------------------------------------------------------------------
272 // Set the lattice Boltzmann dt.
273 //--------------------------------------------------------------------------
274 dt_lb=nevery*(update->dt);
275
276 //--------------------------------------------------------------------------
277 // Set the lattice Boltzmann dx if it wasn't specified in the
278 // input.
279 //--------------------------------------------------------------------------
280 if (setdx == 1) {
281 double dx_lb1 = sqrt(3.0*viscosity*dt_lb/densityinit_real);
282 double mindomain = std::min(std::min(domain->xprd/comm->procgrid[0],domain->yprd/comm->procgrid[1]),domain->zprd/comm->procgrid[2]);
283 dx_lb = mindomain/floor(mindomain/dx_lb1);
284
285 if (comm->me==0) {
286 char str[128];
287 sprintf(str,"Setting the lattice-Boltzmann dx to %10.6f",dx_lb);
288 error->message(FLERR,str);
289 }
290 }
291 //--------------------------------------------------------------------------
292 // If the area per node has not been set by the user, set to the
293 // default value of dx_lb*dx_lb.
294 //--------------------------------------------------------------------------
295 if (setGamma == 0) {
296 if (setArea == 0) {
297 if (comm->me==0) {
298 error->message(FLERR,"Assuming an area per node of dx*dx for all of the MD particles. This should only be used if these all correspond to point particles; otherwise, change using the setArea keyword");
299 }
300 NodeArea = new double[atom->ntypes+1];
301 for (int i=0; i<=atom->ntypes; i++) NodeArea[i] = -1.0;
302 }
303 for (int i=0; i<=atom->ntypes; i++)
304 if (NodeArea[i] < 0.0) NodeArea[i] = dx_lb*dx_lb;
305 }
306 //--------------------------------------------------------------------------
307 // Set a0 if it wasn't specified in the input
308 //--------------------------------------------------------------------------
309 if (seta0 == 1)
310 a_0_real = 0.33333333*dx_lb*dx_lb/dt_lb/dt_lb;
311
312 //--------------------------------------------------------------------------
313 // Check to make sure that the total number of grid points in each direction
314 // divides evenly among the processors in that direction.
315 // Shrink-wrapped boundary conditions (which are not permitted by this fix)
316 // might cause a problem, so check for this. A full check of the boundary
317 // conditions is performed in the init routine, rather than here, as it is
318 // possible to change the BCs between runs.
319 //--------------------------------------------------------------------------
320 double aa;
321 double eps=1.0e-8;
322 aa = (domain->xprd/comm->procgrid[0])/dx_lb;
323 if (fabs(aa - floor(aa+0.5)) > eps) {
324 if (domain->boundary[0][0] != 0) {
325 error->all(FLERR,"the x-direction must be periodic");
326 }
327 char errormessage[200];
328 sprintf(errormessage,"With dx= %f, and the simulation domain divided by %i processors in the x direction, the simulation domain in the x direction must be a multiple of %f",dx_lb,comm->procgrid[0],comm->procgrid[0]*dx_lb);
329 error->all(FLERR,errormessage);
330 }
331 aa = (domain->yprd/comm->procgrid[1])/dx_lb;
332 if (fabs(aa - floor(aa+0.5)) > eps) {
333 if (domain->boundary[1][0] != 0) {
334 error->all(FLERR,"the y-direction must be periodic");
335 }
336 char errormessage[200];
337 sprintf(errormessage,"With dx= %f, and the simulation domain divided by %i processors in the y direction, the simulation domain in the y direction must be a multiple of %f",dx_lb,comm->procgrid[1],comm->procgrid[1]*dx_lb);
338 error->all(FLERR,errormessage);
339 }
340 aa = (domain->zprd/comm->procgrid[2])/dx_lb;
341 if (fabs(aa - floor(aa+0.5)) > eps) {
342 if (domain->boundary[2][0] == 2 || domain->boundary[2][0] == 3) {
343 error->all(FLERR,"the z-direction can not have shrink-wrap boundary conditions");
344 }
345 char errormessage[200];
346 sprintf(errormessage,"With dx= %f, and the simulation domain divided by %i processors in the z direction, the simulation domain in the z direction must be a multiple of %f",dx_lb,comm->procgrid[2],comm->procgrid[2]*dx_lb);
347 error->all(FLERR,errormessage);
348 }
349
350 //--------------------------------------------------------------------------
351 // Set the total number of grid points in each direction.
352 //--------------------------------------------------------------------------
353 Nbx = (int)(domain->xprd/dx_lb + 0.5);
354 Nby = (int)(domain->yprd/dx_lb + 0.5);
355 Nbz = (int)(domain->zprd/dx_lb + 0.5);
356
357 //--------------------------------------------------------------------------
358 // Set the number of grid points in each dimension for the local subgrids.
359 //--------------------------------------------------------------------------
360 subNbx= Nbx/comm->procgrid[0] + 2;
361 subNby= Nby/comm->procgrid[1] + 2;
362 subNbz= Nbz/comm->procgrid[2] + 2;
363
364 //--------------------------------------------------------------------------
365 // In order to calculate the fluid forces correctly, need to have atleast
366 // 5 grid points in each direction per processor.
367 //--------------------------------------------------------------------------
368 if (subNbx<7 || subNby < 7 || subNbz<7)
369 error->all(FLERR,"Need at least 5 grid points in each direction per processor");
370
371 // If there are walls in the z-direction add an extra grid point.
372 if (domain->periodicity[2]==0) {
373 Nbz += 1;
374 if (comm->myloc[2]==comm->procgrid[2]-1)
375 subNbz += 1;
376 }
377
378 if (comm->me==0) {
379 char str[128];
380 if (setdx == 1) {
381 sprintf(str,"Using a lattice-Boltzmann grid of %i by %i by %i total grid points. To change, use the dx keyword",Nbx,Nby,Nbz);
382 } else {
383 sprintf(str,"Using a lattice-Boltzmann grid of %i by %i by %i total grid points.",Nbx,Nby,Nbz);
384 }
385 error->message(FLERR,str);
386 }
387
388 //--------------------------------------------------------------------------
389 // Store the largest value of subNbz, which is needed for allocating the
390 // buf array (since a processor with comm->myloc[2] == comm->procgrid[2]-1
391 // may have an additional subNbz point as compared with the rest).
392 //--------------------------------------------------------------------------
393 int subNbzmax;
394 MPI_Allreduce(&subNbz,&subNbzmax,1,MPI_INT,MPI_MAX,world);
395
396 //--------------------------------------------------------------------------
397 // Create the MPI datatypes used to pass portions of arrays:
398 // datatypes to pass the f and feq arrays.
399 //--------------------------------------------------------------------------
400 MPI_Aint lb,sizeofdouble;
401 MPI_Type_get_extent(MPI_DOUBLE,&lb,&sizeofdouble);
402
403 MPI_Type_vector(subNbz-2,numvel,numvel,MPI_DOUBLE,&oneslice);
404 MPI_Type_commit(&oneslice);
405 MPI_Type_create_hvector(subNby-2,1,numvel*subNbz*sizeofdouble,oneslice,&passxf);
406 MPI_Type_commit(&passxf);
407
408 MPI_Type_create_hvector(subNbx,1,numvel*subNbz*subNby*sizeofdouble,oneslice,&passyf);
409 MPI_Type_commit(&passyf);
410
411 MPI_Type_free(&oneslice);
412 MPI_Type_vector(subNby,numvel,numvel*subNbz,MPI_DOUBLE,&oneslice);
413 MPI_Type_commit(&oneslice);
414 MPI_Type_create_hvector(subNbx,1,numvel*subNbz*subNby*sizeofdouble,oneslice,&passzf);
415 MPI_Type_commit(&passzf);
416
417 // datatypes to pass the u array, and the Ff array.
418 MPI_Type_free(&oneslice);
419 MPI_Type_vector(subNbz+3,3,3,MPI_DOUBLE,&oneslice);
420 MPI_Type_commit(&oneslice);
421 MPI_Type_create_hvector(subNby+3,1,3*(subNbz+3)*sizeofdouble,oneslice,&passxu);
422 MPI_Type_commit(&passxu);
423
424 MPI_Type_create_hvector(subNbx+3,1,3*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passyu);
425 MPI_Type_commit(&passyu);
426
427 MPI_Type_free(&oneslice);
428 MPI_Type_vector(subNby+3,3,3*(subNbz+3),MPI_DOUBLE,&oneslice);
429 MPI_Type_commit(&oneslice);
430 MPI_Type_create_hvector(subNbx+3,1,3*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passzu);
431 MPI_Type_commit(&passzu);
432
433 // datatypes to pass the density array.
434 MPI_Type_free(&oneslice);
435 MPI_Type_vector(subNbz+3,1,1,MPI_DOUBLE,&oneslice);
436 MPI_Type_commit(&oneslice);
437 MPI_Type_create_hvector(subNby+3,1,1*(subNbz+3)*sizeofdouble,oneslice,&passxrho);
438 MPI_Type_commit(&passxrho);
439
440 MPI_Type_create_hvector(subNbx+3,1,1*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passyrho);
441 MPI_Type_commit(&passyrho);
442
443 MPI_Type_free(&oneslice);
444 MPI_Type_vector(subNby+3,1,1*(subNbz+3),MPI_DOUBLE,&oneslice);
445 MPI_Type_commit(&oneslice);
446 MPI_Type_create_hvector(subNbx+3,1,1*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passzrho);
447 MPI_Type_commit(&passzrho);
448
449 // datatypes to receive a portion of the Ff array.
450 MPI_Type_free(&oneslice);
451 MPI_Type_vector(subNbz+3,3,3,MPI_DOUBLE,&oneslice);
452 MPI_Type_commit(&oneslice);
453 MPI_Type_create_hvector(subNby+3,1,3*(subNbz+3)*sizeofdouble,oneslice,&passxtemp);
454 MPI_Type_commit(&passxtemp);
455
456 MPI_Type_create_hvector(subNbx+3,1,3*(subNbz+3)*5*sizeofdouble,oneslice,&passytemp);
457 MPI_Type_commit(&passytemp);
458
459 MPI_Type_free(&oneslice);
460 MPI_Type_vector(subNby+3,3,3*5,MPI_DOUBLE,&oneslice);
461 MPI_Type_commit(&oneslice);
462 MPI_Type_create_hvector(subNbx+3,1,3*5*(subNby+3)*sizeofdouble,oneslice,&passztemp);
463 MPI_Type_commit(&passztemp);
464
465 MPI_Type_free(&oneslice);
466
467 //--------------------------------------------------------------------------
468 // Allocate the necessary arrays.
469 //--------------------------------------------------------------------------
470 memory->create(Ng_lb,numvel,"FixLbFluid:Ng_lb");
471 memory->create(w_lb,numvel,"FixLbFluid:w_lb");
472 memory->create(mg_lb,numvel,numvel,"FixLbFluid:mg_lb");
473 memory->create(e,numvel,3,"FixLbFluid:e");
474 memory->create(feq,subNbx,subNby,subNbz,numvel,"FixLbFluid:feq");
475 if (typeLB == 2) {
476 memory->create(feqold,subNbx,subNby,subNbz,numvel,"FixLbFluid:feqold");
477 memory->create(feqn,subNbx,subNby,subNbz,numvel,"FixLbFluid:feqn");
478 memory->create(feqoldn,subNbx,subNby,subNbz,numvel,"FixLbFluid:feqoldn");
479 }
480 memory->create(f_lb,subNbx,subNby,subNbz,numvel,"FixLbFluid:f_lb");
481 memory->create(fnew,subNbx,subNby,subNbz,numvel,"FixLbFluid:fnew");
482 memory->create(density_lb,subNbx+3,subNby+3,subNbz+3,"FixLbFluid:density_lb");
483 memory->create(u_lb,subNbx+3,subNby+3,subNbz+3,3,"FixLbFluid:u_lb");
484 if (printfluid > 0) {
485 memory->create(buf,subNbx,subNby,subNbzmax,4,"FixLbFluid:buf");
486 if (me==0)
487 memory->create(altogether,Nbx,Nby,Nbz,4,"FixLbFluid:altogether");
488 }
489 memory->create(Ff,subNbx+3,subNby+3,subNbz+3,3,"FixLbFluid:Ff");
490 memory->create(Fftempx,5,subNby+3,subNbz+3,3,"FixLbFluid:Fftempx");
491 memory->create(Fftempy,subNbx+3,5,subNbz+3,3,"FixLbFluid:Fftempy");
492 memory->create(Fftempz,subNbx+3,subNby+3,5,3,"FixLbFluid:Fftempz");
493
494 if (noisestress==1) {
495 random = new RanMars(lmp,seed + comm->me);
496 }
497
498 //--------------------------------------------------------------------------
499 // Rescale the variables to Lattice Boltzmann dimensionless units.
500 //--------------------------------------------------------------------------
501 rescale();
502
503 //--------------------------------------------------------------------------
504 // Initialize the arrays.
505 //--------------------------------------------------------------------------
506 (*this.*initializeLB)();
507 initialize_feq();
508
509 }
510
~FixLbFluid()511 FixLbFluid::~FixLbFluid()
512 {
513
514 atom->delete_callback(id,Atom::GROW);
515 memory->destroy(hydroF);
516
517 memory->destroy(Ng_lb);
518 memory->destroy(w_lb);
519 memory->destroy(mg_lb);
520 memory->destroy(e);
521 memory->destroy(feq);
522 if (typeLB == 2) {
523 memory->destroy(feqold);
524 memory->destroy(feqn);
525 memory->destroy(feqoldn);
526 }
527 memory->destroy(f_lb);
528 memory->destroy(fnew);
529 memory->destroy(density_lb);
530 memory->destroy(u_lb);
531 if (printfluid>0) {
532 if (me==0)
533 memory->destroy(altogether);
534 memory->destroy(buf);
535 }
536 memory->destroy(Ff);
537 memory->destroy(Fftempx);
538 memory->destroy(Fftempy);
539 memory->destroy(Fftempz);
540
541 if (noisestress==1) {
542 delete random;
543 }
544
545 if (setGamma == 1) {
546 delete [] Gamma;
547 } else {
548 delete [] NodeArea;
549 }
550 MPI_Type_free(&passxf);
551 MPI_Type_free(&passyf);
552 MPI_Type_free(&passzf);
553 MPI_Type_free(&passxu);
554 MPI_Type_free(&passyu);
555 MPI_Type_free(&passzu);
556 MPI_Type_free(&passxrho);
557 MPI_Type_free(&passyrho);
558 MPI_Type_free(&passzrho);
559 MPI_Type_free(&passxtemp);
560 MPI_Type_free(&passytemp);
561 MPI_Type_free(&passztemp);
562 }
563
setmask()564 int FixLbFluid::setmask()
565 {
566 int mask =0;
567 mask |= INITIAL_INTEGRATE;
568 mask |= POST_FORCE;
569 mask |= END_OF_STEP;
570 return mask;
571 }
572
init(void)573 void FixLbFluid::init(void)
574 {
575 int i,j;
576
577 if (comm->style != 0)
578 error->universe_all(FLERR,"Fix lb/fluid can only currently be used with "
579 "comm_style brick");
580
581 //--------------------------------------------------------------------------
582 // Check to see if the MD timestep has changed between runs.
583 //--------------------------------------------------------------------------
584 double dt_lb_now;
585 dt_lb_now=nevery*(update->dt);
586
587 if (fabs(dt_lb_now - dt_lb) > 1.0e-12) {
588 error->warning(FLERR,"Timestep has changed between runs with the same lb/fluid. Unphysical results may occur");
589 }
590
591 //--------------------------------------------------------------------------
592 // Make sure the size of the simulation domain has not changed
593 // between runs.
594 //--------------------------------------------------------------------------
595 int Nbx_now,Nby_now,Nbz_now;
596 Nbx_now = (int)(domain->xprd/dx_lb + 0.5);
597 Nby_now = (int)(domain->yprd/dx_lb + 0.5);
598 Nbz_now = (int)(domain->zprd/dx_lb + 0.5);
599 // If there are walls in the z-direction add an extra grid point.
600 if (domain->periodicity[2]==0) {
601 Nbz_now += 1;
602 }
603
604 if (Nbx_now != Nbx || Nby_now != Nby || Nbz_now != Nbz) {
605 error->all(FLERR,"the simulation domain can not change shape between runs with the same lb/fluid");
606 }
607
608 //--------------------------------------------------------------------------
609 // Check to make sure that the chosen LAMMPS boundary types are compatible
610 // with this fix.
611 // shrink-wrap is not compatible in any dimension.
612 // fixed only works in the z-direction.
613 //--------------------------------------------------------------------------
614 if (domain->boundary[0][0] != 0) {
615 error->all(FLERR,"the x-direction must be periodic");
616 }
617 if (domain->boundary[1][0] != 0) {
618 error->all(FLERR,"the y-direction must be periodic");
619 }
620 if (domain->boundary[2][0] == 2 || domain->boundary[2][0] == 3) {
621 error->all(FLERR,"the z-direction can not have shrink-wrap boundary conditions");
622 }
623
624 //--------------------------------------------------------------------------
625 // Check if the lb/viscous fix is also called:
626 //--------------------------------------------------------------------------
627 groupbit_viscouslb = groupbit_pc = groupbit_rigid_pc_sphere = 0;
628 for (i = 0; i < modify->nfix; i++) {
629 if (strcmp(modify->fix[i]->style,"lb/viscous") == 0) {
630 fixviscouslb = 1;
631 groupbit_viscouslb = group->bitmask[modify->fix[i]->igroup];
632 }
633 if (strcmp(modify->fix[i]->style,"lb/pc")==0) {
634 groupbit_pc = group->bitmask[modify->fix[i]->igroup];
635 }
636 if (strcmp(modify->fix[i]->style,"lb/rigid/pc/sphere")==0) {
637 groupbit_rigid_pc_sphere = group->bitmask[modify->fix[i]->igroup];
638 }
639 }
640
641 // Warn if the fluid force is not applied to any of the particles.
642 if (!(groupbit_viscouslb || groupbit_pc || groupbit_rigid_pc_sphere) && comm->me==0) {
643 error->message(FLERR,"Not adding the fluid force to any of the MD particles. To add this force use one of the lb/viscous, lb/pc, or lb/rigid/pc/sphere fixes");
644 }
645
646 // If fix lb/viscous is called for a particular atom, make sure
647 // lb/pc or lb/rigid/pc/sphere are not:
648 if (fixviscouslb == 1) {
649 int *mask = atom->mask;
650 int nlocal = atom->nlocal;
651 for (j=0; j<nlocal; j++) {
652 if ((mask[j] & groupbit) && (mask[j] & groupbit_viscouslb) && (mask[j] & groupbit_pc))
653 error->one(FLERR,"should not use the lb/viscous command when integrating with the lb/pc fix");
654 if ((mask[j] & groupbit) && (mask[j] & groupbit_viscouslb) && (mask[j] & groupbit_rigid_pc_sphere))
655 error->one(FLERR,"should not use the lb/viscous command when integrating with the lb/rigid/pc/sphere fix");
656 }
657 }
658
659 }
660
setup(int)661 void FixLbFluid::setup(int /*vflag*/)
662 {
663 //--------------------------------------------------------------------------
664 // Need to calculate the force on the fluid for a restart run.
665 //--------------------------------------------------------------------------
666 if (step > 0)
667 calc_fluidforce();
668 }
669
initial_integrate(int)670 void FixLbFluid::initial_integrate(int /*vflag*/)
671 {
672 //--------------------------------------------------------------------------
673 // Print a header labelling any output printed to the screen.
674 //--------------------------------------------------------------------------
675 static int printheader = 1;
676
677 if (printheader == 1) {
678 if (force_diagnostic > 0 && me == 0) {
679 printf("-------------------------------------------------------------------------------\n");
680 printf(" F_x F_y F_z T_x T_y T_z\n");
681 printf("-------------------------------------------------------------------------------\n");
682 }
683
684 if (printfluid > 0 && me == 0) {
685 printf("---------------------------------------------------------------------\n");
686 printf(" density u_x u_y u_z \n");
687 printf("---------------------------------------------------------------------\n");
688 }
689 printheader = 0;
690 }
691
692 //--------------------------------------------------------------------------
693 // Determine the equilibrium distribution on the local subgrid.
694 //--------------------------------------------------------------------------
695 (*this.*equilibriumdist)(1,subNbx-1,1,subNby-1,1,subNbz-1);
696
697 //--------------------------------------------------------------------------
698 // Using the equilibrium distribution, calculate the new
699 // distribution function.
700 //--------------------------------------------------------------------------
701 (*this.*update_full)();
702
703 std::swap(f_lb,fnew);
704
705 //--------------------------------------------------------------------------
706 // Calculate moments of the distribution function.
707 //--------------------------------------------------------------------------
708 parametercalc_full();
709
710 //--------------------------------------------------------------------------
711 // Store the equilibrium distribution function, it is needed in
712 // the next time step by the update routine.
713 //--------------------------------------------------------------------------
714 if (typeLB == 2) {
715 std::swap(feqold,feq);
716 std::swap(feqoldn,feqn);
717 }
718
719 //--------------------------------------------------------------------------
720 // Perform diagnostics, and print output for the graphics program
721 //--------------------------------------------------------------------------
722 if (printfluid > 0 && update->ntimestep > 0 && (update->ntimestep % printfluid == 0))
723 streamout();
724
725 }
post_force(int)726 void FixLbFluid::post_force(int /*vflag*/)
727 {
728 if (fixviscouslb==1)
729 calc_fluidforce();
730 }
731
end_of_step()732 void FixLbFluid::end_of_step()
733 {
734 if (fixviscouslb==0)
735 calc_fluidforce();
736
737 if (printrestart>0) {
738 if ((update->ntimestep)%printrestart == 0) {
739 write_restartfile();
740 }
741 }
742
743 }
744
745 //==========================================================================
746 // allocate atom-based array
747 //==========================================================================
grow_arrays(int nmax)748 void FixLbFluid::grow_arrays(int nmax)
749 {
750 memory->grow(hydroF,nmax,3,"FixLbFluid:hydroF");
751 }
752
753 //==========================================================================
754 // copy values within local atom-based array
755 //==========================================================================
copy_arrays(int i,int j,int)756 void FixLbFluid::copy_arrays(int i, int j, int /*delflag*/)
757 {
758 hydroF[j][0] = hydroF[i][0];
759 hydroF[j][1] = hydroF[i][1];
760 hydroF[j][2] = hydroF[i][2];
761 }
762
763 //==========================================================================
764 // pack values in local atom-based array for exchange with another proc
765 //==========================================================================
pack_exchange(int i,double * buf)766 int FixLbFluid::pack_exchange(int i, double *buf)
767 {
768 buf[0] = hydroF[i][0];
769 buf[1] = hydroF[i][1];
770 buf[2] = hydroF[i][2];
771
772 return 3;
773 }
774
775 //==========================================================================
776 // unpack values in local atom-based array from exchange with another proc
777 //==========================================================================
unpack_exchange(int nlocal,double * buf)778 int FixLbFluid::unpack_exchange(int nlocal, double *buf)
779 {
780 hydroF[nlocal][0] = buf[0];
781 hydroF[nlocal][1] = buf[1];
782 hydroF[nlocal][2] = buf[2];
783
784 return 3;
785 }
786
787 //==========================================================================
788 // calculate the force from the local atoms acting on the fluid.
789 //==========================================================================
calc_fluidforce(void)790 void FixLbFluid::calc_fluidforce(void)
791 {
792 int *mask = atom->mask;
793 int nlocal = atom->nlocal;
794 double **x = atom->x;
795 int i,j,k,m;
796 MPI_Request requests[20];
797 double forceloc[3],force[3];
798 double torqueloc[3],torque[3];
799
800 //--------------------------------------------------------------------------
801 // Zero out arrays
802 //--------------------------------------------------------------------------
803 std::fill(&Ff[0][0][0][0],&Ff[0][0][0][0] + (subNbx+3)*(subNby+3)*(subNbz+3)*3,0.0);
804 std::fill(&Fftempx[0][0][0][0],&Fftempx[0][0][0][0] + 5*(subNby+3)*(subNbz+3)*3,0.0);
805 std::fill(&Fftempy[0][0][0][0],&Fftempy[0][0][0][0] + (subNbx+3)*5*(subNbz+3)*3,0.0);
806 std::fill(&Fftempz[0][0][0][0],&Fftempz[0][0][0][0] + (subNbx+3)*(subNby+3)*5*3,0.0);
807
808 forceloc[0] = forceloc[1] = forceloc[2] = 0.0;
809 torqueloc[0] = torqueloc[1] = torqueloc[2] = 0.0;
810
811 for (i=0; i<atom->nmax; i++)
812 for (j=0; j<3; j++)
813 hydroF[i][j] = 0.0;
814
815
816 double unwrap[3];
817 double dx,dy,dz;
818 double massone;
819 imageint *image = atom->image;
820 double *rmass = atom->rmass;
821 double *mass = atom->mass;
822 int *type = atom->type;
823 double sum[4],xcm[4];
824
825 if (force_diagnostic > 0 && update->ntimestep > 0 && (update->ntimestep % force_diagnostic == 0)) {
826 //Calculate the center of mass of the particle group
827 //(needed to calculate the torque).
828 sum[0] = sum[1] = sum[2] = sum[3] = 0.0;
829 for (i=0; i<nlocal; i++) {
830 if (mask[i] & group->bitmask[igroupforce]) {
831
832 domain->unmap(x[i],image[i],unwrap);
833
834 if (rmass) massone = rmass[i];
835 else massone = mass[type[i]];
836
837 sum[0] += unwrap[0]*massone;
838 sum[1] += unwrap[1]*massone;
839 sum[2] += unwrap[2]*massone;
840 sum[3] += massone;
841 }
842 }
843 MPI_Allreduce(&sum[0],&xcm[0],4,MPI_DOUBLE,MPI_SUM,world);
844 xcm[0] = xcm[0]/xcm[3];
845 xcm[1] = xcm[1]/xcm[3];
846 xcm[2] = xcm[2]/xcm[3];
847 }
848
849 //--------------------------------------------------------------------------
850 //Calculate the contribution to the force on the fluid.
851 //--------------------------------------------------------------------------
852 for (i=0; i<nlocal; i++) {
853 if (mask[i] & groupbit) {
854 if (trilinear_stencil==1) {
855 trilinear_interpolation(i);
856 } else {
857 peskin_interpolation(i);
858 }
859
860 if (force_diagnostic > 0 && update->ntimestep > 0 && (update->ntimestep % force_diagnostic == 0)) {
861 if (mask[i] & group->bitmask[igroupforce]) {
862
863 domain->unmap(x[i],image[i],unwrap);
864 dx = unwrap[0] - xcm[0];
865 dy = unwrap[1] - xcm[1];
866 dz = unwrap[2] - xcm[2];
867
868 forceloc[0] += hydroF[i][0];
869 forceloc[1] += hydroF[i][1];
870 forceloc[2] += hydroF[i][2];
871 torqueloc[0] += dy*hydroF[i][2] - dz*hydroF[i][1];
872 torqueloc[1] += dz*hydroF[i][0] - dx*hydroF[i][2];
873 torqueloc[2] += dx*hydroF[i][1] - dy*hydroF[i][0];
874 }
875 }
876 }
877 }
878
879 //--------------------------------------------------------------------------
880 //Communicate the force contributions which lie outside the local processor
881 //sub domain.
882 //--------------------------------------------------------------------------
883 for (i=0; i<10; i++)
884 requests[i]=MPI_REQUEST_NULL;
885 MPI_Isend(&Ff[0][0][0][0],1,passxu,comm->procneigh[0][0],10,world,&requests[0]);
886 MPI_Isend(&Ff[subNbx+2][0][0][0],1,passxu,comm->procneigh[0][0],20,world,&requests[1]);
887 MPI_Isend(&Ff[subNbx-1][0][0][0],1,passxu,comm->procneigh[0][1],30,world,&requests[2]);
888 MPI_Isend(&Ff[subNbx][0][0][0],1,passxu,comm->procneigh[0][1],40,world,&requests[3]);
889 MPI_Isend(&Ff[subNbx+1][0][0][0],1,passxu,comm->procneigh[0][1],50,world,&requests[4]);
890 MPI_Irecv(&Fftempx[0][0][0][0],1,passxtemp,comm->procneigh[0][1],10,world,&requests[5]);
891 MPI_Irecv(&Fftempx[1][0][0][0],1,passxtemp,comm->procneigh[0][1],20,world,&requests[6]);
892 MPI_Irecv(&Fftempx[2][0][0][0],1,passxtemp,comm->procneigh[0][0],30,world,&requests[7]);
893 MPI_Irecv(&Fftempx[3][0][0][0],1,passxtemp,comm->procneigh[0][0],40,world,&requests[8]);
894 MPI_Irecv(&Fftempx[4][0][0][0],1,passxtemp,comm->procneigh[0][0],50,world,&requests[9]);
895 MPI_Waitall(10,requests,MPI_STATUS_IGNORE);
896
897 for (j=0; j<subNby+3; j++) {
898 for (k=0; k<subNbz+3; k++) {
899 for (m=0; m<3; m++) {
900 Ff[subNbx-2][j][k][m] += Fftempx[0][j][k][m];
901 Ff[subNbx-3][j][k][m] += Fftempx[1][j][k][m];
902 Ff[1][j][k][m] += Fftempx[2][j][k][m];
903 Ff[2][j][k][m] += Fftempx[3][j][k][m];
904 Ff[3][j][k][m] += Fftempx[4][j][k][m];
905 }
906 }
907 }
908
909 for (i=0; i<10; i++)
910 requests[i]=MPI_REQUEST_NULL;
911 MPI_Isend(&Ff[0][0][0][0],1,passyu,comm->procneigh[1][0],10,world,&requests[0]);
912 MPI_Isend(&Ff[0][subNby+2][0][0],1,passyu,comm->procneigh[1][0],20,world,&requests[1]);
913 MPI_Isend(&Ff[0][subNby-1][0][0],1,passyu,comm->procneigh[1][1],30,world,&requests[2]);
914 MPI_Isend(&Ff[0][subNby][0][0],1,passyu,comm->procneigh[1][1],40,world,&requests[3]);
915 MPI_Isend(&Ff[0][subNby+1][0][0],1,passyu,comm->procneigh[1][1],50,world,&requests[4]);
916 MPI_Irecv(&Fftempy[0][0][0][0],1,passytemp,comm->procneigh[1][1],10,world,&requests[5]);
917 MPI_Irecv(&Fftempy[0][1][0][0],1,passytemp,comm->procneigh[1][1],20,world,&requests[6]);
918 MPI_Irecv(&Fftempy[0][2][0][0],1,passytemp,comm->procneigh[1][0],30,world,&requests[7]);
919 MPI_Irecv(&Fftempy[0][3][0][0],1,passytemp,comm->procneigh[1][0],40,world,&requests[8]);
920 MPI_Irecv(&Fftempy[0][4][0][0],1,passytemp,comm->procneigh[1][0],50,world,&requests[9]);
921 MPI_Waitall(10,requests,MPI_STATUS_IGNORE);
922
923 for (i=0; i<subNbx+3; i++) {
924 for (k=0; k<subNbz+3; k++) {
925 for (m=0; m<3; m++) {
926 Ff[i][subNby-2][k][m] += Fftempy[i][0][k][m];
927 Ff[i][subNby-3][k][m] += Fftempy[i][1][k][m];
928 Ff[i][1][k][m] += Fftempy[i][2][k][m];
929 Ff[i][2][k][m] += Fftempy[i][3][k][m];
930 Ff[i][3][k][m] += Fftempy[i][4][k][m];
931 }
932 }
933 }
934
935 for (i=0; i<10; i++)
936 requests[i]=MPI_REQUEST_NULL;
937 MPI_Isend(&Ff[0][0][0][0],1,passzu,comm->procneigh[2][0],10,world,&requests[0]);
938 MPI_Isend(&Ff[0][0][subNbz+2][0],1,passzu,comm->procneigh[2][0],20,world,&requests[1]);
939 MPI_Isend(&Ff[0][0][subNbz-1][0],1,passzu,comm->procneigh[2][1],30,world,&requests[2]);
940 MPI_Isend(&Ff[0][0][subNbz][0],1,passzu,comm->procneigh[2][1],40,world,&requests[3]);
941 MPI_Isend(&Ff[0][0][subNbz+1][0],1,passzu,comm->procneigh[2][1],50,world,&requests[4]);
942 MPI_Irecv(&Fftempz[0][0][0][0],1,passztemp,comm->procneigh[2][1],10,world,&requests[5]);
943 MPI_Irecv(&Fftempz[0][0][1][0],1,passztemp,comm->procneigh[2][1],20,world,&requests[6]);
944 MPI_Irecv(&Fftempz[0][0][2][0],1,passztemp,comm->procneigh[2][0],30,world,&requests[7]);
945 MPI_Irecv(&Fftempz[0][0][3][0],1,passztemp,comm->procneigh[2][0],40,world,&requests[8]);
946 MPI_Irecv(&Fftempz[0][0][4][0],1,passztemp,comm->procneigh[2][0],50,world,&requests[9]);
947 MPI_Waitall(10,requests,MPI_STATUS_IGNORE);
948
949 for (i=0; i<subNbx+3; i++) {
950 for (j=0; j<subNby+3; j++) {
951 for (m=0; m<3; m++) {
952 Ff[i][j][subNbz-2][m] += Fftempz[i][j][0][m];
953 Ff[i][j][subNbz-3][m] += Fftempz[i][j][1][m];
954 Ff[i][j][1][m] += Fftempz[i][j][2][m];
955 Ff[i][j][2][m] += Fftempz[i][j][3][m];
956 Ff[i][j][3][m] += Fftempz[i][j][4][m];
957 }
958 }
959 }
960
961 if (force_diagnostic > 0 && update->ntimestep > 0 && (update->ntimestep % force_diagnostic == 0)) {
962 force[0] = force[1] = force[2] = 0.0;
963 torque[0] = torque[1] = torque[2] =0.0;
964
965 MPI_Allreduce(&forceloc[0],&force[0],3,MPI_DOUBLE,MPI_SUM,world);
966 MPI_Allreduce(&torqueloc[0],&torque[0],3,MPI_DOUBLE,MPI_SUM,world);
967
968 if (me==0) {
969 printf("%E %E %E %E %E %E\n",force[0],force[1],force[2],
970 torque[0],torque[1],torque[2]);
971
972 }
973 }
974
975 }
976 //==========================================================================
977 // uses the Peskin stencil to perform the velocity, density and
978 // force interpolations.
979 //==========================================================================
peskin_interpolation(int i)980 void FixLbFluid::peskin_interpolation(int i)
981 {
982 double **x = atom->x;
983 double **v = atom->v;
984 int *type = atom->type;
985 double *rmass = atom->rmass;
986 double *mass = atom->mass;
987 double massone;
988 int ix,iy,iz;
989 int ixp,iyp,izp;
990 double dx1,dy1,dz1;
991 int isten,ii,jj,kk;
992 double r,rsq,weightx,weighty,weightz;
993 double FfP[64];
994 int k;
995 double unode[3];
996 double mnode;
997 double gammavalue;
998
999 //--------------------------------------------------------------------------
1000 //Calculate nearest leftmost grid point.
1001 //Since array indices from 1 to subNb-2 correspond to the
1002 // local subprocessor domain (not indices from 0), use the
1003 // ceiling value.
1004 //--------------------------------------------------------------------------
1005 ix = (int)ceil((x[i][0]-domain->sublo[0])/dx_lb);
1006 iy = (int)ceil((x[i][1]-domain->sublo[1])/dx_lb);
1007 iz = (int)ceil((x[i][2]-domain->sublo[2])/dx_lb);
1008
1009 //--------------------------------------------------------------------------
1010 //Calculate distances to the nearest points.
1011 //--------------------------------------------------------------------------
1012 dx1 = x[i][0] - (domain->sublo[0] + (ix-1)*dx_lb);
1013 dy1 = x[i][1] - (domain->sublo[1] + (iy-1)*dx_lb);
1014 dz1 = x[i][2] - (domain->sublo[2] + (iz-1)*dx_lb);
1015
1016 // Need to convert these to lattice units:
1017 dx1 = dx1/dx_lb;
1018 dy1 = dy1/dx_lb;
1019 dz1 = dz1/dx_lb;
1020
1021 unode[0]=0.0; unode[1]=0.0; unode[2]=0.0;
1022 mnode = 0.0;
1023 isten=0;
1024
1025 //--------------------------------------------------------------------------
1026 // Calculate the interpolation weights, and interpolated values of
1027 // the fluid velocity, and density.
1028 //--------------------------------------------------------------------------
1029 for (ii=-1; ii<3; ii++) {
1030 rsq=(-dx1+ii)*(-dx1+ii);
1031
1032 if (rsq>=4) {
1033 weightx=0.0;
1034 } else {
1035 r=sqrt(rsq);
1036 if (rsq>1) {
1037 weightx=(5.0-2.0*r-sqrt(-7.0+12.0*r-4.0*rsq))/8.;
1038 } else {
1039 weightx=(3.0-2.0*r+sqrt(1.0+4.0*r-4.0*rsq))/8.;
1040 }
1041 }
1042 for (jj=-1; jj<3; jj++) {
1043 rsq=(-dy1+jj)*(-dy1+jj);
1044 if (rsq>=4) {
1045 weighty=0.0;
1046 } else {
1047 r=sqrt(rsq);
1048 if (rsq>1) {
1049 weighty=(5.0-2.0*r-sqrt(-7.0+12.0*r-4.0*rsq))/8.;
1050 } else {
1051 weighty=(3.0-2.0*r+sqrt(1.0+4.0*r-4.0*rsq))/8.;
1052 }
1053 }
1054 for (kk=-1; kk<3; kk++) {
1055 rsq=(-dz1+kk)*(-dz1+kk);
1056 if (rsq>=4) {
1057 weightz=0.0;
1058 } else {
1059 r=sqrt(rsq);
1060 if (rsq>1) {
1061 weightz=(5.0-2.0*r-sqrt(-7.0+12.0*r-4.0*rsq))/8.;
1062 } else {
1063 weightz=(3.0-2.0*r+sqrt(1.0+4.0*r-4.0*rsq))/8.;
1064 }
1065 }
1066 ixp = ix+ii;
1067 iyp = iy+jj;
1068 izp = iz+kk;
1069
1070 //The atom is allowed to be within one lattice grid point outside the
1071 //local processor sub-domain.
1072 if (ixp < -1 || ixp > (subNbx+1) || iyp < -1 || iyp > (subNby+1) || izp < -1 || izp > (subNbz+1))
1073 error->one(FLERR,"Atom outside local processor simulation domain. Either unstable fluid pararmeters, or \
1074 require more frequent neighborlist rebuilds");
1075
1076 if (domain->periodicity[2] == 0 && comm->myloc[2] == 0 && izp < 1)
1077 error->warning(FLERR,"Atom too close to lower z wall. Unphysical results may occur");
1078 if (domain->periodicity[2] == 0 && comm->myloc[2] == (comm->procgrid[2]-1) && (izp > (subNbz-2) ))
1079 error->warning(FLERR,"Atom too close to upper z wall. Unphysical results may occur");
1080
1081 if (ixp==-1) ixp=subNbx+2;
1082 if (iyp==-1) iyp=subNby+2;
1083 if (izp==-1) izp=subNbz+2;
1084
1085 FfP[isten] = weightx*weighty*weightz;
1086 // interpolated velocity based on delta function.
1087 for (k=0; k<3; k++) {
1088 unode[k] += u_lb[ixp][iyp][izp][k]*FfP[isten];
1089 }
1090 if (setGamma==0)
1091 mnode += density_lb[ixp][iyp][izp]*FfP[isten];
1092
1093 isten++;
1094 }
1095 }
1096 }
1097 if (setGamma==0) {
1098 mnode *= NodeArea[type[i]];
1099
1100 if (rmass) massone = rmass[i];
1101 else massone = mass[type[i]];
1102 massone = massone/dm_lb;
1103
1104 gammavalue = 2.0*(mnode*massone)*dtoverdtcollision/(mnode+massone);
1105 } else {
1106 gammavalue = Gamma[type[i]];
1107 }
1108
1109 isten=0;
1110 for (ii=-1; ii<3; ii++)
1111 for (jj=-1; jj<3; jj++)
1112 for (kk=-1; kk<3; kk++) {
1113 ixp = ix+ii;
1114 iyp = iy+jj;
1115 izp = iz+kk;
1116
1117 if (ixp==-1) ixp=subNbx+2;
1118 if (iyp==-1) iyp=subNby+2;
1119 if (izp==-1) izp=subNbz+2;
1120 // Compute the force on the fluid. Need to convert the velocity from
1121 // LAMMPS units to LB units.
1122 for (k=0; k<3; k++) {
1123 Ff[ixp][iyp][izp][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[isten];
1124 }
1125
1126 isten++;
1127 }
1128 for (k=0; k<3; k++)
1129 hydroF[i][k] = -1.0*gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*dm_lb*dx_lb/dt_lb/dt_lb;
1130 }
1131
1132 //==========================================================================
1133 // uses the trilinear stencil to perform the velocity, density and
1134 // force interpolations.
1135 //==========================================================================
trilinear_interpolation(int i)1136 void FixLbFluid::trilinear_interpolation(int i)
1137 {
1138 double **x = atom->x;
1139 double **v = atom->v;
1140 int *type = atom->type;
1141 double *rmass = atom->rmass;
1142 double *mass = atom->mass;
1143 double massone;
1144 int ix,iy,iz;
1145 int ixp,iyp,izp;
1146 double dx1,dy1,dz1;
1147 double FfP[8];
1148 int k;
1149 double unode[3];
1150 double mnode;
1151 double gammavalue;
1152
1153 //--------------------------------------------------------------------------
1154 // Calculate nearest leftmost grid point.
1155 // Since array indices from 1 to subNb-2 correspond to the
1156 // local subprocessor domain (not indices from 0), use the
1157 // ceiling value.
1158 //--------------------------------------------------------------------------
1159 ix = (int)ceil((x[i][0]-domain->sublo[0])/dx_lb);
1160 iy = (int)ceil((x[i][1]-domain->sublo[1])/dx_lb);
1161 iz = (int)ceil((x[i][2]-domain->sublo[2])/dx_lb);
1162
1163 //--------------------------------------------------------------------------
1164 //Calculate distances to the nearest points.
1165 //--------------------------------------------------------------------------
1166 dx1 = x[i][0] - (domain->sublo[0] + (ix-1)*dx_lb);
1167 dy1 = x[i][1] - (domain->sublo[1] + (iy-1)*dx_lb);
1168 dz1 = x[i][2] - (domain->sublo[2] + (iz-1)*dx_lb);
1169
1170 //--------------------------------------------------------------------------
1171 // Need to convert these to lattice units:
1172 //--------------------------------------------------------------------------
1173 dx1 = dx1/dx_lb;
1174 dy1 = dy1/dx_lb;
1175 dz1 = dz1/dx_lb;
1176
1177 //--------------------------------------------------------------------------
1178 // Calculate the interpolation weights
1179 //--------------------------------------------------------------------------
1180 FfP[0] = (1.-dx1)*(1.-dy1)*(1.-dz1);
1181 FfP[1] = (1.-dx1)*(1.-dy1)*dz1;
1182 FfP[2] = (1.-dx1)*dy1*(1.-dz1);
1183 FfP[3] = (1.-dx1)*dy1*dz1;
1184 FfP[4] = dx1*(1.-dy1)*(1.-dz1);
1185 FfP[5] = dx1*(1.-dy1)*dz1;
1186 FfP[6] = dx1*dy1*(1.-dz1);
1187 FfP[7] = dx1*dy1*dz1;
1188
1189 ixp = (ix+1);
1190 iyp = (iy+1);
1191 izp = (iz+1);
1192
1193 //The atom is allowed to be within one lattice grid point outside the
1194 //local processor sub-domain.
1195 if (ix < 0 || ixp > (subNbx+1) || iy < 0 || iyp > (subNby+1) || iz < 0 || izp > (subNbz+1))
1196 error->one(FLERR,"Atom outside local processor simulation domain. Either unstable fluid pararmeters, or \
1197 require more frequent neighborlist rebuilds");
1198
1199 if (domain->periodicity[2] == 0 && comm->myloc[2] == 0 && (iz < 1 || izp < 1))
1200 error->warning(FLERR,"Atom too close to lower z wall. Unphysical results may occur");
1201 if (domain->periodicity[2] == 0 && comm->myloc[2] == (comm->procgrid[2]-1) && (izp > (subNbz-2) || iz > (subNbz-2)))
1202 error->warning(FLERR,"Atom too close to upper z wall. Unphysical results may occur");
1203
1204
1205 for (k=0; k<3; k++) { // tri-linearly interpolated velocity at node
1206 unode[k] = u_lb[ix][iy][iz][k]*FfP[0]
1207 + u_lb[ix][iy][izp][k]*FfP[1]
1208 + u_lb[ix][iyp][iz][k]*FfP[2]
1209 + u_lb[ix][iyp][izp][k]*FfP[3]
1210 + u_lb[ixp][iy][iz][k]*FfP[4]
1211 + u_lb[ixp][iy][izp][k]*FfP[5]
1212 + u_lb[ixp][iyp][iz][k]*FfP[6]
1213 + u_lb[ixp][iyp][izp][k]*FfP[7];
1214 }
1215
1216 if (setGamma==0) {
1217 mnode = density_lb[ix][iy][iz]*FfP[0]
1218 + density_lb[ix][iy][izp]*FfP[1]
1219 + density_lb[ix][iyp][iz]*FfP[2]
1220 + density_lb[ix][iyp][izp]*FfP[3]
1221 + density_lb[ixp][iy][iz]*FfP[4]
1222 + density_lb[ixp][iy][izp]*FfP[5]
1223 + density_lb[ixp][iyp][iz]*FfP[6]
1224 + density_lb[ixp][iyp][izp]*FfP[7];
1225
1226 mnode *= NodeArea[type[i]];
1227
1228 if (rmass) massone = rmass[i];
1229 else massone = mass[type[i]];
1230 massone = massone/dm_lb;
1231
1232 gammavalue = 2.0*(mnode*massone)*dtoverdtcollision/(mnode+massone);
1233 } else {
1234 gammavalue = Gamma[type[i]];
1235 }
1236
1237
1238 for (k=0; k<3; k++) {
1239 Ff[ix][iy][iz][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[0];
1240 Ff[ix][iy][izp][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[1];
1241 Ff[ix][iyp][iz][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[2];
1242 Ff[ix][iyp][izp][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[3];
1243 Ff[ixp][iy][iz][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[4];
1244 Ff[ixp][iy][izp][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[5];
1245 Ff[ixp][iyp][iz][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[6];
1246 Ff[ixp][iyp][izp][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[7];
1247 }
1248
1249 for (k=0; k<3; k++)
1250 hydroF[i][k] = -1.0*gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*dm_lb*dx_lb/dt_lb/dt_lb;
1251
1252 }
1253
1254 //==========================================================================
1255 // read in a fluid restart file. This is only used to restart the
1256 // fluid portion of a LAMMPS simulation.
1257 //==========================================================================
read_restartfile(void)1258 void FixLbFluid::read_restartfile(void)
1259 {
1260 MPI_Status status;
1261 MPI_Datatype realtype;
1262 MPI_Datatype filetype;
1263
1264
1265 int realsizes[4] = {subNbx,subNby,subNbz,numvel};
1266 int realstarts[4] = {1,1,1,0};
1267 int gsizes[4] = {Nbx,Nby,Nbz,numvel};
1268 int lsizes[4] = {subNbx-2,subNby-2,subNbz-2,numvel};
1269 int starts[4] = {comm->myloc[0]*(subNbx-2),comm->myloc[1]*(subNby-2),comm->myloc[2]*(subNbz-2),0};
1270 if (domain->periodicity[2]==0 && comm->myloc[2]==comm->procgrid[2]-1) {
1271 starts[2] = comm->myloc[2]*(subNbz-3);
1272 }
1273
1274 MPI_Type_create_subarray(4,realsizes,lsizes,realstarts,MPI_ORDER_C,MPI_DOUBLE,&realtype);
1275 MPI_Type_commit(&realtype);
1276
1277 MPI_Type_create_subarray(4,gsizes,lsizes,starts,MPI_ORDER_C,MPI_DOUBLE,&filetype);
1278 MPI_Type_commit(&filetype);
1279
1280 MPI_File_set_view(pFileRead,0,MPI_DOUBLE,filetype,(char *) "native",
1281 MPI_INFO_NULL);
1282 MPI_File_seek(pFileRead,0,MPI_SEEK_SET);
1283 MPI_File_read_all(pFileRead,&f_lb[0][0][0][0],1,realtype,&status);
1284 if (typeLB == 2) {
1285 MPI_File_read_all(pFileRead,&feqold[0][0][0][0],1,realtype,&status);
1286 MPI_File_read_all(pFileRead,&feqoldn[0][0][0][0],1,realtype,&status);
1287 }
1288
1289 MPI_Type_free(&realtype);
1290 MPI_Type_free(&filetype);
1291 MPI_File_close(&pFileRead);
1292
1293 }
1294
1295 //==========================================================================
1296 // write a fluid restart file.
1297 //==========================================================================
write_restartfile(void)1298 void FixLbFluid::write_restartfile(void)
1299 {
1300
1301 MPI_File fh;
1302 MPI_Status status;
1303 MPI_Datatype realtype;
1304 MPI_Datatype filetype;
1305
1306 char *hfile;
1307 hfile = new char[32];
1308 sprintf(hfile,"FluidRestart_" BIGINT_FORMAT ".dat",update->ntimestep);
1309
1310 MPI_File_open(world,hfile,MPI_MODE_WRONLY | MPI_MODE_CREATE, MPI_INFO_NULL,&fh);
1311
1312 int realsizes[4] = {subNbx,subNby,subNbz,numvel};
1313 int realstarts[4] = {1,1,1,0};
1314 int gsizes[4] = {Nbx,Nby,Nbz,numvel};
1315 int lsizes[4] = {subNbx-2,subNby-2,subNbz-2,numvel};
1316 int starts[4] = {comm->myloc[0]*(subNbx-2),comm->myloc[1]*(subNby-2),comm->myloc[2]*(subNbz-2),0};
1317 if (domain->periodicity[2]==0 && comm->myloc[2]==comm->procgrid[2]-1) {
1318 starts[2] = comm->myloc[2]*(subNbz-3);
1319 }
1320
1321 MPI_Type_create_subarray(4,realsizes,lsizes,realstarts,MPI_ORDER_C,MPI_DOUBLE,&realtype);
1322 MPI_Type_commit(&realtype);
1323
1324 MPI_Type_create_subarray(4,gsizes,lsizes,starts,MPI_ORDER_C,MPI_DOUBLE,&filetype);
1325 MPI_Type_commit(&filetype);
1326
1327 MPI_File_set_view(fh,0,MPI_DOUBLE,filetype,(char *) "native",MPI_INFO_NULL);
1328 MPI_File_write_all(fh,&f_lb[0][0][0][0],1,realtype,&status);
1329 if (typeLB == 2) {
1330 MPI_File_write_all(fh,&feqold[0][0][0][0],1,realtype,&status);
1331 MPI_File_write_all(fh,&feqoldn[0][0][0][0],1,realtype,&status);
1332 }
1333
1334 MPI_Type_free(&realtype);
1335 MPI_Type_free(&filetype);
1336 MPI_File_close(&fh);
1337 delete [] hfile;
1338
1339 }
1340
1341 //==========================================================================
1342 // rescale the simulation parameters so that dx_lb=dt_lb=dm_lb=1.
1343 // This assumes that all the simulation parameters have been given in
1344 // terms of distance, time and mass units.
1345 //==========================================================================
rescale(void)1346 void FixLbFluid::rescale(void)
1347 {
1348 vwtp = vwtp*dt_lb/dx_lb;
1349 vwbt = vwbt*dt_lb/dx_lb;
1350
1351 bodyforcex = bodyforcex*dt_lb*dt_lb/dx_lb;
1352 bodyforcey = bodyforcey*dt_lb*dt_lb/dx_lb;
1353 bodyforcez = bodyforcez*dt_lb*dt_lb/dx_lb;
1354
1355 tau=(3.0*viscosity/densityinit_real)*dt_lb*dt_lb/dx_lb/dx_lb;
1356 tau /= dt_lb;
1357 if (typeLB==1)
1358 tau = tau + 0.5;
1359
1360 if (setGamma == 0) {
1361 for (int i=0; i<= atom->ntypes; i++) {
1362 NodeArea[i] = NodeArea[i]/dx_lb/dx_lb;
1363 }
1364 } else {
1365 for (int i=0; i<= atom->ntypes; i++) {
1366 Gamma[i] = Gamma[i]*dt_lb/dm_lb;
1367 }
1368 }
1369
1370 densityinit = densityinit_real*dx_lb*dx_lb*dx_lb/dm_lb;
1371
1372 a_0 = a_0_real*dt_lb*dt_lb/(dx_lb*dx_lb);
1373
1374 // Warn if using the D3Q19 model with noise, and a0 is too small.
1375 if (numvel==19 && noisestress==1 && a_0 < 0.2) {
1376 error->warning(FLERR,"Fix lb/fluid WARNING: Chosen value for a0 may be too small. \
1377 Check temperature reproduction.\n");
1378 }
1379
1380 if (noisestress==1) {
1381 if (a_0>0.5555555) {
1382 error->all(FLERR,"Fix lb/fluid ERROR: the Lattice Boltzmann dx and dt need \
1383 to be chosen such that the scaled a_0 < 5/9\n");
1384 }
1385 }
1386
1387 // Courant Condition:
1388 if (a_0 >= 1.0) {
1389 error->all(FLERR,"Fix lb/fluid ERROR: the lattice Boltzmann dx and dt do not \
1390 satisfy the Courant condition.\n");
1391 }
1392
1393 kB = (force->boltz/force->mvv2e)*dt_lb*dt_lb/dx_lb/dx_lb/dm_lb;
1394
1395 if (typeLB==1) {
1396 expminusdtovertau = 0.0;
1397 Dcoeff = 0.0;
1398 namp = 2.0*kB*T*(tau-0.5)/3.0;
1399 noisefactor = 1.0;
1400 if (a_0 <= 0.333333333333333) {
1401 K_0 = 5.17*(0.333333333333333 - a_0);
1402 } else {
1403 K_0 = 2.57*(a_0 - 0.333333333333333);
1404 }
1405 dtoverdtcollision = dt_lb*6.0*viscosity/densityinit_real/dx_lb/dx_lb;
1406 } else if (typeLB==2) {
1407 expminusdtovertau=exp(-1.0/tau);
1408 Dcoeff=(1.0-(1.0-expminusdtovertau)*tau);
1409 namp = 2.0*kB*T/3.;
1410 noisefactor=sqrt((1.0-expminusdtovertau*expminusdtovertau)/
1411 (2.0))/(1.0-expminusdtovertau);
1412 K_0 = 4.5*(1.0/3.0-a_0);
1413 dtoverdtcollision = dt_lb*3.0*viscosity/densityinit_real/dx_lb/dx_lb;
1414 }
1415
1416 }
1417
1418 //==========================================================================
1419 // Set the lattice-Boltzmann velocity vectors and weights for the D3Q15
1420 // model. Initialize the fluid velocity and density.
1421 //==========================================================================
initializeLB15(void)1422 void FixLbFluid::initializeLB15(void)
1423 {
1424 int i,j,k,m;
1425
1426 //velocity vectors.
1427 e[0][0]= 0;
1428 e[0][1]= 0;
1429 e[0][2]= 0;
1430
1431 e[1][0]= 1;
1432 e[1][1]= 0;
1433 e[1][2]= 0;
1434
1435 e[2][0]= 0;
1436 e[2][1]= 1;
1437 e[2][2]= 0;
1438
1439 e[3][0]= -1;
1440 e[3][1]= 0;
1441 e[3][2]= 0;
1442
1443 e[4][0]= 0;
1444 e[4][1]= -1;
1445 e[4][2]= 0;
1446
1447 e[5][0]= 0;
1448 e[5][1]= 0;
1449 e[5][2]= 1;
1450
1451 e[6][0]= 0;
1452 e[6][1]= 0;
1453 e[6][2]= -1;
1454
1455 e[7][0]= 1;
1456 e[7][1]= 1;
1457 e[7][2]= 1;
1458
1459 e[8][0]= -1;
1460 e[8][1]= 1;
1461 e[8][2]= 1;
1462
1463 e[9][0]= -1;
1464 e[9][1]= -1;
1465 e[9][2]= 1;
1466
1467 e[10][0]= 1;
1468 e[10][1]= -1;
1469 e[10][2]= 1;
1470
1471 e[11][0]= 1;
1472 e[11][1]= 1;
1473 e[11][2]= -1;
1474
1475 e[12][0]= -1;
1476 e[12][1]= 1;
1477 e[12][2]= -1;
1478
1479 e[13][0]= -1;
1480 e[13][1]= -1;
1481 e[13][2]= -1;
1482
1483 e[14][0]= 1;
1484 e[14][1]= -1;
1485 e[14][2]= -1;
1486
1487 //weights.
1488 w_lb[0]=2./9.;
1489 w_lb[1]=1./9.;
1490 w_lb[2]=1./9.;
1491 w_lb[3]=1./9.;
1492 w_lb[4]=1./9.;
1493 w_lb[5]=1./9.;
1494 w_lb[6]=1./9.;
1495 w_lb[7]=1./72.;
1496 w_lb[8]=1./72.;
1497 w_lb[9]=1./72.;
1498 w_lb[10]=1./72.;
1499 w_lb[11]=1./72.;
1500 w_lb[12]=1./72.;
1501 w_lb[13]=1./72.;
1502 w_lb[14]=1./72.;
1503
1504 Ng_lb[0]=1.;
1505 Ng_lb[1]=3.;
1506 Ng_lb[2]=3.;
1507 Ng_lb[3]=3.;
1508 Ng_lb[4]=9./2.;
1509 Ng_lb[5]=9./2.;
1510 Ng_lb[6]=9./2.;
1511 Ng_lb[7]=9.;
1512 Ng_lb[8]=9.;
1513 Ng_lb[9]=9.;
1514 Ng_lb[10]=27./2.;
1515 Ng_lb[11]=27./2.;
1516 Ng_lb[12]=27./2.;
1517 Ng_lb[13]=9.;
1518 Ng_lb[14]=1.;
1519
1520 mg_lb[0][0]=1.;mg_lb[0][1]=1.;mg_lb[0][2]=1.;mg_lb[0][3]=1.;mg_lb[0][4]=1.;
1521 mg_lb[0][5]=1.;mg_lb[0][6]=1.;mg_lb[0][7]=1.;mg_lb[0][8]=1.;mg_lb[0][9]=1.;
1522 mg_lb[0][10]=1.;mg_lb[0][11]=1.;mg_lb[0][12]=1.;mg_lb[0][13]=1.;mg_lb[0][14]=1.;
1523 mg_lb[1][0]=0;mg_lb[1][1]=1.;mg_lb[1][2]=0;mg_lb[1][3]=-1.;mg_lb[1][4]=0;
1524 mg_lb[1][5]=0;mg_lb[1][6]=0;mg_lb[1][7]=1.;mg_lb[1][8]=-1.;mg_lb[1][9]=-1.;
1525 mg_lb[1][10]=1.;mg_lb[1][11]=1.;mg_lb[1][12]=-1.;mg_lb[1][13]=-1.;mg_lb[1][14]=1.;
1526 mg_lb[2][0]=0;mg_lb[2][1]=0;mg_lb[2][2]=1.;mg_lb[2][3]=0;mg_lb[2][4]=-1.;
1527 mg_lb[2][5]=0;mg_lb[2][6]=0;mg_lb[2][7]=1.;mg_lb[2][8]=1.;mg_lb[2][9]=-1.;
1528 mg_lb[2][10]=-1.;mg_lb[2][11]=1.;mg_lb[2][12]=1.;mg_lb[2][13]=-1.;mg_lb[2][14]=-1.;
1529 mg_lb[3][0]=0;mg_lb[3][1]=0;mg_lb[3][2]=0;mg_lb[3][3]=0;mg_lb[3][4]=0;
1530 mg_lb[3][5]=1.;mg_lb[3][6]=-1.;mg_lb[3][7]=1.;mg_lb[3][8]=1.;mg_lb[3][9]=1.;
1531 mg_lb[3][10]=1.;mg_lb[3][11]=-1.;mg_lb[3][12]=-1.;mg_lb[3][13]=-1.;mg_lb[3][14]=-1.;
1532 mg_lb[4][0]=-1./3.;mg_lb[4][1]=2./3.;mg_lb[4][2]=-1./3.;mg_lb[4][3]=2./3.;mg_lb[4][4]=-1./3.;
1533 mg_lb[4][5]=-1./3.;mg_lb[4][6]=-1./3.;mg_lb[4][7]=2./3.;mg_lb[4][8]=2./3.;mg_lb[4][9]=2./3.;
1534 mg_lb[4][10]=2./3.;mg_lb[4][11]=2./3.;mg_lb[4][12]=2./3.;mg_lb[4][13]=2./3.;mg_lb[4][14]=2./3.;
1535 mg_lb[5][0]=-1./3.;mg_lb[5][1]=-1./3.;mg_lb[5][2]=2./3.;mg_lb[5][3]=-1./3.;mg_lb[5][4]=2./3.;
1536 mg_lb[5][5]=-1./3.;mg_lb[5][6]=-1./3.;mg_lb[5][7]=2./3.;mg_lb[5][8]=2./3.;mg_lb[5][9]=2./3.;
1537 mg_lb[5][10]=2./3.;mg_lb[5][11]=2./3.;mg_lb[5][12]=2./3.;mg_lb[5][13]=2./3.;mg_lb[5][14]=2./3.;
1538 mg_lb[6][0]=-1./3.;mg_lb[6][1]=-1./3.;mg_lb[6][2]=-1./3.;mg_lb[6][3]=-1./3.;mg_lb[6][4]=-1./3.;
1539 mg_lb[6][5]=2./3.;mg_lb[6][6]=2./3.;mg_lb[6][7]=2./3.;mg_lb[6][8]=2./3.;mg_lb[6][9]=2./3.;
1540 mg_lb[6][10]=2./3.;mg_lb[6][11]=2./3.;mg_lb[6][12]=2./3.;mg_lb[6][13]=2./3.;mg_lb[6][14]=2./3.;
1541 mg_lb[7][0]=0;mg_lb[7][1]=0;mg_lb[7][2]=0;mg_lb[7][3]=0;mg_lb[7][4]=0;
1542 mg_lb[7][5]=0;mg_lb[7][6]=0;mg_lb[7][7]=1;mg_lb[7][8]=-1;mg_lb[7][9]=1;
1543 mg_lb[7][10]=-1;mg_lb[7][11]=1;mg_lb[7][12]=-1;mg_lb[7][13]=1;mg_lb[7][14]=-1;
1544 mg_lb[8][0]=0;mg_lb[8][1]=0;mg_lb[8][2]=0;mg_lb[8][3]=0;mg_lb[8][4]=0;
1545 mg_lb[8][5]=0;mg_lb[8][6]=0;mg_lb[8][7]=1;mg_lb[8][8]=1;mg_lb[8][9]=-1;
1546 mg_lb[8][10]=-1;mg_lb[8][11]=-1;mg_lb[8][12]=-1;mg_lb[8][13]=1;mg_lb[8][14]=1;
1547 mg_lb[9][0]=0;mg_lb[9][1]=0;mg_lb[9][2]=0;mg_lb[9][3]=0;mg_lb[9][4]=0;
1548 mg_lb[9][5]=0;mg_lb[9][6]=0;mg_lb[9][7]=1;mg_lb[9][8]=-1;mg_lb[9][9]=-1;
1549 mg_lb[9][10]=1;mg_lb[9][11]=-1;mg_lb[9][12]=1;mg_lb[9][13]=1;mg_lb[9][14]=-1;
1550 mg_lb[10][0]=0;mg_lb[10][1]=0;mg_lb[10][2]=-1./3.;mg_lb[10][3]=0;mg_lb[10][4]=1./3.;
1551 mg_lb[10][5]=0;mg_lb[10][6]=0;mg_lb[10][7]=2./3.;mg_lb[10][8]=2./3.;mg_lb[10][9]=-2./3.;
1552 mg_lb[10][10]=-2./3.;mg_lb[10][11]=2./3.;mg_lb[10][12]=2./3.;mg_lb[10][13]=-2./3.;mg_lb[10][14]=-2./3.;
1553 mg_lb[11][0]=0;mg_lb[11][1]=0;mg_lb[11][2]=0;mg_lb[11][3]=0;mg_lb[11][4]=0;
1554 mg_lb[11][5]=-1./3.;mg_lb[11][6]=1./3.;mg_lb[11][7]=2./3.;mg_lb[11][8]=2./3.;mg_lb[11][9]=2./3.;
1555 mg_lb[11][10]=2./3.;mg_lb[11][11]=-2./3.;mg_lb[11][12]=-2./3.;mg_lb[11][13]=-2./3.;mg_lb[11][14]=-2./3.;
1556 mg_lb[12][0]=0;mg_lb[12][1]=-1./3.;mg_lb[12][2]=0;mg_lb[12][3]=1./3.;mg_lb[12][4]=0;
1557 mg_lb[12][5]=0;mg_lb[12][6]=0;mg_lb[12][7]=2./3.;mg_lb[12][8]=-2./3.;mg_lb[12][9]=-2./3.;
1558 mg_lb[12][10]=2./3.;mg_lb[12][11]=2./3.;mg_lb[12][12]=-2./3.;mg_lb[12][13]=-2./3.;mg_lb[12][14]=2./3.;
1559 mg_lb[13][0]=0;mg_lb[13][1]=0;mg_lb[13][2]=0;mg_lb[13][3]=0;mg_lb[13][4]=0;
1560 mg_lb[13][5]=0;mg_lb[13][6]=0;mg_lb[13][7]=1;mg_lb[13][8]=-1;mg_lb[13][9]=1;
1561 mg_lb[13][10]=-1;mg_lb[13][11]=-1;mg_lb[13][12]=1;mg_lb[13][13]=-1;mg_lb[13][14]=1;
1562 mg_lb[14][0]=sqrt(2.);mg_lb[14][1]=-1./sqrt(2.);mg_lb[14][2]=-1./sqrt(2.);
1563 mg_lb[14][3]=-1./sqrt(2.);mg_lb[14][4]=-1./sqrt(2.);
1564 mg_lb[14][5]=-1./sqrt(2.);mg_lb[14][6]=-1./sqrt(2.);mg_lb[14][7]=sqrt(2.);
1565 mg_lb[14][8]=sqrt(2.);mg_lb[14][9]=sqrt(2.);
1566 mg_lb[14][10]=sqrt(2.);mg_lb[14][11]=sqrt(2.);mg_lb[14][12]=sqrt(2.);
1567 mg_lb[14][13]=sqrt(2.);mg_lb[14][14]=sqrt(2.);
1568
1569 for (i=0; i<subNbx+3; i++)
1570 for (j=0; j<subNby+3; j++)
1571 for (k=0; k<subNbz+3; k++) {
1572 u_lb[i][j][k][0]=0.0;
1573 u_lb[i][j][k][1]=0.0;
1574 u_lb[i][j][k][2]=0.0;
1575 density_lb[i][j][k] = densityinit;
1576 }
1577 for (i=0; i<subNbx; i++)
1578 for (j=0; j<subNby; j++)
1579 for (k=0; k<subNbz; k++)
1580 for (m=0; m<15; m++)
1581 f_lb[i][j][k][m] = density_lb[i][j][k]/15.0;
1582
1583 }
1584
1585 //==========================================================================
1586 // Set the lattice-Boltzmann velocity vectors and weights for the D3Q19
1587 // model. Initialize the fluid velocity and density.
1588 //==========================================================================
initializeLB19(void)1589 void FixLbFluid::initializeLB19(void)
1590 {
1591 int i,j,k,m;
1592
1593 //velocity vectors.
1594 e[0][0]= 0;
1595 e[0][1]= 0;
1596 e[0][2]= 0;
1597
1598 e[1][0]= 1;
1599 e[1][1]= 0;
1600 e[1][2]= 0;
1601
1602 e[2][0]= 0;
1603 e[2][1]= 1;
1604 e[2][2]= 0;
1605
1606 e[3][0]= -1;
1607 e[3][1]= 0;
1608 e[3][2]= 0;
1609
1610 e[4][0]= 0;
1611 e[4][1]= -1;
1612 e[4][2]= 0;
1613
1614 e[5][0]= 0;
1615 e[5][1]= 0;
1616 e[5][2]= 1;
1617
1618 e[6][0]= 0;
1619 e[6][1]= 0;
1620 e[6][2]= -1;
1621
1622 e[7][0] = 1;
1623 e[7][1] = 1;
1624 e[7][2] = 0;
1625
1626 e[8][0] = 1;
1627 e[8][1] = -1;
1628 e[8][2] = 0;
1629
1630 e[9][0] = -1;
1631 e[9][1] = 1;
1632 e[9][2] = 0;
1633
1634 e[10][0] = -1;
1635 e[10][1] = -1;
1636 e[10][2] = 0;
1637
1638 e[11][0] = 1;
1639 e[11][1] = 0;
1640 e[11][2] = 1;
1641
1642 e[12][0] = 1;
1643 e[12][1] = 0;
1644 e[12][2] = -1;
1645
1646 e[13][0] = -1;
1647 e[13][1] = 0;
1648 e[13][2] = 1;
1649
1650 e[14][0] = -1;
1651 e[14][1] = 0;
1652 e[14][2] = -1;
1653
1654 e[15][0] = 0;
1655 e[15][1] = 1;
1656 e[15][2] = 1;
1657
1658 e[16][0] = 0;
1659 e[16][1] = 1;
1660 e[16][2] = -1;
1661
1662 e[17][0] = 0;
1663 e[17][1] = -1;
1664 e[17][2] = 1;
1665
1666 e[18][0] = 0;
1667 e[18][1] = -1;
1668 e[18][2] = -1;
1669
1670 //weights.
1671 w_lb[0]=1./3.;
1672 w_lb[1]=1./18.;
1673 w_lb[2]=1./18.;
1674 w_lb[3]=1./18.;
1675 w_lb[4]=1./18.;
1676 w_lb[5]=1./18.;
1677 w_lb[6]=1./18.;
1678 w_lb[7]=1./36.;
1679 w_lb[8]=1./36.;
1680 w_lb[9]=1./36.;
1681 w_lb[10]=1./36.;
1682 w_lb[11]=1./36.;
1683 w_lb[12]=1./36.;
1684 w_lb[13]=1./36.;
1685 w_lb[14]=1./36.;
1686 w_lb[15]=1./36.;
1687 w_lb[16]=1./36.;
1688 w_lb[17]=1./36.;
1689 w_lb[18]=1./36.;
1690
1691 Ng_lb[0]=1.;
1692 Ng_lb[1]=3.;
1693 Ng_lb[2]=3.;
1694 Ng_lb[3]=3.;
1695 Ng_lb[4]=9./2.;
1696 Ng_lb[5]=9./2.;
1697 Ng_lb[6]=9./2.;
1698 Ng_lb[7]=9.;
1699 Ng_lb[8]=9.;
1700 Ng_lb[9]=9.;
1701 Ng_lb[10]=27./2.;
1702 Ng_lb[11]=27./2.;
1703 Ng_lb[12]=27./2.;
1704 Ng_lb[13]=18.;
1705 Ng_lb[14]=18.;
1706 Ng_lb[15]=18.;
1707 Ng_lb[16]=162./7.;
1708 Ng_lb[17]=126./5.;
1709 Ng_lb[18]=30.;
1710
1711 mg_lb[0][0] = 1.; mg_lb[0][1] = 1.; mg_lb[0][2] = 1.; mg_lb[0][3] = 1.; mg_lb[0][4] = 1.;
1712 mg_lb[0][5] = 1.; mg_lb[0][6] = 1.; mg_lb[0][7] = 1.; mg_lb[0][8] = 1.; mg_lb[0][9] = 1.;
1713 mg_lb[0][10]= 1.; mg_lb[0][11]= 1.; mg_lb[0][12]= 1.; mg_lb[0][13]= 1.; mg_lb[0][14]= 1.;
1714 mg_lb[0][15]= 1.; mg_lb[0][16]= 1.; mg_lb[0][17]= 1.; mg_lb[0][18]= 1.;
1715
1716 mg_lb[1][0] = 0.; mg_lb[1][1] = 1.; mg_lb[1][2] = 0.; mg_lb[1][3] =-1.; mg_lb[1][4] = 0.;
1717 mg_lb[1][5] = 0.; mg_lb[1][6] = 0.; mg_lb[1][7] = 1.; mg_lb[1][8] = 1.; mg_lb[1][9] =-1.;
1718 mg_lb[1][10]=-1.; mg_lb[1][11]= 1.; mg_lb[1][12]= 1.; mg_lb[1][13]=-1.; mg_lb[1][14]=-1.;
1719 mg_lb[1][15]= 0.; mg_lb[1][16]= 0.; mg_lb[1][17]= 0.; mg_lb[1][18]= 0.;
1720
1721 mg_lb[2][0] = 0.; mg_lb[2][1] = 0.; mg_lb[2][2] = 1.; mg_lb[2][3] = 0.; mg_lb[2][4] =-1.;
1722 mg_lb[2][5] = 0.; mg_lb[2][6] = 0.; mg_lb[2][7] = 1.; mg_lb[2][8] =-1.; mg_lb[2][9] = 1.;
1723 mg_lb[2][10]=-1.; mg_lb[2][11]= 0.; mg_lb[2][12]= 0.; mg_lb[2][13]= 0.; mg_lb[2][14]= 0.;
1724 mg_lb[2][15]= 1.; mg_lb[2][16]= 1.; mg_lb[2][17]=-1.; mg_lb[2][18]=-1.;
1725
1726 mg_lb[3][0] = 0.; mg_lb[3][1] = 0.; mg_lb[3][2] = 0.; mg_lb[3][3] = 0.; mg_lb[3][4] = 0.;
1727 mg_lb[3][5] = 1.; mg_lb[3][6] =-1.; mg_lb[3][7] = 0.; mg_lb[3][8] = 0.; mg_lb[3][9] = 0.;
1728 mg_lb[3][10]= 0.; mg_lb[3][11]= 1.; mg_lb[3][12]=-1.; mg_lb[3][13]= 1.; mg_lb[3][14]=-1.;
1729 mg_lb[3][15]= 1.; mg_lb[3][16]=-1.; mg_lb[3][17]= 1.; mg_lb[3][18]=-1.;
1730
1731 mg_lb[4][0] =-1./3.; mg_lb[4][1] = 2./3.; mg_lb[4][2] =-1./3.; mg_lb[4][3] = 2./3.; mg_lb[4][4] =-1./3.;
1732 mg_lb[4][5] =-1./3.; mg_lb[4][6] =-1./3.; mg_lb[4][7] = 2./3.; mg_lb[4][8] = 2./3.; mg_lb[4][9] = 2./3.;
1733 mg_lb[4][10]= 2./3.; mg_lb[4][11]= 2./3.; mg_lb[4][12]= 2./3.; mg_lb[4][13]= 2./3.; mg_lb[4][14]= 2./3.;
1734 mg_lb[4][15]=-1./3.; mg_lb[4][16]=-1./3.; mg_lb[4][17]=-1./3.; mg_lb[4][18]=-1./3.;
1735
1736 mg_lb[5][0] =-1./3.; mg_lb[5][1] =-1./3.; mg_lb[5][2] = 2./3.; mg_lb[5][3] =-1./3.; mg_lb[5][4] = 2./3.;
1737 mg_lb[5][5] =-1./3.; mg_lb[5][6] =-1./3.; mg_lb[5][7] = 2./3.; mg_lb[5][8] = 2./3.; mg_lb[5][9] = 2./3.;
1738 mg_lb[5][10]= 2./3.; mg_lb[5][11]=-1./3.; mg_lb[5][12]=-1./3.; mg_lb[5][13]=-1./3.; mg_lb[5][14]=-1./3.;
1739 mg_lb[5][15]= 2./3.; mg_lb[5][16]= 2./3.; mg_lb[5][17]= 2./3.; mg_lb[5][18]= 2./3.;
1740
1741 mg_lb[6][0] =-1./3.; mg_lb[6][1] =-1./3.; mg_lb[6][2] =-1./3.; mg_lb[6][3] =-1./3.; mg_lb[6][4] =-1./3.;
1742 mg_lb[6][5] = 2./3.; mg_lb[6][6] = 2./3.; mg_lb[6][7] =-1./3.; mg_lb[6][8] =-1./3.; mg_lb[6][9] =-1./3.;
1743 mg_lb[6][10]=-1./3.; mg_lb[6][11]= 2./3.; mg_lb[6][12]= 2./3.; mg_lb[6][13]= 2./3.; mg_lb[6][14]= 2./3.;
1744 mg_lb[6][15]= 2./3.; mg_lb[6][16]= 2./3.; mg_lb[6][17]= 2./3.; mg_lb[6][18]= 2./3.;
1745
1746 mg_lb[7][0] = 0.; mg_lb[7][1] = 0.; mg_lb[7][2] = 0.; mg_lb[7][3] = 0.; mg_lb[7][4] = 0.;
1747 mg_lb[7][5] = 0.; mg_lb[7][6] = 0.; mg_lb[7][7] = 1.; mg_lb[7][8] =-1.; mg_lb[7][9] =-1.;
1748 mg_lb[7][10]= 1.; mg_lb[7][11]= 0.; mg_lb[7][12]= 0.; mg_lb[7][13]= 0.; mg_lb[7][14]= 0.;
1749 mg_lb[7][15]= 0.; mg_lb[7][16]= 0.; mg_lb[7][17]= 0.; mg_lb[7][18]= 0.;
1750
1751 mg_lb[8][0] = 0.; mg_lb[8][1] = 0.; mg_lb[8][2] = 0.; mg_lb[8][3] = 0.; mg_lb[8][4] = 0.;
1752 mg_lb[8][5] = 0.; mg_lb[8][6] = 0.; mg_lb[8][7] = 0.; mg_lb[8][8] = 0.; mg_lb[8][9] = 0.;
1753 mg_lb[8][10]= 0.; mg_lb[8][11]= 1.; mg_lb[8][12]=-1.; mg_lb[8][13]=-1.; mg_lb[8][14]= 1.;
1754 mg_lb[8][15]= 0.; mg_lb[8][16]= 0.; mg_lb[8][17]= 0.; mg_lb[8][18]= 0.;
1755
1756 mg_lb[9][0] = 0.; mg_lb[9][1] = 0.; mg_lb[9][2] = 0.; mg_lb[9][3] = 0.; mg_lb[9][4] = 0.;
1757 mg_lb[9][5] = 0.; mg_lb[9][6] = 0.; mg_lb[9][7] = 0.; mg_lb[9][8] = 0.; mg_lb[9][9] = 0.;
1758 mg_lb[9][10]= 0.; mg_lb[9][11]= 0.; mg_lb[9][12]= 0.; mg_lb[9][13]= 0.; mg_lb[9][14]= 0.;
1759 mg_lb[9][15]= 1.; mg_lb[9][16]=-1.; mg_lb[9][17]=-1.; mg_lb[9][18]= 1.;
1760
1761 mg_lb[10][0] = 0.; mg_lb[10][1] =-1./3.; mg_lb[10][2] = 0.; mg_lb[10][3] = 1./3.; mg_lb[10][4] = 0.;
1762 mg_lb[10][5] = 0.; mg_lb[10][6] = 0.; mg_lb[10][7] = 2./3.; mg_lb[10][8] = 2./3.; mg_lb[10][9] =-2./3.;
1763 mg_lb[10][10]=-2./3.; mg_lb[10][11]=-1./3.; mg_lb[10][12]=-1./3.; mg_lb[10][13]= 1./3.; mg_lb[10][14]= 1./3.;
1764 mg_lb[10][15]= 0.; mg_lb[10][16]= 0.; mg_lb[10][17]= 0.; mg_lb[10][18]= 0.;
1765
1766 mg_lb[11][0] = 0.; mg_lb[11][1] = 0.; mg_lb[11][2] =-1./3.; mg_lb[11][3] = 0.; mg_lb[11][4] = 1./3.;
1767 mg_lb[11][5] = 0.; mg_lb[11][6] = 0.; mg_lb[11][7] = 2./3.; mg_lb[11][8] =-2./3.; mg_lb[11][9] = 2./3.;
1768 mg_lb[11][10]=-2./3.; mg_lb[11][11]= 0.; mg_lb[11][12]= 0.; mg_lb[11][13]= 0.; mg_lb[11][14]= 0.;
1769 mg_lb[11][15]=-1./3.; mg_lb[11][16]=-1./3.; mg_lb[11][17]= 1./3.; mg_lb[11][18]= 1./3.;
1770
1771 mg_lb[12][0] = 0.; mg_lb[12][1] = 0.; mg_lb[12][2] = 0.; mg_lb[12][3] = 0.; mg_lb[12][4] = 0.;
1772 mg_lb[12][5] =-1./3.; mg_lb[12][6] = 1./3.; mg_lb[12][7] = 0.; mg_lb[12][8] = 0.; mg_lb[12][9] = 0.;
1773 mg_lb[12][10]= 0.; mg_lb[12][11]= 2./3.; mg_lb[12][12]=-2./3.; mg_lb[12][13]= 2./3.; mg_lb[12][14]=-2./3.;
1774 mg_lb[12][15]=-1./3.; mg_lb[12][16]= 1./3.; mg_lb[12][17]=-1./3.; mg_lb[12][18]= 1./3.;
1775
1776 mg_lb[13][0] = 0.; mg_lb[13][1] =-0.5; mg_lb[13][2] = 0.; mg_lb[13][3] = 0.5; mg_lb[13][4] = 0.;
1777 mg_lb[13][5] = 0.; mg_lb[13][6] = 0.; mg_lb[13][7] = 0.; mg_lb[13][8] = 0.; mg_lb[13][9] = 0.;
1778 mg_lb[13][10]= 0.; mg_lb[13][11]= 0.5; mg_lb[13][12]= 0.5; mg_lb[13][13]=-0.5; mg_lb[13][14]=-0.5;
1779 mg_lb[13][15]= 0.; mg_lb[13][16]= 0.; mg_lb[13][17]= 0.; mg_lb[13][18]= 0.;
1780
1781 mg_lb[14][0] = 0.; mg_lb[14][1] = 0.; mg_lb[14][2] = 0.; mg_lb[14][3] = 0.; mg_lb[14][4] = 0.;
1782 mg_lb[14][5] =-0.5; mg_lb[14][6] = 0.5; mg_lb[14][7] = 0.; mg_lb[14][8] = 0.; mg_lb[14][9] = 0.;
1783 mg_lb[14][10]= 0.; mg_lb[14][11]= 0.; mg_lb[14][12]= 0.; mg_lb[14][13]= 0.; mg_lb[14][14]= 0.;
1784 mg_lb[14][15]= 0.5; mg_lb[14][16]=-0.5; mg_lb[14][17]= 0.5; mg_lb[14][18]=-0.5;
1785
1786 mg_lb[15][0] = 0.; mg_lb[15][1] = 0.; mg_lb[15][2] =-0.5; mg_lb[15][3] = 0.; mg_lb[15][4] = 0.5;
1787 mg_lb[15][5] = 0.; mg_lb[15][6] = 0.; mg_lb[15][7] = 0.; mg_lb[15][8] = 0.; mg_lb[15][9] = 0.;
1788 mg_lb[15][10]= 0.; mg_lb[15][11]= 0.; mg_lb[15][12]= 0.; mg_lb[15][13]= 0.; mg_lb[15][14]= 0.;
1789 mg_lb[15][15]= 0.5; mg_lb[15][16]= 0.5; mg_lb[15][17]=-0.5; mg_lb[15][18]=-0.5;
1790
1791 mg_lb[16][0] = 1./18.; mg_lb[16][1] =-5./18.; mg_lb[16][2] =-5./18.; mg_lb[16][3] =-5./18.; mg_lb[16][4] =-5./18.;
1792 mg_lb[16][5] = 2./9.; mg_lb[16][6] = 2./9.; mg_lb[16][7] = 7./18.; mg_lb[16][8] = 7./18.; mg_lb[16][9] = 7./18.;
1793 mg_lb[16][10]= 7./18.; mg_lb[16][11]=-1./9.; mg_lb[16][12]=-1./9.; mg_lb[16][13]=-1./9.; mg_lb[16][14]=-1./9.;
1794 mg_lb[16][15]=-1./9.; mg_lb[16][16]=-1./9.; mg_lb[16][17]=-1./9.; mg_lb[16][18]=-1./9.;
1795
1796 mg_lb[17][0] = 1./14.; mg_lb[17][1] =-5./14.; mg_lb[17][2] = 1./7.; mg_lb[17][3] =-5./14.; mg_lb[17][4] = 1./7.;
1797 mg_lb[17][5] =-3./14.; mg_lb[17][6] =-3./14.; mg_lb[17][7] = 0.; mg_lb[17][8] = 0.; mg_lb[17][9] = 0.;
1798 mg_lb[17][10]= 0.; mg_lb[17][11]= 5./14.; mg_lb[17][12]= 5./14.; mg_lb[17][13]= 5./14.; mg_lb[17][14]= 5./14.;
1799 mg_lb[17][15]=-1./7.; mg_lb[17][16]=-1./7.; mg_lb[17][17]=-1./7.; mg_lb[17][18]=-1./7.;
1800
1801 mg_lb[18][0] = 1./10.; mg_lb[18][1] = 0.; mg_lb[18][2] =-3./10.; mg_lb[18][3] = 0.; mg_lb[18][4] =-3./10.;
1802 mg_lb[18][5] =-3./10.; mg_lb[18][6] =-3./10.; mg_lb[18][7] = 0.; mg_lb[18][8] = 0.; mg_lb[18][9] = 0.;
1803 mg_lb[18][10]= 0.; mg_lb[18][11]= 0.; mg_lb[18][12]= 0.; mg_lb[18][13]= 0.; mg_lb[18][14]= 0.;
1804 mg_lb[18][15]= 3./10.; mg_lb[18][16]= 3./10.; mg_lb[18][17]= 3./10.; mg_lb[18][18]= 3./10.;
1805
1806 for (i=0; i<subNbx+3; i++)
1807 for (j=0; j<subNby+3; j++)
1808 for (k=0; k<subNbz+3; k++) {
1809 u_lb[i][j][k][0]=0.0;
1810 u_lb[i][j][k][1]=0.0;
1811 u_lb[i][j][k][2]=0.0;
1812 density_lb[i][j][k] = densityinit;
1813 }
1814 for (i=0; i<subNbx; i++)
1815 for (j=0; j<subNby; j++)
1816 for (k=0; k<subNbz; k++)
1817 for (m=0; m<19; m++)
1818 f_lb[i][j][k][m] = density_lb[i][j][k]/19.0;
1819
1820 }
1821
1822 //==========================================================================
1823 // Initialize the equilibrium distribution functions
1824 // (this just uses the initial fluid parameters, and assumes no forces).
1825 //==========================================================================
initialize_feq(void)1826 void FixLbFluid::initialize_feq(void)
1827 {
1828 int i,j,k,p;
1829 MPI_Request requests[8];
1830 int numrequests;
1831
1832 // If using the standary LB integrator, do not need to send feqn.
1833 if (typeLB == 1) {
1834 numrequests = 4;
1835 } else {
1836 numrequests = 8;
1837 }
1838
1839 std::fill(&Ff[0][0][0][0],&Ff[0][0][0][0] + (subNbx+3)*(subNby+3)*(subNbz+3)*3,0.0);
1840 std::fill(&Fftempx[0][0][0][0],&Fftempx[0][0][0][0] + 5*(subNby+3)*(subNbz+3)*3,0.0);
1841 std::fill(&Fftempy[0][0][0][0],&Fftempy[0][0][0][0] + (subNbx+3)*5*(subNbz+3)*3,0.0);
1842 std::fill(&Fftempz[0][0][0][0],&Fftempz[0][0][0][0] + (subNbx+3)*(subNby+3)*5*3,0.0);
1843
1844 if (readrestart == 0) {
1845 step=0;
1846
1847 parametercalc_full();
1848 (*this.*equilibriumdist)(1,subNbx-1,1,subNby-1,1,subNbz-1);
1849
1850 for (i=0; i<numrequests; i++)
1851 requests[i]=MPI_REQUEST_NULL;
1852 MPI_Isend(&feq[1][1][1][0],1,passxf,comm->procneigh[0][0],15,world,&requests[0]);
1853 MPI_Irecv(&feq[0][1][1][0],1,passxf,comm->procneigh[0][0],25,world,&requests[1]);
1854 MPI_Isend(&feq[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],25,world,&requests[2]);
1855 MPI_Irecv(&feq[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],15,world,&requests[3]);
1856 if (typeLB == 2) {
1857 MPI_Isend(&feqn[1][1][1][0],1,passxf,comm->procneigh[0][0],10,world,&requests[4]);
1858 MPI_Irecv(&feqn[0][1][1][0],1,passxf,comm->procneigh[0][0],20,world,&requests[5]);
1859 MPI_Isend(&feqn[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],20,world,&requests[6]);
1860 MPI_Irecv(&feqn[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],10,world,&requests[7]);
1861 }
1862 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
1863
1864 for (i=0; i<numrequests; i++)
1865 requests[i]=MPI_REQUEST_NULL;
1866 MPI_Isend(&feq[0][1][1][0],1,passyf,comm->procneigh[1][0],15,world,&requests[0]);
1867 MPI_Irecv(&feq[0][0][1][0],1,passyf,comm->procneigh[1][0],25,world,&requests[1]);
1868 MPI_Isend(&feq[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],25,world,&requests[2]);
1869 MPI_Irecv(&feq[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],15,world,&requests[3]);
1870 if (typeLB == 2) {
1871 MPI_Isend(&feqn[0][1][1][0],1,passyf,comm->procneigh[1][0],10,world,&requests[4]);
1872 MPI_Irecv(&feqn[0][0][1][0],1,passyf,comm->procneigh[1][0],20,world,&requests[5]);
1873 MPI_Isend(&feqn[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],20,world,&requests[6]);
1874 MPI_Irecv(&feqn[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],10,world,&requests[7]);
1875 }
1876 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
1877
1878 for (i=0; i<numrequests; i++)
1879 requests[i]=MPI_REQUEST_NULL;
1880 MPI_Isend(&feq[0][0][1][0],1,passzf,comm->procneigh[2][0],15,world,&requests[0]);
1881 MPI_Irecv(&feq[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&requests[1]);
1882 MPI_Isend(&feq[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&requests[2]);
1883 MPI_Irecv(&feq[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&requests[3]);
1884 if (typeLB == 2) {
1885 MPI_Isend(&feqn[0][0][1][0],1,passzf,comm->procneigh[2][0],10,world,&requests[4]);
1886 MPI_Irecv(&feqn[0][0][0][0],1,passzf,comm->procneigh[2][0],20,world,&requests[5]);
1887 MPI_Isend(&feqn[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],20,world,&requests[6]);
1888 MPI_Irecv(&feqn[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],10,world,&requests[7]);
1889 }
1890 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
1891
1892 //Save feqold.
1893 if (typeLB == 2) {
1894 for (i=0; i<subNbx; i++)
1895 for (j=0; j<subNby; j++)
1896 for (k=0; k<subNbz; k++)
1897 for (p=0; p<numvel; p++) {
1898 feqold[i][j][k][p] = feq[i][j][k][p];
1899 feqoldn[i][j][k][p] = feqn[i][j][k][p];
1900 }
1901 }
1902 } else {
1903 step = 1;
1904
1905 read_restartfile();
1906
1907 if (typeLB == 2) {
1908 for (i=0; i<8; i++)
1909 requests[i]=MPI_REQUEST_NULL;
1910 MPI_Isend(&feqold[1][1][1][0],1,passxf,comm->procneigh[0][0],15,world,&requests[0]);
1911 MPI_Irecv(&feqold[0][1][1][0],1,passxf,comm->procneigh[0][0],25,world,&requests[1]);
1912 MPI_Isend(&feqold[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],25,world,&requests[2]);
1913 MPI_Irecv(&feqold[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],15,world,&requests[3]);
1914 MPI_Isend(&feqoldn[1][1][1][0],1,passxf,comm->procneigh[0][0],10,world,&requests[4]);
1915 MPI_Irecv(&feqoldn[0][1][1][0],1,passxf,comm->procneigh[0][0],20,world,&requests[5]);
1916 MPI_Isend(&feqoldn[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],20,world,&requests[6]);
1917 MPI_Irecv(&feqoldn[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],10,world,&requests[7]);
1918 MPI_Waitall(8,requests,MPI_STATUS_IGNORE);
1919
1920 for (i=0; i<8; i++)
1921 requests[i]=MPI_REQUEST_NULL;
1922 MPI_Isend(&feqold[0][1][1][0],1,passyf,comm->procneigh[1][0],15,world,&requests[0]);
1923 MPI_Irecv(&feqold[0][0][1][0],1,passyf,comm->procneigh[1][0],25,world,&requests[1]);
1924 MPI_Isend(&feqold[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],25,world,&requests[2]);
1925 MPI_Irecv(&feqold[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],15,world,&requests[3]);
1926 MPI_Isend(&feqoldn[0][1][1][0],1,passyf,comm->procneigh[1][0],10,world,&requests[4]);
1927 MPI_Irecv(&feqoldn[0][0][1][0],1,passyf,comm->procneigh[1][0],20,world,&requests[5]);
1928 MPI_Isend(&feqoldn[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],20,world,&requests[6]);
1929 MPI_Irecv(&feqoldn[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],10,world,&requests[7]);
1930 MPI_Waitall(8,requests,MPI_STATUS_IGNORE);
1931
1932 for (i=0; i<8; i++)
1933 requests[i]=MPI_REQUEST_NULL;
1934 MPI_Isend(&feqold[0][0][1][0],1,passzf,comm->procneigh[2][0],15,world,&requests[0]);
1935 MPI_Irecv(&feqold[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&requests[1]);
1936 MPI_Isend(&feqold[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&requests[2]);
1937 MPI_Irecv(&feqold[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&requests[3]);
1938 MPI_Isend(&feqoldn[0][0][1][0],1,passzf,comm->procneigh[2][0],10,world,&requests[4]);
1939 MPI_Irecv(&feqoldn[0][0][0][0],1,passzf,comm->procneigh[2][0],20,world,&requests[5]);
1940 MPI_Isend(&feqoldn[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],20,world,&requests[6]);
1941 MPI_Irecv(&feqoldn[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],10,world,&requests[7]);
1942 MPI_Waitall(8,requests,MPI_STATUS_IGNORE);
1943 }
1944 parametercalc_full();
1945 }
1946 }
1947
1948 //==========================================================================
1949 // Compute the lattice Boltzmann equilibrium distribution functions for
1950 // the D3Q15 model.
1951 //==========================================================================
equilibriumdist15(int xstart,int xend,int ystart,int yend,int zstart,int zend)1952 void FixLbFluid::equilibriumdist15(int xstart, int xend, int ystart, int yend, int zstart, int zend) {
1953
1954 double rho;
1955 int i, j, k, l, iup, idwn, jup, jdwn, kup, kdwn;
1956 double Fx_w, Fy_w, Fz_w;
1957
1958 double total_density(0.0);
1959 double drhox, drhoy, drhoz, drhoxx, drhoyy, drhozz;
1960 double Pxx, Pyy, Pzz, Pxy, Pxz, Pyz;
1961 double grs, p0;
1962 double dPdrho;
1963
1964 double S[2][3],std;
1965 int jj;
1966
1967 double etacov[15],ghostnoise;
1968
1969
1970 for (i=xstart; i<xend; i++) {
1971 iup=i+1;
1972 idwn=i-1;
1973 for (j=ystart; j<yend; j++) {
1974 jup=j+1;
1975 jdwn=j-1;
1976 for (k=zstart; k<zend; k++) {
1977 kup=k+1;
1978 kdwn=k-1;
1979
1980 rho=density_lb[i][j][k];
1981 total_density += rho;
1982
1983 // Derivatives.
1984 drhox = (density_lb[iup][j][k] - density_lb[idwn][j][k])/2.0;
1985 drhoxx = (density_lb[iup][j][k] - 2.0*density_lb[i][j][k] +
1986 density_lb[idwn][j][k]);
1987
1988 drhoy = (density_lb[i][jup][k] - density_lb[i][jdwn][k])/2.0;
1989 drhoyy = (density_lb[i][jup][k] - 2.0*density_lb[i][j][k] +
1990 density_lb[i][jdwn][k]);
1991
1992 drhoz = (density_lb[i][j][kup] - density_lb[i][j][kdwn])/2.0;
1993 drhozz = (density_lb[i][j][kup] - 2.0*density_lb[i][j][k] +
1994 density_lb[i][j][kdwn]);
1995
1996 // Need one-sided derivatives for the boundary of the domain, if fixed boundary
1997 // conditions are used.
1998 if (domain->periodicity[2]==0) {
1999 if (comm->myloc[2]==0 && k==1) {
2000 drhoz = (-3.0*density_lb[i][j][k] + 4.0*density_lb[i][j][k+1] -
2001 density_lb[i][j][k+2])/2.0;
2002 drhozz = (-density_lb[i][j][k+3] + 4.0*density_lb[i][j][k+2] -
2003 5.0*density_lb[i][j][k+1] + 2.0*rho);
2004 }
2005 if (comm->myloc[2]==comm->procgrid[2]-1 && k==subNbz-2) {
2006 drhoz = -(-3.0*density_lb[i][j][k] + 4.0*density_lb[i][j][k-1] -
2007 density_lb[i][j][k-2])/2.0;
2008 drhozz = (-density_lb[i][j][k-3] + 4.0*density_lb[i][j][k-2] -
2009 5.0*density_lb[i][j][k-1] + 2.0*rho);
2010 }
2011 }
2012
2013 grs = drhox*drhox + drhoy*drhoy + drhoz*drhoz;
2014
2015 p0 = rho*a_0-kappa_lb*rho*(drhoxx + drhoyy + drhozz);
2016 // kappa_lb is the square gradient coeff in the pressure tensor
2017
2018 dPdrho = a_0; //assuming here that kappa_lb = 0.
2019
2020
2021 if (typeLB==1) {
2022 Pxx = p0 + kappa_lb*(drhox*drhox - 0.5*grs)+(tau-0.5)*(1.0/3.0-dPdrho)*
2023 (3.0*u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz);
2024 Pyy = p0 + kappa_lb*(drhoy*drhoy - 0.5*grs)+(tau-0.5)*(1.0/3.0-dPdrho)*
2025 (u_lb[i][j][k][0]*drhox+3.0*u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz);
2026 Pzz = p0 + kappa_lb*(drhoz*drhoz - 0.5*grs)+(tau-0.5)*(1.0/3.0-dPdrho)*
2027 (u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+3.0*u_lb[i][j][k][2]*drhoz);
2028 Pxy = kappa_lb*drhox*drhoy+(tau-0.5)*(1.0/3.0-dPdrho)*
2029 (u_lb[i][j][k][0]*drhoy+u_lb[i][j][k][1]*drhox);
2030 Pxz = kappa_lb*drhox*drhoz+(tau-0.5)*(1.0/3.0-dPdrho)*
2031 (u_lb[i][j][k][0]*drhoz+u_lb[i][j][k][2]*drhox);
2032 Pyz = kappa_lb*drhoy*drhoz+(tau-0.5)*(1.0/3.0-dPdrho)*
2033 (u_lb[i][j][k][1]*drhoz+u_lb[i][j][k][2]*drhoy);
2034 } else if (typeLB==2) {
2035 Pxx = p0 + kappa_lb*(drhox*drhox - 0.5*grs)+tau*(1.0/3.0-dPdrho)*
2036 (3.0*u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz);
2037 Pyy = p0 + kappa_lb*(drhoy*drhoy - 0.5*grs)+tau*(1.0/3.0-dPdrho)*
2038 (u_lb[i][j][k][0]*drhox+3.0*u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz);
2039 Pzz = p0 + kappa_lb*(drhoz*drhoz - 0.5*grs)+tau*(1.0/3.0-dPdrho)*
2040 (u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+3.0*u_lb[i][j][k][2]*drhoz);
2041 Pxy = kappa_lb*drhox*drhoy+tau*(1.0/3.0-dPdrho)*
2042 (u_lb[i][j][k][0]*drhoy+u_lb[i][j][k][1]*drhox);
2043 Pxz = kappa_lb*drhox*drhoz+tau*(1.0/3.0-dPdrho)*
2044 (u_lb[i][j][k][0]*drhoz+u_lb[i][j][k][2]*drhox);
2045 Pyz = kappa_lb*drhoy*drhoz+tau*(1.0/3.0-dPdrho)*
2046 (u_lb[i][j][k][1]*drhoz+u_lb[i][j][k][2]*drhoy);
2047 }
2048
2049 Fx_w = Ff[i][j][k][0];
2050 Fy_w = Ff[i][j][k][1];
2051 Fz_w = Ff[i][j][k][2];
2052
2053 etacov[0] = rho;
2054 etacov[1] = rho*u_lb[i][j][k][0] + Fx_w*tau + rho*bodyforcex*tau;
2055 etacov[2] = rho*u_lb[i][j][k][1] + Fy_w*tau + rho*bodyforcey*tau;
2056 etacov[3] = rho*u_lb[i][j][k][2] + Fz_w*tau + rho*bodyforcez*tau;
2057
2058 etacov[4] = Pxx + rho*u_lb[i][j][k][0]*u_lb[i][j][k][0] -rho/3. +
2059 tau*(2.0*u_lb[i][j][k][0]*(Fx_w+rho*bodyforcex));
2060 etacov[5] = Pyy + rho*u_lb[i][j][k][1]*u_lb[i][j][k][1] -rho/3. +
2061 tau*(2.0*u_lb[i][j][k][1]*(Fy_w+rho*bodyforcey));
2062 etacov[6] = Pzz + rho*u_lb[i][j][k][2]*u_lb[i][j][k][2] -rho/3. +
2063 tau*(2.0*u_lb[i][j][k][2]*(Fz_w+rho*bodyforcez));
2064 etacov[7] = Pxy + rho*u_lb[i][j][k][0]*u_lb[i][j][k][1] +
2065 tau*(u_lb[i][j][k][0]*(Fy_w+rho*bodyforcey) + (Fx_w+rho*bodyforcex)*u_lb[i][j][k][1]);
2066 etacov[8] = Pyz + rho*u_lb[i][j][k][1]*u_lb[i][j][k][2] +
2067 tau*(u_lb[i][j][k][1]*(Fz_w+rho*bodyforcez) + (Fy_w+rho*bodyforcey)*u_lb[i][j][k][2]);
2068 etacov[9] = Pxz + rho*u_lb[i][j][k][0]*u_lb[i][j][k][2] +
2069 tau*(u_lb[i][j][k][0]*(Fz_w+rho*bodyforcez) + (Fx_w+rho*bodyforcex)*u_lb[i][j][k][2]);
2070 etacov[10] = 0.0;
2071 etacov[11] = 0.0;
2072 etacov[12] = 0.0;
2073 etacov[13] = rho*u_lb[i][j][k][0]*u_lb[i][j][k][1]*u_lb[i][j][k][2];
2074 const double TrP = Pxx+Pyy+Pzz;
2075 etacov[14] = K_0*(rho-TrP);
2076
2077 for (l=0; l<15; l++) {
2078
2079 feq[i][j][k][l] = 0.0;
2080 for (int ii=0; ii<15; ii++)
2081 feq[i][j][k][l] += w_lb[l]*mg_lb[ii][l]*etacov[ii]*Ng_lb[ii];
2082
2083 if (typeLB == 2) {
2084 feqn[i][j][k][l] = feq[i][j][k][l];
2085 }
2086 }
2087
2088 if (noisestress==1) {
2089 std = sqrt(namp*rho);
2090
2091 for (jj=0; jj<3; jj++)
2092 S[0][jj] = std*random->gaussian();
2093 for (jj=0; jj<3; jj++)
2094 S[1][jj] = std*random->gaussian();
2095
2096 etacov[4] = (S[0][0]*sqrt(3.0-3.0*a_0));
2097 etacov[5] = ((1.0-3.0*a_0)*S[0][0]/sqrt(3.0-3.0*a_0)+
2098 sqrt((8.0-12.0*a_0)/(3.0-3.0*a_0))*S[0][1]);
2099 etacov[6] = ((1.0-3.0*a_0)*S[0][0]/sqrt(3.0-3.0*a_0)+
2100 (2.0-6.0*a_0)*S[0][1]/sqrt((8.0-12.0*a_0)*(3.0-3.0*a_0))+
2101 sqrt((5.0-9.0*a_0)/(2.0-3.0*a_0))*S[0][2]);
2102 etacov[7] = S[1][0];
2103 etacov[8] = S[1][1];
2104 etacov[9] = S[1][2];
2105
2106 for (l=10; l<15; l++) {
2107 etacov[l] = sqrt(9.0*namp*rho/Ng_lb[l])*random->gaussian();
2108 }
2109 etacov[14] += -K_0*(etacov[4]+etacov[5]+etacov[6]); //correction from noise to TrP
2110
2111 for (l=0; l<15; l++) {
2112 ghostnoise = w_lb[l]*
2113 (mg_lb[4][l]*etacov[4]*Ng_lb[4] + mg_lb[5][l]*etacov[5]*Ng_lb[5] +
2114 mg_lb[6][l]*etacov[6]*Ng_lb[6] + mg_lb[7][l]*etacov[7]*Ng_lb[7] +
2115 mg_lb[8][l]*etacov[8]*Ng_lb[8] + mg_lb[9][l]*etacov[9]*Ng_lb[9] +
2116 mg_lb[10][l]*etacov[10]*Ng_lb[10] + mg_lb[11][l]*etacov[11]*Ng_lb[11]
2117 + mg_lb[12][l]*etacov[12]*Ng_lb[12] + mg_lb[13][l]*etacov[13]*Ng_lb[13]
2118 + mg_lb[14][l]*etacov[14]*Ng_lb[14]);
2119 feq[i][j][k][l] += ghostnoise*noisefactor;
2120 }
2121 }
2122 }
2123 }
2124 }
2125 }
2126
2127 //==========================================================================
2128 // Compute the lattice Boltzmann equilibrium distribution functions for
2129 // the D3Q19 model.
2130 //==========================================================================
equilibriumdist19(int xstart,int xend,int ystart,int yend,int zstart,int zend)2131 void FixLbFluid::equilibriumdist19(int xstart, int xend, int ystart, int yend, int zstart, int zend) {
2132
2133 double rho;
2134 int i, j, k, l, iup, idwn, jup, jdwn, kup, kdwn;
2135 double Fx_w, Fy_w, Fz_w;
2136
2137 double total_density(0.0);
2138 double drhox, drhoy, drhoz, drhoxx, drhoyy, drhozz;
2139 double Pxx, Pyy, Pzz, Pxy, Pxz, Pyz;
2140 double grs, p0;
2141 double dPdrho;
2142
2143 double S[2][3],std;
2144 int jj;
2145
2146 double etacov[19],ghostnoise;
2147
2148 for (i=xstart; i<xend; i++) {
2149 iup=i+1;
2150 idwn=i-1;
2151 for (j=ystart; j<yend; j++) {
2152 jup=j+1;
2153 jdwn=j-1;
2154 for (k=zstart; k<zend; k++) {
2155 kup=k+1;
2156 kdwn=k-1;
2157
2158 rho=density_lb[i][j][k];
2159 total_density += rho;
2160
2161 // Derivatives.
2162 drhox = (density_lb[iup][j][k] - density_lb[idwn][j][k])/2.0;
2163 drhoxx = (density_lb[iup][j][k] - 2.0*density_lb[i][j][k] +
2164 density_lb[idwn][j][k]);
2165
2166 drhoy = (density_lb[i][jup][k] - density_lb[i][jdwn][k])/2.0;
2167 drhoyy = (density_lb[i][jup][k] - 2.0*density_lb[i][j][k] +
2168 density_lb[i][jdwn][k]);
2169
2170 drhoz = (density_lb[i][j][kup] - density_lb[i][j][kdwn])/2.0;
2171 drhozz = (density_lb[i][j][kup] - 2.0*density_lb[i][j][k] +
2172 density_lb[i][j][kdwn]);
2173
2174 // Need one-sided derivatives for the boundary of the domain, if fixed boundary
2175 // conditions are used.
2176 if (domain->periodicity[2]==0) {
2177 if (comm->myloc[2]==0 && k==1) {
2178 drhoz = (-3.0*density_lb[i][j][k] + 4.0*density_lb[i][j][k+1] -
2179 density_lb[i][j][k+2])/2.0;
2180 drhozz = (-density_lb[i][j][k+3] + 4.0*density_lb[i][j][k+2] -
2181 5.0*density_lb[i][j][k+1] + 2.0*rho);
2182 }
2183 if (comm->myloc[2]==comm->procgrid[2]-1 && k==subNbz-2) {
2184 drhoz = -(-3.0*density_lb[i][j][k] + 4.0*density_lb[i][j][k-1] -
2185 density_lb[i][j][k-2])/2.0;
2186 drhozz = (-density_lb[i][j][k-3] + 4.0*density_lb[i][j][k-2] -
2187 5.0*density_lb[i][j][k-1] + 2.0*rho);
2188 }
2189 }
2190
2191 grs = drhox*drhox + drhoy*drhoy + drhoz*drhoz;
2192
2193 p0 = rho*a_0-kappa_lb*rho*(drhoxx + drhoyy + drhozz);
2194 // kappa_lb is the square gradient coeff in the pressure tensor
2195
2196 dPdrho = a_0; //assuming here that kappa_lb = 0.
2197
2198
2199 if (typeLB==1) {
2200 Pxx = p0 + kappa_lb*(drhox*drhox - 0.5*grs)+(tau-0.5)*(1.0/3.0-dPdrho)*
2201 (3.0*u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz);
2202 Pyy = p0 + kappa_lb*(drhoy*drhoy - 0.5*grs)+(tau-0.5)*(1.0/3.0-dPdrho)*
2203 (u_lb[i][j][k][0]*drhox+3.0*u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz);
2204 Pzz = p0 + kappa_lb*(drhoz*drhoz - 0.5*grs)+(tau-0.5)*(1.0/3.0-dPdrho)*
2205 (u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+3.0*u_lb[i][j][k][2]*drhoz);
2206 Pxy = kappa_lb*drhox*drhoy+(tau-0.5)*(1.0/3.0-dPdrho)*
2207 (u_lb[i][j][k][0]*drhoy+u_lb[i][j][k][1]*drhox);
2208 Pxz = kappa_lb*drhox*drhoz+(tau-0.5)*(1.0/3.0-dPdrho)*
2209 (u_lb[i][j][k][0]*drhoz+u_lb[i][j][k][2]*drhox);
2210 Pyz = kappa_lb*drhoy*drhoz+(tau-0.5)*(1.0/3.0-dPdrho)*
2211 (u_lb[i][j][k][1]*drhoz+u_lb[i][j][k][2]*drhoy);
2212 } else if (typeLB==2) {
2213 Pxx = p0 + kappa_lb*(drhox*drhox - 0.5*grs)+tau*(1.0/3.0-dPdrho)*
2214 (3.0*u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz);
2215 Pyy = p0 + kappa_lb*(drhoy*drhoy - 0.5*grs)+tau*(1.0/3.0-dPdrho)*
2216 (u_lb[i][j][k][0]*drhox+3.0*u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz);
2217 Pzz = p0 + kappa_lb*(drhoz*drhoz - 0.5*grs)+tau*(1.0/3.0-dPdrho)*
2218 (u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+3.0*u_lb[i][j][k][2]*drhoz);
2219 Pxy = kappa_lb*drhox*drhoy+tau*(1.0/3.0-dPdrho)*
2220 (u_lb[i][j][k][0]*drhoy+u_lb[i][j][k][1]*drhox);
2221 Pxz = kappa_lb*drhox*drhoz+tau*(1.0/3.0-dPdrho)*
2222 (u_lb[i][j][k][0]*drhoz+u_lb[i][j][k][2]*drhox);
2223 Pyz = kappa_lb*drhoy*drhoz+tau*(1.0/3.0-dPdrho)*
2224 (u_lb[i][j][k][1]*drhoz+u_lb[i][j][k][2]*drhoy);
2225 }
2226
2227 Fx_w = Ff[i][j][k][0];
2228 Fy_w = Ff[i][j][k][1];
2229 Fz_w = Ff[i][j][k][2];
2230
2231 etacov[0] = rho;
2232 etacov[1] = rho*u_lb[i][j][k][0] + Fx_w*tau + rho*bodyforcex*tau;
2233 etacov[2] = rho*u_lb[i][j][k][1] + Fy_w*tau + rho*bodyforcey*tau;
2234 etacov[3] = rho*u_lb[i][j][k][2] + Fz_w*tau + rho*bodyforcez*tau;
2235
2236 etacov[4] = Pxx + rho*u_lb[i][j][k][0]*u_lb[i][j][k][0] -rho/3. +
2237 tau*(2.0*u_lb[i][j][k][0]*(Fx_w+rho*bodyforcex));
2238 etacov[5] = Pyy + rho*u_lb[i][j][k][1]*u_lb[i][j][k][1] -rho/3. +
2239 tau*(2.0*u_lb[i][j][k][1]*(Fy_w+rho*bodyforcey));
2240 etacov[6] = Pzz + rho*u_lb[i][j][k][2]*u_lb[i][j][k][2] -rho/3. +
2241 tau*(2.0*u_lb[i][j][k][2]*(Fz_w+rho*bodyforcez));
2242 etacov[7] = Pxy + rho*u_lb[i][j][k][0]*u_lb[i][j][k][1] +
2243 tau*(u_lb[i][j][k][0]*(Fy_w+rho*bodyforcey) + (Fx_w+rho*bodyforcex)*u_lb[i][j][k][1]);
2244 etacov[8] = Pxz + rho*u_lb[i][j][k][0]*u_lb[i][j][k][2] +
2245 tau*(u_lb[i][j][k][0]*(Fz_w+rho*bodyforcez) + (Fx_w+rho*bodyforcex)*u_lb[i][j][k][2]);
2246 etacov[9] = Pyz + rho*u_lb[i][j][k][1]*u_lb[i][j][k][2] +
2247 tau*(u_lb[i][j][k][1]*(Fz_w+rho*bodyforcez) + (Fy_w+rho*bodyforcey)*u_lb[i][j][k][2]);
2248 etacov[10] = 0.0;
2249 etacov[11] = 0.0;
2250 etacov[12] = 0.0;
2251 etacov[13] = 0.0;
2252 etacov[14] = 0.0;
2253 etacov[15] = 0.0;
2254 etacov[16] = 0.0;
2255 etacov[17] = 0.0;
2256 etacov[18] = 0.0;
2257
2258 for (l=0; l<19; l++) {
2259
2260 feq[i][j][k][l] = 0.0;
2261 for (int ii=0; ii<19; ii++)
2262 feq[i][j][k][l] += w_lb[l]*mg_lb[ii][l]*etacov[ii]*Ng_lb[ii];
2263
2264 if (typeLB == 2) {
2265 feqn[i][j][k][l] = feq[i][j][k][l];
2266 }
2267 }
2268
2269 if (noisestress==1) {
2270 std = sqrt(namp*rho);
2271
2272 for (jj=0; jj<3; jj++)
2273 S[0][jj] = std*random->gaussian();
2274 for (jj=0; jj<3; jj++)
2275 S[1][jj] = std*random->gaussian();
2276
2277 etacov[4] = (S[0][0]*sqrt(3.0-3.0*a_0));
2278 etacov[5] = ((1.0-3.0*a_0)*S[0][0]/sqrt(3.0-3.0*a_0)+
2279 sqrt((8.0-12.0*a_0)/(3.0-3.0*a_0))*S[0][1]);
2280 etacov[6] = ((1.0-3.0*a_0)*S[0][0]/sqrt(3.0-3.0*a_0)+
2281 (2.0-6.0*a_0)*S[0][1]/sqrt((8.0-12.0*a_0)*(3.0-3.0*a_0))+
2282 sqrt((5.0-9.0*a_0)/(2.0-3.0*a_0))*S[0][2]);
2283 etacov[7] = S[1][0];
2284 etacov[8] = S[1][1];
2285 etacov[9] = S[1][2];
2286
2287 for (l=10; l<19; l++) {
2288 etacov[l] = sqrt(9.0*namp*rho/Ng_lb[l])*random->gaussian();
2289 }
2290
2291 for (l=0; l<19; l++) {
2292 ghostnoise = w_lb[l]*
2293 (mg_lb[4][l]*etacov[4]*Ng_lb[4] + mg_lb[5][l]*etacov[5]*Ng_lb[5] +
2294 mg_lb[6][l]*etacov[6]*Ng_lb[6] + mg_lb[7][l]*etacov[7]*Ng_lb[7] +
2295 mg_lb[8][l]*etacov[8]*Ng_lb[8] + mg_lb[9][l]*etacov[9]*Ng_lb[9] +
2296 mg_lb[10][l]*etacov[10]*Ng_lb[10] + mg_lb[11][l]*etacov[11]*Ng_lb[11]
2297 + mg_lb[12][l]*etacov[12]*Ng_lb[12] + mg_lb[13][l]*etacov[13]*Ng_lb[13]
2298 + mg_lb[14][l]*etacov[14]*Ng_lb[14] + mg_lb[15][l]*etacov[15]*Ng_lb[15]
2299 + mg_lb[16][l]*etacov[16]*Ng_lb[16] + mg_lb[17][l]*etacov[17]*Ng_lb[17]
2300 + mg_lb[18][l]*etacov[18]*Ng_lb[18]);
2301 feq[i][j][k][l] += ghostnoise*noisefactor;
2302 }
2303 }
2304
2305 }
2306 }
2307 }
2308 }
2309
2310 //==========================================================================
2311 // Calculate the fluid density and velocity over the entire simulation
2312 // domain.
2313 //==========================================================================
parametercalc_full(void)2314 void FixLbFluid::parametercalc_full(void)
2315 {
2316 MPI_Request requests[4];
2317 MPI_Request requests2[12];
2318 int numrequests;
2319 int i;
2320
2321 //--------------------------------------------------------------------------
2322 // send the boundaries of f_lb, as they will be needed later by the update
2323 // routine, and use these to calculate the density and velocity on the
2324 // boundary.
2325 //--------------------------------------------------------------------------
2326 for (i=0; i<4; i++)
2327 requests[i]=MPI_REQUEST_NULL;
2328 MPI_Isend(&f_lb[1][1][1][0],1,passxf,comm->procneigh[0][0],10,world,&requests[0]);
2329 MPI_Irecv(&f_lb[0][1][1][0],1,passxf,comm->procneigh[0][0],20,world,&requests[1]);
2330 MPI_Isend(&f_lb[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],20,world,&requests[2]);
2331 MPI_Irecv(&f_lb[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],10,world,&requests[3]);
2332 parametercalc_part(1,subNbx-1,1,subNby-1,1,subNbz-1);
2333 MPI_Waitall(4,requests,MPI_STATUS_IGNORE);
2334
2335 for (i=0; i<4; i++)
2336 requests[i]=MPI_REQUEST_NULL;
2337 MPI_Isend(&f_lb[0][1][1][0],1,passyf,comm->procneigh[1][0],10,world,&requests[0]);
2338 MPI_Irecv(&f_lb[0][0][1][0],1,passyf,comm->procneigh[1][0],20,world,&requests[1]);
2339 MPI_Isend(&f_lb[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],20,world,&requests[2]);
2340 MPI_Irecv(&f_lb[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],10,world,&requests[3]);
2341 parametercalc_part(0,1,1,subNby-1,1,subNbz-1);
2342 parametercalc_part(subNbx-1,subNbx,1,subNby-1,1,subNbz-1);
2343 MPI_Waitall(4,requests,MPI_STATUS_IGNORE);
2344
2345 for (i=0; i<4; i++)
2346 requests[i]=MPI_REQUEST_NULL;
2347 MPI_Isend(&f_lb[0][0][1][0],1,passzf,comm->procneigh[2][0],10,world,&requests[0]);
2348 MPI_Irecv(&f_lb[0][0][0][0],1,passzf,comm->procneigh[2][0],20,world,&requests[1]);
2349 MPI_Isend(&f_lb[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],20,world,&requests[2]);
2350 MPI_Irecv(&f_lb[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],10,world,&requests[3]);
2351 parametercalc_part(0,subNbx,0,1,1,subNbz-1);
2352 parametercalc_part(0,subNbx,subNby-1,subNby,1,subNbz-1);
2353 MPI_Waitall(4,requests,MPI_STATUS_IGNORE);
2354
2355 parametercalc_part(0,subNbx,0,subNby,0,1);
2356 parametercalc_part(0,subNbx,0,subNby,subNbz-1,subNbz);
2357
2358 //--------------------------------------------------------------------------
2359 // Send the remaining portions of the u array (and density array if Gamma
2360 // is set the default way).
2361 //--------------------------------------------------------------------------
2362 if (setGamma == 0) numrequests = 12;
2363 else numrequests = 6;
2364
2365 for (i=0; i<numrequests; i++)
2366 requests2[i]=MPI_REQUEST_NULL;
2367 MPI_Isend(&u_lb[2][0][0][0],1,passxu,comm->procneigh[0][0],10,world,&requests2[0]);
2368 MPI_Isend(&u_lb[3][0][0][0],1,passxu,comm->procneigh[0][0],20,world,&requests2[1]);
2369 MPI_Isend(&u_lb[subNbx-3][0][0][0],1,passxu,comm->procneigh[0][1],30,world,&requests2[2]);
2370 MPI_Irecv(&u_lb[subNbx][0][0][0],1,passxu,comm->procneigh[0][1],10,world,&requests2[3]);
2371 MPI_Irecv(&u_lb[subNbx+1][0][0][0],1,passxu,comm->procneigh[0][1],20,world,&requests2[4]);
2372 MPI_Irecv(&u_lb[subNbx+2][0][0][0],1,passxu,comm->procneigh[0][0],30,world,&requests2[5]);
2373 if (setGamma==0) {
2374 MPI_Isend(&density_lb[2][0][0],1,passxrho,comm->procneigh[0][0],40,world,&requests2[6]);
2375 MPI_Isend(&density_lb[3][0][0],1,passxrho,comm->procneigh[0][0],50,world,&requests2[7]);
2376 MPI_Isend(&density_lb[subNbx-3][0][0],1,passxrho,comm->procneigh[0][1],60,world,&requests2[8]);
2377 MPI_Irecv(&density_lb[subNbx][0][0],1,passxrho,comm->procneigh[0][1],40,world,&requests2[9]);
2378 MPI_Irecv(&density_lb[subNbx+1][0][0],1,passxrho,comm->procneigh[0][1],50,world,&requests2[10]);
2379 MPI_Irecv(&density_lb[subNbx+2][0][0],1,passxrho,comm->procneigh[0][0],60,world,&requests2[11]);
2380 }
2381 MPI_Waitall(numrequests,requests2,MPI_STATUS_IGNORE);
2382
2383 for (i=0; i<numrequests; i++)
2384 requests2[i]=MPI_REQUEST_NULL;
2385 MPI_Isend(&u_lb[0][2][0][0],1,passyu,comm->procneigh[1][0],10,world,&requests2[0]);
2386 MPI_Isend(&u_lb[0][3][0][0],1,passyu,comm->procneigh[1][0],20,world,&requests2[1]);
2387 MPI_Isend(&u_lb[0][subNby-3][0][0],1,passyu,comm->procneigh[1][1],30,world,&requests2[2]);
2388 MPI_Irecv(&u_lb[0][subNby][0][0],1,passyu,comm->procneigh[1][1],10,world,&requests2[3]);
2389 MPI_Irecv(&u_lb[0][subNby+1][0][0],1,passyu,comm->procneigh[1][1],20,world,&requests2[4]);
2390 MPI_Irecv(&u_lb[0][subNby+2][0][0],1,passyu,comm->procneigh[1][0],30,world,&requests2[5]);
2391 if (setGamma==0) {
2392 MPI_Isend(&density_lb[0][2][0],1,passyrho,comm->procneigh[1][0],40,world,&requests2[6]);
2393 MPI_Isend(&density_lb[0][3][0],1,passyrho,comm->procneigh[1][0],50,world,&requests2[7]);
2394 MPI_Isend(&density_lb[0][subNby-3][0],1,passyrho,comm->procneigh[1][1],60,world,&requests2[8]);
2395 MPI_Irecv(&density_lb[0][subNby][0],1,passyrho,comm->procneigh[1][1],40,world,&requests2[9]);
2396 MPI_Irecv(&density_lb[0][subNby+1][0],1,passyrho,comm->procneigh[1][1],50,world,&requests2[10]);
2397 MPI_Irecv(&density_lb[0][subNby+2][0],1,passyrho,comm->procneigh[1][0],60,world,&requests2[11]);
2398 }
2399 MPI_Waitall(numrequests,requests2,MPI_STATUS_IGNORE);
2400
2401 for (i=0; i<12; i++)
2402 requests2[i]=MPI_REQUEST_NULL;
2403 int requestcount=0;
2404 if (domain->periodicity[2]!=0 || comm->myloc[2] != 0) {
2405 MPI_Isend(&u_lb[0][0][2][0],1,passzu,comm->procneigh[2][0],10,world,&requests2[requestcount]);
2406 MPI_Isend(&u_lb[0][0][3][0],1,passzu,comm->procneigh[2][0],20,world,&requests2[requestcount+1]);
2407 MPI_Irecv(&u_lb[0][0][subNbz+2][0],1,passzu,comm->procneigh[2][0],30,world,&requests2[requestcount+2]);
2408 requestcount=requestcount+3;
2409 if (setGamma==0) {
2410 MPI_Isend(&density_lb[0][0][2],1,passzrho,comm->procneigh[2][0],40,world,&requests2[requestcount]);
2411 MPI_Isend(&density_lb[0][0][3],1,passzrho,comm->procneigh[2][0],50,world,&requests2[requestcount+1]);
2412 MPI_Irecv(&density_lb[0][0][subNbz+2],1,passzrho,comm->procneigh[2][0],60,world,&requests2[requestcount+2]);
2413 requestcount=requestcount+3;
2414 }
2415 }
2416 if (domain->periodicity[2]!=0 || comm->myloc[2] != (comm->procgrid[2]-1)) {
2417 MPI_Isend(&u_lb[0][0][subNbz-3][0],1,passzu,comm->procneigh[2][1],30,world,&requests2[requestcount]);
2418 MPI_Irecv(&u_lb[0][0][subNbz][0],1,passzu,comm->procneigh[2][1],10,world,&requests2[requestcount+1]);
2419 MPI_Irecv(&u_lb[0][0][subNbz+1][0],1,passzu,comm->procneigh[2][1],20,world,&requests2[requestcount+2]);
2420 requestcount=requestcount+3;
2421 if (setGamma==0) {
2422 MPI_Isend(&density_lb[0][0][subNbz-3],1,passzrho,comm->procneigh[2][1],60,world,&requests2[requestcount]);
2423 MPI_Irecv(&density_lb[0][0][subNbz],1,passzrho,comm->procneigh[2][1],40,world,&requests2[requestcount+1]);
2424 MPI_Irecv(&density_lb[0][0][subNbz+1],1,passzrho,comm->procneigh[2][1],50,world,&requests2[requestcount+2]);
2425 requestcount=requestcount+3;
2426 }
2427 }
2428 MPI_Waitall(requestcount,requests2,MPI_STATUS_IGNORE);
2429
2430 }
2431
2432 //==========================================================================
2433 // Calculate the fluid density and velocity over a simulation volume
2434 // specified by xstart,xend; ystart,yend; zstart,zend.
2435 //==========================================================================
parametercalc_part(int xstart,int xend,int ystart,int yend,int zstart,int zend)2436 void FixLbFluid::parametercalc_part(int xstart, int xend, int ystart, int yend, int zstart, int zend)
2437 {
2438 int i,j,k,m;
2439
2440 for (i=xstart; i<xend; i++) {
2441 for (j=ystart; j<yend; j++) {
2442 for (k=zstart; k<zend; k++) {
2443
2444 density_lb[i][j][k]=0.0;
2445 u_lb[i][j][k][0]=0.0;
2446 u_lb[i][j][k][1]=0.0;
2447 u_lb[i][j][k][2]=0.0;
2448 for (m=0; m<numvel; m++) {
2449
2450 density_lb[i][j][k] += f_lb[i][j][k][m];
2451
2452 u_lb[i][j][k][0] += f_lb[i][j][k][m]*e[m][0];
2453 u_lb[i][j][k][1] += f_lb[i][j][k][m]*e[m][1];
2454 u_lb[i][j][k][2] += f_lb[i][j][k][m]*e[m][2];
2455
2456 }
2457
2458 //For the on-lattice wall scheme, need to set this velocity to zero.
2459 if (domain->periodicity[2]==0) {
2460 if (comm->myloc[2]==0) {
2461 if (k==1) {
2462 u_lb[i][j][k][2]=0.0;
2463 }
2464 }
2465 if (comm->myloc[2]==comm->procgrid[2]-1) {
2466 if (k==subNbz-2) {
2467 u_lb[i][j][k][2]=0.0;
2468 }
2469 }
2470
2471 }
2472
2473 u_lb[i][j][k][0]=u_lb[i][j][k][0]/density_lb[i][j][k];
2474 u_lb[i][j][k][1]=u_lb[i][j][k][1]/density_lb[i][j][k];
2475 u_lb[i][j][k][2]=u_lb[i][j][k][2]/density_lb[i][j][k];
2476 }
2477 }
2478 }
2479
2480 }
2481
2482 //==========================================================================
2483 // Update the distribution function over a simulation volume specified
2484 // by xstart,xend; ystart,yend; zstart,zend.
2485 //==========================================================================
update_periodic(int xstart,int xend,int ystart,int yend,int zstart,int zend)2486 void FixLbFluid::update_periodic(int xstart, int xend, int ystart, int yend, int zstart, int zend)
2487 {
2488 int i,j,k,m;
2489 int imod,jmod,kmod,imodm,jmodm,kmodm;
2490
2491 for (i=xstart; i<xend; i++)
2492 for (j=ystart; j<yend; j++)
2493 for (k=zstart; k<zend; k++) {
2494
2495 if (typeLB==1) {
2496 for (m=0; m<numvel; m++) {
2497 imod = i-e[m][0];
2498 jmod = j-e[m][1];
2499 kmod = k-e[m][2];
2500
2501 fnew[i][j][k][m] = f_lb[imod][jmod][kmod][m] + (feq[imod][jmod][kmod][m]-f_lb[imod][jmod][kmod][m])/tau;
2502 }
2503 } else if (typeLB==2) {
2504 for (m=0; m<numvel; m++) {
2505 imod = i-e[m][0];
2506 jmod = j-e[m][1];
2507 kmod = k-e[m][2];
2508
2509 fnew[i][j][k][m] = feq[imod][jmod][kmod][m] + (f_lb[imod][jmod][kmod][m] - feq[imod][jmod][kmod][m])*expminusdtovertau;
2510 }
2511
2512 fnew[i][j][k][0]+=Dcoeff*(feq[i][j][k][0]-feqold[i][j][k][0]);
2513 for (m=1; m<numvel; m++) {
2514 imod = i-e[m][0];
2515 jmod = j-e[m][1];
2516 kmod = k-e[m][2];
2517 imodm = i+e[m][0];
2518 jmodm = j+e[m][1];
2519 kmodm = k+e[m][2];
2520
2521 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) + (0.5-Dcoeff*(tau+0.5))*
2522 (feqn[imodm][jmodm][kmodm][m] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
2523
2524 }
2525 }
2526 }
2527 }
2528
2529 //==========================================================================
2530 // Print the fluid properties to the screen.
2531 //==========================================================================
streamout(void)2532 void FixLbFluid::streamout(void)
2533 {
2534 int i,j,k;
2535 int istart,jstart,kstart;
2536 int iend,jend,kend;
2537 int w,iproc;
2538 int size,sizeloc;
2539 MPI_Request request_send,request_recv;
2540 MPI_Status status;
2541
2542 //--------------------------------------------------------------------------
2543 // **Uncomment in order to test conservation of mass and momentum.
2544 //--------------------------------------------------------------------------
2545 // massloc=0.0;
2546 // momentumloc[0]=momentumloc[1]=momentumloc[2]=0.0;
2547 // for (i=1; i<subNbx-1; i++) {
2548 // for (j=1; j<subNby-1; j++) {
2549 // for (k=1; k<subNbz-1; k++) {
2550 // massloc += density_lb[i][j][k];
2551 // momentumloc[0] += density_lb[i][j][k]*u_lb[i][j][k][0];
2552 // momentumloc[1] += density_lb[i][j][k]*u_lb[i][j][k][1];
2553 // momentumloc[2] += density_lb[i][j][k]*u_lb[i][j][k][2];
2554 // }
2555 // }
2556 // }
2557
2558 // MPI_Allreduce(&massloc,&mass,1,MPI_DOUBLE,MPI_SUM,world);
2559 // MPI_Allreduce(&momentumloc[0],&momentum[0],3,MPI_DOUBLE,MPI_SUM,world);
2560
2561 // if (comm->me==0) {
2562 // printf("%16.12f %16.12f %16.12f %16.12f\n",mass*dm_lb,momentum[0]*dm_lb*dx_lb/dt_lb,momentum[1]*dm_lb*dx_lb/dt_lb,momentum[2]*dm_lb*dx_lb/dt_lb);
2563 // }
2564
2565 sizeloc=(subNbx*subNby*subNbz*4);
2566 MPI_Allreduce(&sizeloc,&size,1,MPI_INT,MPI_MAX,world);
2567
2568 if (me==0) {
2569 for (iproc=0; iproc < comm->nprocs; iproc++) {
2570 if (iproc) {
2571 MPI_Irecv(&buf[0][0][0][0],size,MPI_DOUBLE,iproc,0,world,&request_recv);
2572 MPI_Wait(&request_recv,&status);
2573
2574 istart=static_cast<int> (buf[0][0][0][0]);
2575 jstart=static_cast<int> (buf[0][0][0][1]);
2576 kstart=static_cast<int> (buf[0][0][0][2]);
2577 iend=static_cast<int> (buf[0][0][1][0]);
2578 jend=static_cast<int> (buf[0][0][1][1]);
2579 kend=static_cast<int> (buf[0][0][1][2]);
2580
2581 for (i=istart; i<iend; i++) {
2582 for (j=jstart; j<jend; j++) {
2583 for (k=kstart; k<kend; k++) {
2584 for (w=0; w<4; w++) {
2585 altogether[i][j][k][w]=buf[i-istart+1][j-jstart+1][k-kstart+1][w];
2586 }
2587 }
2588 }
2589 }
2590 } else {
2591 for (i=1; i<subNbx-1; i++) {
2592 for (j=1; j<subNby-1; j++) {
2593 for (k=1; k<subNbz-1; k++) {
2594 altogether[i-1][j-1][k-1][0]=density_lb[i][j][k];
2595 altogether[i-1][j-1][k-1][1]=u_lb[i][j][k][0];
2596 altogether[i-1][j-1][k-1][2]=u_lb[i][j][k][1];
2597 altogether[i-1][j-1][k-1][3]=u_lb[i][j][k][2];
2598 }
2599 }
2600 }
2601 }
2602 }
2603 //i = Nbx/2;
2604 //j = Nby/2;
2605 for (i=0; i<Nbx; i++)
2606 for (j=0; j<Nby; j++)
2607 for (k=0; k<Nbz; k++) {
2608 printf("%16.12f %16.12f %16.12f %16.12f\n",altogether[i][j][k][0]*dm_lb/dx_lb/dx_lb/dx_lb,altogether[i][j][k][1]*dx_lb/dt_lb,altogether[i][j][k][2]*dx_lb/dt_lb,altogether[i][j][k][3]*dx_lb/dt_lb);
2609 }
2610
2611
2612 } else {
2613 istart=comm->myloc[0]*(subNbx-2);
2614 jstart=comm->myloc[1]*(subNby-2);
2615 if (domain->periodicity[2]==0) {
2616 if (comm->myloc[2]==comm->procgrid[2]-1) {
2617 kstart=comm->myloc[2]*(subNbz-3);
2618 } else {
2619 kstart=comm->myloc[2]*(subNbz-2);
2620 }
2621 } else {
2622 kstart=comm->myloc[2]*(subNbz-2);
2623 }
2624 iend=istart+subNbx-2;
2625 jend=jstart+subNby-2;
2626 kend=kstart+subNbz-2;
2627 for (i=0; i<subNbx; i++) {
2628 for (j=0; j<subNby; j++) {
2629 for (k=0; k<subNbz; k++) {
2630 buf[i][j][k][0]=density_lb[i][j][k];
2631 buf[i][j][k][1]=u_lb[i][j][k][0];
2632 buf[i][j][k][2]=u_lb[i][j][k][1];
2633 buf[i][j][k][3]=u_lb[i][j][k][2];
2634 }
2635 }
2636 }
2637 buf[0][0][0][0]=istart;
2638 buf[0][0][0][1]=jstart;
2639 buf[0][0][0][2]=kstart;
2640 buf[0][0][1][0]=iend;
2641 buf[0][0][1][1]=jend;
2642 buf[0][0][1][2]=kend;
2643
2644 MPI_Isend(&buf[0][0][0][0],size,MPI_DOUBLE,0,0,world,&request_send);
2645 MPI_Wait(&request_send,&status);
2646 }
2647
2648 }
2649
2650 //==========================================================================
2651 // Update the distribution functions over the entire simulation domain for
2652 // the D3Q15 model.
2653 //==========================================================================
update_full15(void)2654 void FixLbFluid::update_full15(void)
2655 {
2656
2657 MPI_Request req_send15,req_recv15;
2658 MPI_Request req_send25,req_recv25;
2659 MPI_Request requests[8];
2660 int numrequests;
2661 double tmp1;
2662 MPI_Status status;
2663 double rb;
2664 int i,j,k,m;
2665 int imod,jmod,kmod;
2666 int imodm,jmodm,kmodm;
2667
2668 //--------------------------------------------------------------------------
2669 // If using the standard LB integrator, do not need to send info about feqn.
2670 //--------------------------------------------------------------------------
2671 if (typeLB == 1) {
2672 numrequests = 4;
2673 } else {
2674 numrequests = 8;
2675 }
2676
2677 //--------------------------------------------------------------------------
2678 // Fixed z boundary conditions.
2679 //--------------------------------------------------------------------------
2680 if (domain->periodicity[2]==0) {
2681
2682 for (i=0; i<numrequests; i++)
2683 requests[i]=MPI_REQUEST_NULL;
2684 MPI_Isend(&feq[1][1][1][0],1,passxf,comm->procneigh[0][0],15,world,&requests[0]);
2685 MPI_Irecv(&feq[0][1][1][0],1,passxf,comm->procneigh[0][0],25,world,&requests[1]);
2686 MPI_Isend(&feq[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],25,world,&requests[2]);
2687 MPI_Irecv(&feq[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],15,world,&requests[3]);
2688 if (typeLB == 2) {
2689 MPI_Isend(&feqn[1][1][1][0],1,passxf,comm->procneigh[0][0],10,world,&requests[4]);
2690 MPI_Irecv(&feqn[0][1][1][0],1,passxf,comm->procneigh[0][0],20,world,&requests[5]);
2691 MPI_Isend(&feqn[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],20,world,&requests[6]);
2692 MPI_Irecv(&feqn[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],10,world,&requests[7]);
2693 }
2694 update_periodic(2,subNbx-2,2,subNby-2,2,subNbz-2);
2695 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
2696
2697
2698 for (i=0; i<numrequests; i++)
2699 requests[i]=MPI_REQUEST_NULL;
2700 MPI_Isend(&feq[0][1][1][0],1,passyf,comm->procneigh[1][0],15,world,&requests[0]);
2701 MPI_Irecv(&feq[0][0][1][0],1,passyf,comm->procneigh[1][0],25,world,&requests[1]);
2702 MPI_Isend(&feq[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],25,world,&requests[2]);
2703 MPI_Irecv(&feq[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],15,world,&requests[3]);
2704 if (typeLB == 2) {
2705 MPI_Isend(&feqn[0][1][1][0],1,passyf,comm->procneigh[1][0],10,world,&requests[4]);
2706 MPI_Irecv(&feqn[0][0][1][0],1,passyf,comm->procneigh[1][0],20,world,&requests[5]);
2707 MPI_Isend(&feqn[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],20,world,&requests[6]);
2708 MPI_Irecv(&feqn[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],10,world,&requests[7]);
2709 }
2710 update_periodic(1,2,2,subNby-2,2,subNbz-2);
2711 update_periodic(subNbx-2,subNbx-1,2,subNby-2,2,subNbz-2);
2712 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
2713
2714 for (i=0; i<numrequests; i++)
2715 requests[i]=MPI_REQUEST_NULL;
2716 MPI_Isend(&feq[0][0][1][0],1,passzf,comm->procneigh[2][0],15,world,&requests[0]);
2717 MPI_Irecv(&feq[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&requests[1]);
2718 MPI_Isend(&feq[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&requests[2]);
2719 MPI_Irecv(&feq[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&requests[3]);
2720 if (typeLB == 2) {
2721 MPI_Isend(&feqn[0][0][1][0],1,passzf,comm->procneigh[2][0],10,world,&requests[4]);
2722 MPI_Irecv(&feqn[0][0][0][0],1,passzf,comm->procneigh[2][0],20,world,&requests[5]);
2723 MPI_Isend(&feqn[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],20,world,&requests[6]);
2724 MPI_Irecv(&feqn[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],10,world,&requests[7]);
2725 }
2726 update_periodic(1,subNbx-1,1,2,2,subNbz-2);
2727 update_periodic(1,subNbx-1,subNby-2,subNby-1,2,subNbz-2);
2728 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
2729
2730 if (typeLB==1) {
2731 update_periodic(1,subNbx-1,1,subNby-1,1,2);
2732 update_periodic(1,subNbx-1,1,subNby-1,subNbz-2,subNbz-1);
2733 } else if (typeLB==2) {
2734 if (comm->myloc[2]==0) {
2735 for (i=1; i<subNbx-1; i++) {
2736 for (j=1;j<subNby-1;j++) {
2737 k=1;
2738 for (m=0; m<15; m++) {
2739 imod = i-e[m][0];
2740 jmod = j-e[m][1];
2741 kmod = k-e[m][2];
2742
2743 fnew[i][j][k][m] = feq[imod][jmod][kmod][m] + (f_lb[imod][jmod][kmod][m]-feq[imod][jmod][kmod][m])*expminusdtovertau;
2744 }
2745
2746 for (m=0; m<15; m++) {
2747 imod = i-e[m][0];
2748 jmod = j-e[m][1];
2749 kmod = k-e[m][2];
2750 imodm = i+e[m][0];
2751 jmodm = j+e[m][1];
2752 kmodm = k+e[m][2];
2753
2754 if (m==5)
2755 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][6] - feqold[imod][jmod][kmod][m]) +
2756 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][6] - feqn[imod][jmod][kmod][6]);
2757 else if (m==7)
2758 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][11] - feqold[imod][jmod][kmod][m]) +
2759 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][11] - feqn[imod][jmod][kmod][11]);
2760 else if (m==8)
2761 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][12] - feqold[imod][jmod][kmod][m]) +
2762 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][12] - feqn[imod][jmod][kmod][12]);
2763 else if (m==9)
2764 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][13] - feqold[imod][jmod][kmod][m]) +
2765 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][13] - feqn[imod][jmod][kmod][13]);
2766 else if (m==10)
2767 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][14] - feqold[imod][jmod][kmod][m]) +
2768 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][14] - feqn[imod][jmod][kmod][14]);
2769 else if (m==6)
2770 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m]-feqold[imod][jmod][kmod][m]) +
2771 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][5] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
2772 else if (m==11)
2773 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m]-feqold[imod][jmod][kmod][m]) +
2774 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][7] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
2775 else if (m==12)
2776 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m]-feqold[imod][jmod][kmod][m]) +
2777 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][8] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
2778 else if (m==13)
2779 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m]-feqold[imod][jmod][kmod][m]) +
2780 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][9] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
2781 else if (m==14)
2782 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m]-feqold[imod][jmod][kmod][m]) +
2783 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][10] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
2784 else
2785 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
2786 (0.5-Dcoeff*(tau+0.5))*(feqn[imodm][jmodm][kmodm][m] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
2787
2788 }
2789 }
2790 }
2791 } else {
2792 update_periodic(1,subNbx-1,1,subNby-1,1,2);
2793 }
2794 if (comm->myloc[2]==comm->procgrid[2]-1) {
2795 for (i=1;i<subNbx-1;i++) {
2796 for (j=1;j<subNby-1;j++) {
2797 k=subNbz-2;
2798 for (m=0; m<15; m++) {
2799 imod = i-e[m][0];
2800 jmod = j-e[m][1];
2801 kmod = k-e[m][2];
2802
2803 fnew[i][j][k][m] = feq[imod][jmod][kmod][m] + (f_lb[imod][jmod][kmod][m]-feq[imod][jmod][kmod][m])*expminusdtovertau;
2804 }
2805 for (m=0; m<15; m++) {
2806 imod = i-e[m][0];
2807 jmod = j-e[m][1];
2808 kmod = k-e[m][2];
2809 imodm = i+e[m][0];
2810 jmodm = j+e[m][1];
2811 kmodm = k+e[m][2];
2812
2813 if (m==6)
2814 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][5] - feqold[imod][jmod][kmod][m]) +
2815 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][5] - feqn[imod][jmod][kmod][5]);
2816 else if (m==11)
2817 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][7] - feqold[imod][jmod][kmod][m]) +
2818 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][7] - feqn[imod][jmod][kmod][7]);
2819 else if (m==12)
2820 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][8] - feqold[imod][jmod][kmod][m]) +
2821 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][8] - feqn[imod][jmod][kmod][8]);
2822 else if (m==13)
2823 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][9] - feqold[imod][jmod][kmod][m]) +
2824 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][9] - feqn[imod][jmod][kmod][9]);
2825 else if (m==14)
2826 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][10] - feqold[imod][jmod][kmod][m]) +
2827 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][10] - feqn[imod][jmod][kmod][10]);
2828 else if (m==5)
2829 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
2830 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][6] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
2831 else if (m==7)
2832 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
2833 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][11] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
2834 else if (m==8)
2835 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
2836 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][12] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
2837 else if (m==9)
2838 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
2839 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][13] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
2840 else if (m==10)
2841 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
2842 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][14] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
2843 else
2844 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
2845 (0.5-Dcoeff*(tau+0.5))*(feqn[imodm][jmodm][kmodm][m] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
2846
2847 }
2848 }
2849 }
2850 } else {
2851 update_periodic(1,subNbx-1,1,subNby-1,subNbz-2,subNbz-1);
2852 }
2853 }
2854
2855 req_send15=MPI_REQUEST_NULL;
2856 req_recv25=MPI_REQUEST_NULL;
2857 req_send25=MPI_REQUEST_NULL;
2858 req_recv15=MPI_REQUEST_NULL;
2859
2860 if (comm->myloc[2]==0) {
2861 MPI_Isend(&fnew[0][0][1][0],1,passzf,comm->procneigh[2][0],15,world,&req_send15);
2862 MPI_Irecv(&fnew[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&req_recv25);
2863 }
2864
2865 if (comm->myloc[2]==comm->procgrid[2]-1) {
2866 MPI_Isend(&fnew[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&req_send25);
2867 MPI_Irecv(&fnew[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&req_recv15);
2868 }
2869 if (comm->myloc[2]==0) {
2870 MPI_Wait(&req_send15,&status);
2871 MPI_Wait(&req_recv25,&status);
2872
2873 for (i=1;i<subNbx-1;i++) {
2874 for (j=1;j<subNby-1;j++) {
2875 k=1;
2876 if (typeLB == 1) {
2877 fnew[i][j][k][5]=fnew[i][j][k-1][6];
2878 tmp1=fnew[i][j][k-1][11]+fnew[i][j][k-1][12]+fnew[i][j][k-1][13]+fnew[i][j][k-1][14];
2879 } else {
2880 fnew[i][j][k][5]=fnew[i][j][k-1][6] + (0.5-Dcoeff*(tau+0.5))*feqn[i][j][k+1][5];
2881 tmp1=fnew[i][j][k-1][11]+fnew[i][j][k-1][12]+fnew[i][j][k-1][13]+fnew[i][j][k-1][14] +
2882 (0.5-Dcoeff*(tau+0.5))*(feqn[i-1][j-1][k+1][7] + feqn[i+1][j-1][k+1][8] +
2883 feqn[i+1][j+1][k+1][9] + feqn[i-1][j+1][k+1][10]);
2884 }
2885
2886 fnew[i][j][k][7]=-0.25*(fnew[i][j][k][1]+fnew[i][j][k][2]-fnew[i][j][k][3]-
2887 fnew[i][j][k][4]+2.0*fnew[i][j][k][11]-2.0*fnew[i][j][k][13]-tmp1);
2888 fnew[i][j][k][8]=0.25*(fnew[i][j][k][1]-fnew[i][j][k][2]-fnew[i][j][k][3]+
2889 fnew[i][j][k][4]+2.0*fnew[i][j][k][14]-2.0*fnew[i][j][k][12]+tmp1);
2890 fnew[i][j][k][9]=0.25*(fnew[i][j][k][1]+fnew[i][j][k][2]-fnew[i][j][k][3]-
2891 fnew[i][j][k][4]+2.0*fnew[i][j][k][11]-2.0*fnew[i][j][k][13]+tmp1);
2892 fnew[i][j][k][10]=-0.25*(fnew[i][j][k][1]-fnew[i][j][k][2]-fnew[i][j][k][3]+
2893 fnew[i][j][k][4]+2.0*fnew[i][j][k][14]-2.0*fnew[i][j][k][12]-tmp1);
2894
2895
2896
2897 rb=fnew[i][j][k][0]+fnew[i][j][k][1]+fnew[i][j][k][2]+fnew[i][j][k][3]+fnew[i][j][k][4]+
2898 fnew[i][j][k][5]+fnew[i][j][k][6]+tmp1+fnew[i][j][k][11]+fnew[i][j][k][12]+
2899 fnew[i][j][k][13]+fnew[i][j][k][14];
2900
2901 fnew[i][j][k][7] += 0.25*rb*vwbt;
2902 fnew[i][j][k][8] += 0.25*rb*vwbt;
2903 fnew[i][j][k][9] += -0.25*rb*vwbt;
2904 fnew[i][j][k][10] += -0.25*rb*vwbt;
2905 }
2906 }
2907
2908 }
2909 if (comm->myloc[2]==comm->procgrid[2]-1) {
2910 MPI_Wait(&req_send25,&status);
2911 MPI_Wait(&req_recv15,&status);
2912
2913 for (i=1;i<subNbx-1;i++) {
2914 for (j=1;j<subNby-1;j++) {
2915 k=subNbz-2;
2916
2917 if (typeLB == 1) {
2918 fnew[i][j][k][6]=fnew[i][j][k+1][5];
2919 tmp1=fnew[i][j][k+1][7]+fnew[i][j][k+1][8]+fnew[i][j][k+1][9]+fnew[i][j][k+1][10];
2920 } else {
2921 fnew[i][j][k][6]=fnew[i][j][k+1][5] + (0.5-Dcoeff*(tau+0.5))*feqn[i][j][k-1][6];
2922 tmp1=fnew[i][j][k+1][7]+fnew[i][j][k+1][8]+fnew[i][j][k+1][9]+fnew[i][j][k+1][10] +
2923 (0.5-Dcoeff*(tau+0.5))*(feqn[i-1][j-1][k-1][11] + feqn[i+1][j-1][k-1][12] +
2924 feqn[i+1][j+1][k-1][13] + feqn[i-1][j+1][k-1][14]);
2925 }
2926
2927 fnew[i][j][k][11]=-0.25*(fnew[i][j][k][1]+fnew[i][j][k][2]-fnew[i][j][k][3]-
2928 fnew[i][j][k][4]+2.0*fnew[i][j][k][7]-2.0*fnew[i][j][k][9]-tmp1);
2929 fnew[i][j][k][12]=0.25*(fnew[i][j][k][1]-fnew[i][j][k][2]-fnew[i][j][k][3]+
2930 fnew[i][j][k][4]-2.0*fnew[i][j][k][8]+2.0*fnew[i][j][k][10]+tmp1);
2931 fnew[i][j][k][13]=0.25*(fnew[i][j][k][1]+fnew[i][j][k][2]-fnew[i][j][k][3]-
2932 fnew[i][j][k][4]+2.0*fnew[i][j][k][7]-2.0*fnew[i][j][k][9]+tmp1);
2933 fnew[i][j][k][14]=-0.25*(fnew[i][j][k][1]-fnew[i][j][k][2]-fnew[i][j][k][3]+
2934 fnew[i][j][k][4]-2.0*fnew[i][j][k][8]+2.0*fnew[i][j][k][10]-tmp1);
2935
2936
2937 rb=fnew[i][j][k][0]+fnew[i][j][k][1]+fnew[i][j][k][2]+fnew[i][j][k][3]+fnew[i][j][k][4]+
2938 fnew[i][j][k][5]+fnew[i][j][k][6]+fnew[i][j][k][7]+fnew[i][j][k][8]+fnew[i][j][k][9]+
2939 fnew[i][j][k][10]+tmp1;
2940
2941 fnew[i][j][k][11] += 0.25*rb*vwtp;
2942 fnew[i][j][k][12] += 0.25*rb*vwtp;
2943 fnew[i][j][k][13] += -0.25*rb*vwtp;
2944 fnew[i][j][k][14] += -0.25*rb*vwtp;
2945 }
2946 }
2947 }
2948
2949 //--------------------------------------------------------------------------
2950 // Periodic z boundary conditions.
2951 //--------------------------------------------------------------------------
2952 } else {
2953
2954 for (i=0; i<numrequests; i++)
2955 requests[i]=MPI_REQUEST_NULL;
2956 MPI_Isend(&feq[1][1][1][0],1,passxf,comm->procneigh[0][0],15,world,&requests[0]);
2957 MPI_Irecv(&feq[0][1][1][0],1,passxf,comm->procneigh[0][0],25,world,&requests[1]);
2958 MPI_Isend(&feq[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],25,world,&requests[2]);
2959 MPI_Irecv(&feq[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],15,world,&requests[3]);
2960 if (typeLB == 2) {
2961 MPI_Isend(&feqn[1][1][1][0],1,passxf,comm->procneigh[0][0],10,world,&requests[4]);
2962 MPI_Irecv(&feqn[0][1][1][0],1,passxf,comm->procneigh[0][0],20,world,&requests[5]);
2963 MPI_Isend(&feqn[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],20,world,&requests[6]);
2964 MPI_Irecv(&feqn[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],10,world,&requests[7]);
2965 }
2966 update_periodic(2,subNbx-2,2,subNby-2,2,subNbz-2);
2967 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
2968
2969 for (i=0; i<numrequests; i++)
2970 requests[i]=MPI_REQUEST_NULL;
2971 MPI_Isend(&feq[0][1][1][0],1,passyf,comm->procneigh[1][0],15,world,&requests[0]);
2972 MPI_Irecv(&feq[0][0][1][0],1,passyf,comm->procneigh[1][0],25,world,&requests[1]);
2973 MPI_Isend(&feq[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],25,world,&requests[2]);
2974 MPI_Irecv(&feq[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],15,world,&requests[3]);
2975 if (typeLB == 2) {
2976 MPI_Isend(&feqn[0][1][1][0],1,passyf,comm->procneigh[1][0],10,world,&requests[4]);
2977 MPI_Irecv(&feqn[0][0][1][0],1,passyf,comm->procneigh[1][0],20,world,&requests[5]);
2978 MPI_Isend(&feqn[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],20,world,&requests[6]);
2979 MPI_Irecv(&feqn[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],10,world,&requests[7]);
2980 }
2981 update_periodic(1,2,2,subNby-2,2,subNbz-2);
2982 update_periodic(subNbx-2,subNbx-1,2,subNby-2,2,subNbz-2);
2983 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
2984
2985 for (i=0; i<numrequests; i++)
2986 requests[i]=MPI_REQUEST_NULL;
2987 MPI_Isend(&feq[0][0][1][0],1,passzf,comm->procneigh[2][0],15,world,&requests[0]);
2988 MPI_Irecv(&feq[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&requests[1]);
2989 MPI_Isend(&feq[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&requests[2]);
2990 MPI_Irecv(&feq[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&requests[3]);
2991 if (typeLB == 2) {
2992 MPI_Isend(&feqn[0][0][1][0],1,passzf,comm->procneigh[2][0],10,world,&requests[4]);
2993 MPI_Irecv(&feqn[0][0][0][0],1,passzf,comm->procneigh[2][0],20,world,&requests[5]);
2994 MPI_Isend(&feqn[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],20,world,&requests[6]);
2995 MPI_Irecv(&feqn[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],10,world,&requests[7]);
2996 }
2997 update_periodic(1,subNbx-1,1,2,2,subNbz-2);
2998 update_periodic(1,subNbx-1,subNby-2,subNby-1,2,subNbz-2);
2999 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
3000
3001 update_periodic(1,subNbx-1,1,subNby-1,1,2);
3002 update_periodic(1,subNbx-1,1,subNby-1,subNbz-2,subNbz-1);
3003 }
3004
3005 }
3006
3007 //==========================================================================
3008 // Update the distribution functions over the entire simulation domain for
3009 // the D3Q19 model.
3010 //==========================================================================
update_full19(void)3011 void FixLbFluid::update_full19(void)
3012 {
3013
3014 MPI_Request req_send15,req_recv15;
3015 MPI_Request req_send25,req_recv25;
3016 MPI_Request requests[8];
3017 int numrequests;
3018 double tmp1,tmp2,tmp3;
3019 MPI_Status status;
3020 double rb;
3021 int i,j,k,m;
3022 int imod,jmod,kmod;
3023 int imodm,jmodm,kmodm;
3024
3025 //--------------------------------------------------------------------------
3026 // If using the standard LB integrator, do not need to send info about feqn.
3027 //--------------------------------------------------------------------------
3028 if (typeLB == 1) {
3029 numrequests = 4;
3030 } else {
3031 numrequests = 8;
3032 }
3033
3034 //--------------------------------------------------------------------------
3035 // Fixed z boundary conditions.
3036 //--------------------------------------------------------------------------
3037 if (domain->periodicity[2]==0) {
3038
3039 for (i=0; i<numrequests; i++)
3040 requests[i]=MPI_REQUEST_NULL;
3041 MPI_Isend(&feq[1][1][1][0],1,passxf,comm->procneigh[0][0],15,world,&requests[0]);
3042 MPI_Irecv(&feq[0][1][1][0],1,passxf,comm->procneigh[0][0],25,world,&requests[1]);
3043 MPI_Isend(&feq[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],25,world,&requests[2]);
3044 MPI_Irecv(&feq[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],15,world,&requests[3]);
3045 if (typeLB == 2) {
3046 MPI_Isend(&feqn[1][1][1][0],1,passxf,comm->procneigh[0][0],10,world,&requests[4]);
3047 MPI_Irecv(&feqn[0][1][1][0],1,passxf,comm->procneigh[0][0],20,world,&requests[5]);
3048 MPI_Isend(&feqn[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],20,world,&requests[6]);
3049 MPI_Irecv(&feqn[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],10,world,&requests[7]);
3050 }
3051 update_periodic(2,subNbx-2,2,subNby-2,2,subNbz-2);
3052 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
3053
3054 for (i=0; i<numrequests; i++)
3055 requests[i]=MPI_REQUEST_NULL;
3056 MPI_Isend(&feq[0][1][1][0],1,passyf,comm->procneigh[1][0],15,world,&requests[0]);
3057 MPI_Irecv(&feq[0][0][1][0],1,passyf,comm->procneigh[1][0],25,world,&requests[1]);
3058 MPI_Isend(&feq[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],25,world,&requests[2]);
3059 MPI_Irecv(&feq[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],15,world,&requests[3]);
3060 if (typeLB == 2) {
3061 MPI_Isend(&feqn[0][1][1][0],1,passyf,comm->procneigh[1][0],10,world,&requests[4]);
3062 MPI_Irecv(&feqn[0][0][1][0],1,passyf,comm->procneigh[1][0],20,world,&requests[5]);
3063 MPI_Isend(&feqn[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],20,world,&requests[6]);
3064 MPI_Irecv(&feqn[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],10,world,&requests[7]);
3065 }
3066 update_periodic(1,2,2,subNby-2,2,subNbz-2);
3067 update_periodic(subNbx-2,subNbx-1,2,subNby-2,2,subNbz-2);
3068 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
3069
3070 for (i=0; i<numrequests; i++)
3071 requests[i]=MPI_REQUEST_NULL;
3072 MPI_Isend(&feq[0][0][1][0],1,passzf,comm->procneigh[2][0],15,world,&requests[0]);
3073 MPI_Irecv(&feq[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&requests[1]);
3074 MPI_Isend(&feq[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&requests[2]);
3075 MPI_Irecv(&feq[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&requests[3]);
3076 if (typeLB == 2) {
3077 MPI_Isend(&feqn[0][0][1][0],1,passzf,comm->procneigh[2][0],10,world,&requests[4]);
3078 MPI_Irecv(&feqn[0][0][0][0],1,passzf,comm->procneigh[2][0],20,world,&requests[5]);
3079 MPI_Isend(&feqn[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],20,world,&requests[6]);
3080 MPI_Irecv(&feqn[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],10,world,&requests[7]);
3081 }
3082 update_periodic(1,subNbx-1,1,2,2,subNbz-2);
3083 update_periodic(1,subNbx-1,subNby-2,subNby-1,2,subNbz-2);
3084 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
3085
3086 if (typeLB==1) {
3087 update_periodic(1,subNbx-1,1,subNby-1,1,2);
3088 update_periodic(1,subNbx-1,1,subNby-1,subNbz-2,subNbz-1);
3089 } else if (typeLB==2) {
3090 if (comm->myloc[2]==0) {
3091 for (i=1; i<subNbx-1; i++) {
3092 for (j=1;j<subNby-1;j++) {
3093 k=1;
3094 for (m=0; m<19; m++) {
3095 imod = i-e[m][0];
3096 jmod = j-e[m][1];
3097 kmod = k-e[m][2];
3098
3099 fnew[i][j][k][m] = feq[imod][jmod][kmod][m] + (f_lb[imod][jmod][kmod][m]-feq[imod][jmod][kmod][m])*expminusdtovertau;
3100 }
3101
3102 for (m=0; m<19; m++) {
3103 imod = i-e[m][0];
3104 jmod = j-e[m][1];
3105 kmod = k-e[m][2];
3106 imodm = i+e[m][0];
3107 jmodm = j+e[m][1];
3108 kmodm = k+e[m][2];
3109
3110 if (m==5)
3111 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][6] - feqold[imod][jmod][kmod][m]) +
3112 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][6] - feqn[imod][jmod][kmod][6]);
3113 else if (m==11)
3114 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][12] - feqold[imod][jmod][kmod][m]) +
3115 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][12] - feqn[imod][jmod][kmod][12]);
3116 else if (m==13)
3117 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][14] - feqold[imod][jmod][kmod][m]) +
3118 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][14] - feqn[imod][jmod][kmod][14]);
3119 else if (m==15)
3120 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][16] - feqold[imod][jmod][kmod][m]) +
3121 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][16] - feqn[imod][jmod][kmod][16]);
3122 else if (m==17)
3123 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][18] - feqold[imod][jmod][kmod][m]) +
3124 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][18] - feqn[imod][jmod][kmod][18]);
3125 else if (m==6)
3126 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
3127 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][5] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
3128 else if (m==12)
3129 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
3130 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][11] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
3131 else if (m==14)
3132 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
3133 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][13] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
3134 else if (m==16)
3135 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
3136 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][15] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
3137 else if (m==18)
3138 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
3139 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][17] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
3140 else
3141 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
3142 (0.5-Dcoeff*(tau+0.5))*(feqn[imodm][jmodm][kmodm][m] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
3143 }
3144 }
3145 }
3146 } else {
3147 update_periodic(1,subNbx-1,1,subNby-1,1,2);
3148 }
3149 if (comm->myloc[2]==comm->procgrid[2]-1) {
3150 for (i=1;i<subNbx-1;i++) {
3151 for (j=1;j<subNby-1;j++) {
3152 k=subNbz-2;
3153 for (m=0; m<19; m++) {
3154 imod = i-e[m][0];
3155 jmod = j-e[m][1];
3156 kmod = k-e[m][2];
3157
3158 fnew[i][j][k][m] = feq[imod][jmod][kmod][m] + (f_lb[imod][jmod][kmod][m]-feq[imod][jmod][kmod][m])*expminusdtovertau;
3159 }
3160 for (m=0; m<19; m++) {
3161 imod = i-e[m][0];
3162 jmod = j-e[m][1];
3163 kmod = k-e[m][2];
3164 imodm = i+e[m][0];
3165 jmodm = j+e[m][1];
3166 kmodm = k+e[m][2];
3167
3168 if (m==6)
3169 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][5] - feqold[imod][jmod][kmod][m]) +
3170 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][5] - feqn[imod][jmod][kmod][5]);
3171 else if (m==12)
3172 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][11] - feqold[imod][jmod][kmod][m]) +
3173 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][11] - feqn[imod][jmod][kmod][11]);
3174 else if (m==14)
3175 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][13] - feqold[imod][jmod][kmod][m]) +
3176 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][13] - feqn[imod][jmod][kmod][13]);
3177 else if (m==16)
3178 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][15] - feqold[imod][jmod][kmod][m]) +
3179 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][15] - feqn[imod][jmod][kmod][15]);
3180 else if (m==18)
3181 fnew[i][j][k][m] += Dcoeff*(feq[imod][jmod][kmod][17] - feqold[imod][jmod][kmod][m]) +
3182 (0.5-Dcoeff*(tau+0.5))*(feqoldn[imod][jmod][kmod][m] - feqoldn[imod][jmod][kmod][17] - feqn[imod][jmod][kmod][17]);
3183 else if (m==5)
3184 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
3185 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][6] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
3186 else if (m==11)
3187 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
3188 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][12] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
3189 else if (m==13)
3190 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
3191 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][14] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
3192 else if (m==15)
3193 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
3194 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][16] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
3195 else if (m==17)
3196 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
3197 (0.5-Dcoeff*(tau+0.5))*(feqn[i][j][k][18] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
3198 else
3199 fnew[i][j][k][m] += Dcoeff*(feq[i][j][k][m] - feqold[imod][jmod][kmod][m]) +
3200 (0.5-Dcoeff*(tau+0.5))*(feqn[imodm][jmodm][kmodm][m] - feqoldn[i][j][k][m] - feqn[i][j][k][m] + feqoldn[imod][jmod][kmod][m]);
3201 }
3202 }
3203 }
3204 } else {
3205 update_periodic(1,subNbx-1,1,subNby-1,subNbz-2,subNbz-1);
3206 }
3207 }
3208
3209 req_send15=MPI_REQUEST_NULL;
3210 req_recv25=MPI_REQUEST_NULL;
3211 req_send25=MPI_REQUEST_NULL;
3212 req_recv15=MPI_REQUEST_NULL;
3213
3214 if (comm->myloc[2]==0) {
3215 MPI_Isend(&fnew[0][0][1][0],1,passzf,comm->procneigh[2][0],15,world,&req_send15);
3216 MPI_Irecv(&fnew[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&req_recv25);
3217 }
3218
3219 if (comm->myloc[2]==comm->procgrid[2]-1) {
3220 MPI_Isend(&fnew[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&req_send25);
3221 MPI_Irecv(&fnew[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&req_recv15);
3222 }
3223 if (comm->myloc[2]==0) {
3224 MPI_Wait(&req_send15,&status);
3225 MPI_Wait(&req_recv25,&status);
3226
3227 for (i=1;i<subNbx-1;i++) {
3228 for (j=1;j<subNby-1;j++) {
3229 k=1;
3230
3231 if (typeLB == 1) {
3232 fnew[i][j][k][5]=fnew[i][j][k-1][6];
3233 tmp1=fnew[i][j][k-1][12]+fnew[i][j][k-1][14]+fnew[i][j][k-1][16]+fnew[i][j][k-1][18];
3234 } else {
3235 fnew[i][j][k][5]=fnew[i][j][k-1][6] + (0.5-Dcoeff*(tau+0.5))*feqn[i][j][k+1][5];
3236 tmp1=fnew[i][j][k-1][12]+fnew[i][j][k-1][14]+fnew[i][j][k-1][16]+fnew[i][j][k-1][18] +
3237 (0.5-Dcoeff*(tau+0.5))*(feqn[i-1][j][k+1][11] + feqn[i+1][j][k+1][13] +
3238 feqn[i][j-1][k+1][15] + feqn[i][j+1][k+1][17]);
3239 }
3240
3241 tmp2=fnew[i][j][k][3]+fnew[i][j][k][9]+fnew[i][j][k][10]+fnew[i][j][k][14]-
3242 fnew[i][j][k][1]-fnew[i][j][k][7]-fnew[i][j][k][8]-fnew[i][j][k][12];
3243
3244 rb=fnew[i][j][k][0]+fnew[i][j][k][1]+fnew[i][j][k][2]+fnew[i][j][k][3]+fnew[i][j][k][4]+
3245 fnew[i][j][k][5]+fnew[i][j][k][6]+fnew[i][j][k][7]+fnew[i][j][k][8]+fnew[i][j][k][9]+
3246 fnew[i][j][k][10]+fnew[i][j][k][12]+fnew[i][j][k][14]+fnew[i][j][k][16]+fnew[i][j][k][18]+tmp1;
3247
3248 tmp3=rb*vwbt-fnew[i][j][k][2]+fnew[i][j][k][4]-fnew[i][j][k][7]+fnew[i][j][k][8]-fnew[i][j][k][9]+
3249 fnew[i][j][k][10]-fnew[i][j][k][16]+fnew[i][j][k][18];
3250
3251 fnew[i][j][k][11] = 0.25*(tmp1+2.0*tmp2);
3252 fnew[i][j][k][13] = 0.25*(tmp1-2.0*tmp2);
3253 fnew[i][j][k][15] = 0.25*(tmp1+2.0*tmp3);
3254 fnew[i][j][k][17] = 0.25*(tmp1-2.0*tmp3);
3255 }
3256 }
3257
3258 }
3259 if (comm->myloc[2]==comm->procgrid[2]-1) {
3260 MPI_Wait(&req_send25,&status);
3261 MPI_Wait(&req_recv15,&status);
3262
3263 for (i=1;i<subNbx-1;i++) {
3264 for (j=1;j<subNby-1;j++) {
3265 k=subNbz-2;
3266
3267 if (typeLB == 1) {
3268 fnew[i][j][k][6]=fnew[i][j][k+1][5];
3269 tmp1=fnew[i][j][k+1][11]+fnew[i][j][k+1][13]+fnew[i][j][k+1][15]+fnew[i][j][k+1][17];
3270 } else {
3271 fnew[i][j][k][6]=fnew[i][j][k+1][5] + (0.5-Dcoeff*(tau+0.5))*feqn[i][j][k-1][6];
3272 tmp1=fnew[i][j][k+1][11]+fnew[i][j][k+1][13]+fnew[i][j][k+1][15]+fnew[i][j][k+1][17] +
3273 (0.5-Dcoeff*(tau+0.5))*(feqn[i-1][j][k-1][12] + feqn[i+1][j][k-1][14] +
3274 feqn[i][j-1][k-1][16] + feqn[i][j+1][k-1][18]);
3275 }
3276
3277 tmp2=fnew[i][j][k][3]+fnew[i][j][k][9]+fnew[i][j][k][10]+fnew[i][j][k][13]-fnew[i][j][k][1]-
3278 fnew[i][j][k][7]-fnew[i][j][k][8]-fnew[i][j][k][11];
3279
3280 rb=fnew[i][j][k][0]+fnew[i][j][k][1]+fnew[i][j][k][2]+fnew[i][j][k][3]+fnew[i][j][k][4]+
3281 fnew[i][j][k][5]+fnew[i][j][k][6]+fnew[i][j][k][7]+fnew[i][j][k][8]+fnew[i][j][k][9]+
3282 fnew[i][j][k][10]+fnew[i][j][k][11]+fnew[i][j][k][13]+fnew[i][j][k][15]+fnew[i][j][k][17]+tmp1;
3283
3284 tmp3=rb*vwtp-fnew[i][j][k][2]+fnew[i][j][k][4]-fnew[i][j][k][7]+fnew[i][j][k][8]-fnew[i][j][k][9]+
3285 fnew[i][j][k][10]-fnew[i][j][k][15]+fnew[i][j][k][17];
3286
3287 fnew[i][j][k][12] = 0.25*(tmp1+2.0*tmp2);
3288 fnew[i][j][k][14] = 0.25*(tmp1-2.0*tmp2);
3289 fnew[i][j][k][16] = 0.25*(tmp1+2.0*tmp3);
3290 fnew[i][j][k][18] = 0.25*(tmp1-2.0*tmp3);
3291 }
3292 }
3293 }
3294
3295 //--------------------------------------------------------------------------
3296 // Periodic z boundary conditions.
3297 //--------------------------------------------------------------------------
3298 } else {
3299
3300 for (i=0; i<numrequests; i++)
3301 requests[i]=MPI_REQUEST_NULL;
3302 MPI_Isend(&feq[1][1][1][0],1,passxf,comm->procneigh[0][0],15,world,&requests[0]);
3303 MPI_Irecv(&feq[0][1][1][0],1,passxf,comm->procneigh[0][0],25,world,&requests[1]);
3304 MPI_Isend(&feq[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],25,world,&requests[2]);
3305 MPI_Irecv(&feq[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],15,world,&requests[3]);
3306 if (typeLB == 2) {
3307 MPI_Isend(&feqn[1][1][1][0],1,passxf,comm->procneigh[0][0],10,world,&requests[4]);
3308 MPI_Irecv(&feqn[0][1][1][0],1,passxf,comm->procneigh[0][0],20,world,&requests[5]);
3309 MPI_Isend(&feqn[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],20,world,&requests[6]);
3310 MPI_Irecv(&feqn[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],10,world,&requests[7]);
3311 }
3312 update_periodic(2,subNbx-2,2,subNby-2,2,subNbz-2);
3313 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
3314
3315 for (i=0; i<numrequests; i++)
3316 requests[i]=MPI_REQUEST_NULL;
3317 MPI_Isend(&feq[0][1][1][0],1,passyf,comm->procneigh[1][0],15,world,&requests[0]);
3318 MPI_Irecv(&feq[0][0][1][0],1,passyf,comm->procneigh[1][0],25,world,&requests[1]);
3319 MPI_Isend(&feq[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],25,world,&requests[2]);
3320 MPI_Irecv(&feq[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],15,world,&requests[3]);
3321 if (typeLB == 2) {
3322 MPI_Isend(&feqn[0][1][1][0],1,passyf,comm->procneigh[1][0],10,world,&requests[4]);
3323 MPI_Irecv(&feqn[0][0][1][0],1,passyf,comm->procneigh[1][0],20,world,&requests[5]);
3324 MPI_Isend(&feqn[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],20,world,&requests[6]);
3325 MPI_Irecv(&feqn[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],10,world,&requests[7]);
3326 }
3327 update_periodic(1,2,2,subNby-2,2,subNbz-2);
3328 update_periodic(subNbx-2,subNbx-1,2,subNby-2,2,subNbz-2);
3329 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
3330
3331 for (i=0; i<numrequests; i++)
3332 requests[i]=MPI_REQUEST_NULL;
3333 MPI_Isend(&feq[0][0][1][0],1,passzf,comm->procneigh[2][0],15,world,&requests[0]);
3334 MPI_Irecv(&feq[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&requests[1]);
3335 MPI_Isend(&feq[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&requests[2]);
3336 MPI_Irecv(&feq[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&requests[3]);
3337 if (typeLB == 2) {
3338 MPI_Isend(&feqn[0][0][1][0],1,passzf,comm->procneigh[2][0],10,world,&requests[4]);
3339 MPI_Irecv(&feqn[0][0][0][0],1,passzf,comm->procneigh[2][0],20,world,&requests[5]);
3340 MPI_Isend(&feqn[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],20,world,&requests[6]);
3341 MPI_Irecv(&feqn[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],10,world,&requests[7]);
3342 }
3343 update_periodic(1,subNbx-1,1,2,2,subNbz-2);
3344 update_periodic(1,subNbx-1,subNby-2,subNby-1,2,subNbz-2);
3345 MPI_Waitall(numrequests,requests,MPI_STATUS_IGNORE);
3346
3347 update_periodic(1,subNbx-1,1,subNby-1,1,2);
3348 update_periodic(1,subNbx-1,1,subNby-1,subNbz-2,subNbz-1);
3349 }
3350
3351 }
3352
3353
3354