Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

This commit was manufactured by cvs2svn to create branch 'lb_experime…

…ntal'.

Sprout from master 2010-09-23 02:00:31 UTC CVS Automatic 'RELEASE_NOTES updated by cronjob'
Cherrypick from lb_merge_and_boundaries 2010-09-29 09:21:39 UTC Stefan Kesselheim <Stefan.Kesselheim@icp.uni-stuttgart.de> 'some changes :-)':
    lb-boundaries.h
    thermostat.c
    initialize.c
    lb.c
    communication.h
    lattice.c
    statistics_fluid.c
    lb.h
    communication.c
    halo.c
    lb-boundaries.c
    config.c
    random.h
    lb-d3q19.h
    statistics_fluid.h
    integrate.c
    halo.h
    lattice.h
  • Loading branch information...
commit 335896787a9df83bc8f97ebd8e5d6a9001ec1d71 1 parent 888685d
cvs2svn authored
View
139 communication.c
@@ -32,6 +32,7 @@
#include "random.h"
#include "lj.h"
#include "lb.h"
+#include "lb-boundaries.h"
#include "morse.h"
#include "buckingham.h"
#include "tab.h"
@@ -181,11 +182,13 @@ typedef void (SlaveCallback)(int node, int param);
#define REQ_ICCP3M_INIT 53
/** Action number for \ref mpi_send_rotational_inertia. */
#define REQ_SET_RINERTIA 54
-/** Action number for \ref mpi_send_mu_E. */
-#define REQ_SET_MU_E 55
+/** Action number for \ref mpi_bcast_constraint */
+#define REQ_BCAST_LBBOUNDARY 55
+/** Action number for \ref mpi_recv_fluid_border_flag */
+#define REQ_LB_GET_BORDER_FLAG 56
/** Total number of action numbers. */
-#define REQ_MAXIMUM 56
+#define REQ_MAXIMUM 57
/*@}*/
@@ -203,7 +206,6 @@ void mpi_place_particle_slave(int node, int parm);
void mpi_send_v_slave(int node, int parm);
void mpi_send_f_slave(int node, int parm);
void mpi_send_q_slave(int node, int parm);
-void mpi_send_mu_E_slave(int node, int parm);
void mpi_send_type_slave(int node, int parm);
void mpi_send_bond_slave(int node, int parm);
void mpi_recv_part_slave(int node, int parm);
@@ -219,6 +221,7 @@ void mpi_bcast_coulomb_params_slave(int node, int parm);
void mpi_send_ext_slave(int node, int parm);
void mpi_remove_particle_slave(int node, int parm);
void mpi_bcast_constraint_slave(int node, int parm);
+void mpi_bcast_lb_boundary_slave(int node, int parm);
void mpi_random_seed_slave(int node, int parm);
void mpi_random_stat_slave(int node, int parm);
void mpi_lj_cap_forces_slave(int node, int parm);
@@ -245,6 +248,7 @@ void mpi_send_dip_slave(int node, int parm);
void mpi_send_dipm_slave(int node, int parm);
void mpi_send_fluid_slave(int node, int parm);
void mpi_recv_fluid_slave(int node, int parm);
+void mpi_recv_fluid_border_flag_slave(int node, int parm);
void mpi_local_stress_tensor_slave(int node, int parm);
void mpi_ljangle_cap_forces_slave(int node, int parm);
void mpi_send_isVirtual_slave(int node, int parm);
@@ -310,7 +314,8 @@ static SlaveCallback *slave_callbacks[] = {
mpi_iccp3m_iteration_slave, /* 52: REQ_ICCP3M_ITERATION */
mpi_iccp3m_init_slave, /* 53: REQ_ICCP3M_INIT */
mpi_send_rotational_inertia_slave,/* 54: REQ_SET_RINERTIA */
- mpi_send_mu_E_slave, /* 55: REQ_SET_MU_E */
+ mpi_bcast_lb_boundary_slave, /* 55: REQ_BCAST_LBBOUNDARY */
+ mpi_recv_fluid_border_flag_slave /* 56: REQ_LB_GET_BORDER_FLAG */
};
/** Names to be printed when communication debugging is on. */
@@ -380,6 +385,8 @@ char *names[] = {
"REQ_ICCP3M_ITERATION", /* 52 */
"REQ_ICCP3M_INIT", /* 53 */
"SET_RINERTIA", /* 54 */
+ "REQ_BCAST_LBBOUNDARY", /* 55 */
+ "REQ_LB_GET_BORDER_FLAG" /* 56 */
};
/** the requests are compiled here. So after a crash you get the last issued request */
@@ -731,39 +738,6 @@ void mpi_send_q_slave(int pnode, int part)
#endif
}
-/********************* REQ_SET_MU_E ********/
-void mpi_send_mu_E(int pnode, int part, double mu_E[3])
-{
-#ifdef LB_ELECTROHYDRODYNAMICS
- mpi_issue(REQ_SET_MU_E, pnode, part);
-
- if (pnode == this_node) {
- Particle *p = local_particles[part];
- p->p.mu_E[0] = mu_E[0];
- p->p.mu_E[1] = mu_E[1];
- p->p.mu_E[2] = mu_E[2];
- }
- else {
- MPI_Send(&mu_E, 3, MPI_DOUBLE, pnode, REQ_SET_MU_E, MPI_COMM_WORLD);
- }
-
- on_particle_change();
-#endif
-}
-
-void mpi_send_mu_E_slave(int pnode, int part)
-{
-#ifdef LB_ELECTROHYDRODYNAMICS
- if (pnode == this_node) {
- Particle *p = local_particles[part];
- MPI_Status status;
- MPI_Recv(&p->p.mu_E, 3, MPI_DOUBLE, 0, REQ_SET_MU_E,
- MPI_COMM_WORLD, &status);
- }
-
- on_particle_change();
-#endif
-}
/********************* REQ_SET_M ********/
void mpi_send_mass(int pnode, int part, double mass)
@@ -1599,10 +1573,10 @@ void mpi_local_stress_tensor(DoubleList *TensorInBin, int bins[3], int periodic[
}
void mpi_local_stress_tensor_slave(int ana_num, int job) {
- int bins[3] = {0,0,0};
- int periodic[3]= {0,0,0};
- double range_start[3]= {0,0,0};
- double range[3]= {0,0,0};
+ int bins[3];
+ int periodic[3];
+ double range_start[3];
+ double range[3];
DoubleList *TensorInBin;
int i, j;
@@ -2027,6 +2001,54 @@ void mpi_bcast_constraint_slave(int node, int parm)
#endif
}
+/*************** REQ_BCAST_LBBOUNDARY ************/
+void mpi_bcast_lb_boundary(int del_num)
+{
+#ifdef LB_BOUNDARIES
+ mpi_issue(REQ_BCAST_LBBOUNDARY, 0, del_num);
+
+ if (del_num == -1) {
+ /* bcast new boundaries */
+ MPI_Bcast(&lb_boundaries[n_lb_boundaries-1], sizeof(LB_Boundary), MPI_BYTE, 0, MPI_COMM_WORLD);
+ }
+ else if (del_num == -2) {
+ /* delete all boundaries */
+ n_lb_boundaries = 0;
+ lb_boundaries = realloc(lb_boundaries,n_lb_boundaries*sizeof(LB_Boundary));
+ }
+ else {
+ memcpy(&lb_boundaries[del_num],&lb_boundaries[n_lb_boundaries-1],sizeof(LB_Boundary));
+ n_lb_boundaries--;
+ lb_boundaries = realloc(lb_boundaries,n_lb_boundaries*sizeof(LB_Boundary));
+ }
+
+ on_lb_boundary_change();
+#endif
+}
+
+void mpi_bcast_lb_boundary_slave(int node, int parm)
+{
+#ifdef LB_BOUNDARIES
+ if(parm == -1) {
+ n_lb_boundaries++;
+ lb_boundaries = realloc(lb_boundaries,n_lb_boundaries*sizeof(LB_Boundary));
+ MPI_Bcast(&lb_boundaries[n_lb_boundaries-1], sizeof(LB_Boundary), MPI_BYTE, 0, MPI_COMM_WORLD);
+ }
+ else if (parm == -2) {
+ /* delete all boundaries */
+ n_lb_boundaries = 0;
+ lb_boundaries = realloc(lb_boundaries,n_lb_boundaries*sizeof(LB_Boundary));
+ }
+ else {
+ memcpy(&lb_boundaries[parm],&lb_boundaries[n_lb_boundaries-1],sizeof(LB_Boundary));
+ n_lb_boundaries--;
+ lb_boundaries = realloc(lb_boundaries,n_lb_boundaries*sizeof(LB_Boundary));
+ }
+
+ on_lb_boundary_change();
+#endif
+}
+
/*************** REQ_RANDOM_SEED ************/
void mpi_random_seed(int cnt, long *seed) {
long this_idum = print_random_seed();
@@ -2520,7 +2542,7 @@ void mpi_send_exclusion_slave(int part1, int part2)
void mpi_send_fluid(int node, int index, double rho, double *j, double *pi) {
#ifdef LB
if (node==this_node) {
- lb_set_local_fields(index, rho, j, pi);
+ lb_calc_n_equilibrium(index, rho, j, pi);
} else {
double data[10] = { rho, j[0], j[1], j[2], pi[0], pi[1], pi[2], pi[3], pi[4], pi[5] };
mpi_issue(REQ_SET_FLUID, node, index);
@@ -2535,7 +2557,7 @@ void mpi_send_fluid_slave(int node, int index) {
double data[10];
MPI_Status status;
MPI_Recv(data, 10, MPI_DOUBLE, 0, REQ_SET_FLUID, MPI_COMM_WORLD, &status);
- lb_set_local_fields(index, data[0], &data[1], &data[4]);
+ lb_calc_n_equilibrium(index, data[0], &data[1], &data[4]);
}
#endif
}
@@ -2544,7 +2566,7 @@ void mpi_send_fluid_slave(int node, int index) {
void mpi_recv_fluid(int node, int index, double *rho, double *j, double *pi) {
#ifdef LB
if (node==this_node) {
- lb_get_local_fields(index, rho, j, pi);
+ lb_calc_local_fields(index, rho, j, pi);
} else {
double data[10];
mpi_issue(REQ_GET_FLUID, node, index);
@@ -2569,12 +2591,37 @@ void mpi_recv_fluid_slave(int node, int index) {
#ifdef LB
if (node==this_node) {
double data[10];
- lb_get_local_fields(index, &data[0], &data[1], &data[4]);
+ lb_calc_local_fields(index, &data[0], &data[1], &data[4]);
MPI_Send(data, 10, MPI_DOUBLE, 0, REQ_GET_FLUID, MPI_COMM_WORLD);
}
#endif
}
+/************** REQ_LB_GET_BORDER_FLAG **************/
+void mpi_recv_fluid_border_flag(int node, int index, int *border) {
+#ifdef LB
+ if (node==this_node) {
+ lb_local_fields_get_border_flag(index, border);
+ } else {
+ int data;
+ mpi_issue(REQ_LB_GET_BORDER_FLAG, node, index);
+ MPI_Status status;
+ MPI_Recv(&data, 1, MPI_INTEGER, node, REQ_LB_GET_BORDER_FLAG, MPI_COMM_WORLD, &status);
+ *border = data;
+ }
+#endif
+}
+
+void mpi_recv_fluid_border_flag_slave(int node, int index) {
+#ifdef LB
+ if (node==this_node) {
+ double data;
+ lb_local_fields_get_border_flag(index, &data);
+ MPI_Send(&data, 1, MPI_INTEGER, 0, REQ_LB_GET_BORDER_FLAG, MPI_COMM_WORLD);
+ }
+#endif
+}
+
/********************* REQ_ICCP3M_ITERATION ********/
int mpi_iccp3m_iteration(int dummy)
{
View
20 communication.h
@@ -142,14 +142,6 @@ void mpi_send_mass(int node, int part, double mass);
*/
void mpi_send_q(int node, int part, double q);
-/** Issue REQ_SET_MU_E: send particle electrophoretic mobility.
- Also calls \ref on_particle_change.
- \param part the particle.
- \param node the node it is attached to.
- \param mu_E its new mobility.
-*/
-void mpi_send_mu_E(int node, int part, double mu_E[3]);
-
#ifdef ROTATIONAL_INERTIA
/** Issue REQ_SET_ROTATIONAL_INERTIA: send particle rotational inertia.
Also calls \ref on_particle_change.
@@ -367,6 +359,11 @@ void mpi_send_ext(int pnode, int part, int flag, int mask, double force[3]);
/** Issue REQ_BCAST_COULOMB: send new coulomb parameters. */
void mpi_bcast_constraint(int del_num);
+#ifdef LB_BOUNDARIES
+/** Issue REQ_LB_BOUNDARY: set up walls for lb fluid */
+void mpi_bcast_lb_boundary(int del_num);
+#endif
+
/** Issue REQ_RANDOM_SEED: read/set seed of random number generators on each node. */
void mpi_random_seed(int cnt, long *seed);
@@ -437,6 +434,13 @@ void mpi_send_fluid(int node, int index, double rho, double *j, double *pi);
*/
void mpi_recv_fluid(int node, int index, double *rho, double *j, double *pi);
+/** Issue REQ_GET_FLUID: Receive a single lattice site from a processor.
+ * @param node processor to send to
+ * @param index index of the lattice site
+ * @param border local border flag
+ */
+void mpi_recv_fluid_border_flag(int node, int index, int *boundary);
+
/** Issue REQ_ICCP3M_ITERATION: performs iccp3m iteration.
@return nonzero on error
*/
View
3  config.c
@@ -33,9 +33,6 @@
#error MODES requires the fftw
#endif
-#ifdef LB
-#error LB requires the fftw
-#endif
#endif
View
137 halo.c
@@ -7,7 +7,7 @@
* You should have received a copy of that license along with this program;
* if not, refer to http://www.espresso.mpg.de/license.html where its current version can be found, or
* write to Max-Planck-Institute for Polymer Research, Theory Group, PO Box 3148, 55021 Mainz, Germany.
- * Copyright (c) 2002-2009; all rights reserved unless otherwise stated.
+ * Copyright (c) 2002-2006; all rights reserved unless otherwise stated.
*/
/** \file halo.c
@@ -24,6 +24,9 @@
#ifdef LATTICE
+/** Primitive fieldtypes and their initializers */
+struct _Fieldtype fieldtype_double = { 0, NULL, NULL, sizeof(double), 0, 0, 0, 0, NULL };
+
/** Creates a fieldtype describing the data layout
* @param count number of subtypes (Input)
* @param lengths array of lenghts of the subtytpes (Input)
@@ -36,6 +39,9 @@ void halo_create_fieldtype(int count, int* lengths, int *disps, int extent, Fiel
Fieldtype ntype = *newtype = malloc(sizeof(*ntype));
+ ntype->subtype = NULL;
+ ntype->vflag = 0;
+
ntype->vblocks = 1;
ntype->vstride = 1;
ntype->vskip = 1;
@@ -69,17 +75,45 @@ void halo_create_field_vector(int vblocks, int vstride, int vskip, Fieldtype old
int i;
Fieldtype ntype = *newtype = malloc(sizeof(*ntype));
+
+ ntype->subtype = oldtype;
+ ntype->vflag = 1;
ntype->vblocks = vblocks;
ntype->vstride = vstride;
- ntype->vskip = vskip;
-
- ntype->extent = oldtype->extent;
+ ntype->vskip = vskip;
+ ntype->extent = oldtype->extent * ((vblocks-1)*vskip + vstride);
+
int count = ntype->count = oldtype->count;
ntype->lengths = malloc(count*2*sizeof(int));
ntype->disps = (int *)((char *)ntype->lengths + count*sizeof(int));
+
+ for (i=0;i<count;i++) {
+ ntype->disps[i] = oldtype->disps[i];
+ ntype->lengths[i] = oldtype->lengths[i];
+ }
+
+}
+
+void halo_create_field_hvector(int vblocks, int vstride, int vskip, Fieldtype oldtype, Fieldtype *newtype) {
+ int i;
+ Fieldtype ntype = *newtype = malloc(sizeof(*ntype));
+
+ ntype->subtype = oldtype;
+ ntype->vflag = 0;
+
+ ntype->vblocks = vblocks;
+ ntype->vstride = vstride;
+ ntype->vskip = vskip;
+
+ ntype->extent = oldtype->extent*vstride + (vblocks-1)*vskip;
+
+ int count = ntype->count = oldtype->count;
+ ntype->lengths = malloc(count*2*sizeof(int));
+ ntype->disps = (int *)((char *)ntype->lengths + count*sizeof(int));
+
for (i=0;i<count;i++) {
ntype->disps[i] = oldtype->disps[i];
ntype->lengths[i] = oldtype->lengths[i];
@@ -125,31 +159,58 @@ MDINLINE void halo_dtset(void *dest, int value, Fieldtype type) {
}
+MDINLINE void halo_dtcopy(void *r_buffer, void *s_buffer, int count, Fieldtype type);
+
+MDINLINE void halo_copy_vector(void *r_buffer, void *s_buffer, int count, Fieldtype type, int vflag) {
+ int i, j;
+ void *dest, *src;
+
+ int vblocks = type->vblocks;
+ int vstride = type->vstride;
+ int vskip = type->vskip;
+ int extent = type->extent;
+
+ HALO_TRACE(fprintf(stderr, "%d: halo_copy_vector %p %p vblocks=%d vstride=%d vskip=%d extent=%d subtype_extent=%d\n",this_node,r_buffer,s_buffer,vblocks,vstride,vskip,extent,type->subtype->extent));
+
+ if (vflag){
+ vskip *= type->subtype->extent;
+ }
+
+ for (i=0; i<count; i++, s_buffer+=extent, r_buffer+=extent) {
+ for (j=0, dest=r_buffer, src=s_buffer; j<vblocks; j++, dest+=vskip, src+=vskip) {
+ halo_dtcopy(dest,src,vstride,type->subtype);
+ }
+ }
+
+}
+
/** Copy lattice data with layout described by fieldtype.
* @param r_buffer data destination
* @param s_buffer data source
* @param type field layout type
*/
-MDINLINE void halo_dtcopy(void *r_buffer, void *s_buffer, Fieldtype type) {
- int i, j, k;
- void *dest, *src;
-
- int vblocks = type->vblocks;
- int vstride = type->vstride;
- int vskip = type->vskip;
- int count = type->count;
- int *lens = type->lengths;
- int *disps = type->disps;
- int extent = type->extent;
-
- HALO_TRACE(fprintf(stderr, "%d: halo comm local copy r_buffer=%p s_buffer=%p\n",this_node,r_buffer,s_buffer));
-
- for (i=0; i<vblocks; i++, r_buffer+=vskip*extent, s_buffer+=vskip*extent) {
- for (j=0, dest=r_buffer, src=s_buffer; j<vstride; j++, dest+=extent, src+=extent) {
- for (k=0; k<count; k++) {
- memcpy(dest+disps[k],src+disps[k],lens[k]);
- }
+MDINLINE void halo_dtcopy(void *r_buffer, void *s_buffer, int count, Fieldtype type) {
+ int i, j;
+
+ HALO_TRACE(fprintf(stderr, "%d: halo_dtcopy r_buffer=%p s_buffer=%p blocks=%d stride=%d skip=%d\n",this_node,r_buffer,s_buffer,type->vblocks,type->vstride,type->vskip));
+
+ if (type->subtype) {
+ halo_copy_vector(r_buffer, s_buffer, count, type, type->vflag);
+ } else {
+
+ for (i=0; i<count; i++, s_buffer+=type->extent, r_buffer+=type->extent) {
+ if (!type->count) {
+ memcpy(r_buffer,s_buffer,type->extent);
+ } else {
+
+ for (j=0; j<type->count; j++) {
+ memcpy(r_buffer+type->disps[j],s_buffer+type->disps[j],type->lengths[j]);
+ }
+
}
+
+ }
+
}
}
@@ -165,13 +226,12 @@ void prepare_halo_communication(HaloCommunicator *hc, Lattice *lattice, Fieldtyp
int k, n, dir, lr, cnt, num = 0 ;
int *grid = lattice->grid ;
int *period = lattice->halo_grid ;
- void *data = lattice->data;
for (n=0; n<hc->num; n++) {
MPI_Type_free(&(hc->halo_info[n].datatype));
}
- num = 2*3; // two communications in each space direction
+ num = 2*3; /* two communications in each space direction */
hc->num = num ;
hc->halo_info = realloc(hc->halo_info,num*sizeof(HaloInfo)) ;
@@ -199,12 +259,13 @@ void prepare_halo_communication(HaloCommunicator *hc, Lattice *lattice, Fieldtyp
if (lr==0) {
/* send to left, recv from right */
- hinfo->send_buffer = &(((char *)data)[extent * stride * 1]);
- hinfo->recv_buffer = &(((char *)data)[extent * stride * (grid[dir]+1)]);
+ hinfo->s_offset = extent * stride * 1;
+ hinfo->r_offset = extent * stride * (grid[dir]+1);
} else {
/* send to right, recv from left */
- hinfo->send_buffer = &(((char *)data)[extent * stride * grid[dir]]);
- hinfo->recv_buffer = &(((char *)data)[extent * stride * 0]);
+ hinfo->s_offset = extent * stride * grid[dir];
+ hinfo->r_offset = extent * stride * 0;
+
}
hinfo->source_node = node_neighbors[2*dir+1-lr];
@@ -214,7 +275,7 @@ void prepare_halo_communication(HaloCommunicator *hc, Lattice *lattice, Fieldtyp
MPI_Type_vector(nblocks, stride, skip, datatype, &hinfo->datatype);
MPI_Type_commit(&hinfo->datatype);
-
+
#ifdef PARTIAL_PERIODIC
if ( !PERIODIC(dir) && (boundary[2*dir+lr] != 0 || boundary[2*dir+1-lr] != 0) ) {
if (node_grid[dir] == 1) {
@@ -245,7 +306,7 @@ void prepare_halo_communication(HaloCommunicator *hc, Lattice *lattice, Fieldtyp
}
}
- HALO_TRACE(fprintf(stderr,"%d: prepare_halo_communication dir=%d lr=%d s_buffer=%p r_buffer=%p, s_node=%d d_node=%d extent=%d type=%d\n",this_node,dir,lr,hinfo->send_buffer,hinfo->recv_buffer,hinfo->source_node,hinfo->dest_node,(int)extent,hinfo->type)) ;
+ HALO_TRACE(fprintf(stderr,"%d: prepare_halo_communication dir=%d lr=%d s_offset=%ld r_offset=%ld s_node=%d d_node=%d type=%d\n",this_node,dir,lr,hinfo->s_offset,hinfo->r_offset,hinfo->source_node,hinfo->dest_node,hinfo->type)) ;
cnt++;
@@ -261,7 +322,6 @@ void release_halo_communication(HaloCommunicator *hc) {
int n;
for (n=0; n<hc->num; n++) {
- HALO_TRACE(fprintf(stderr,"%d: freeing %p\n",this_node,&(hc->halo_info[n].datatype)));
MPI_Type_free(&(hc->halo_info[n].datatype));
}
@@ -273,7 +333,7 @@ void release_halo_communication(HaloCommunicator *hc) {
* described by the halo communicator
* @param hc halo communicator describing the parallelization scheme
*/
-void halo_communication(HaloCommunicator *hc) {
+void halo_communication(HaloCommunicator *hc, void *base) {
int n, comm_type, s_node, r_node;
void *s_buffer, *r_buffer ;
@@ -282,21 +342,21 @@ void halo_communication(HaloCommunicator *hc) {
MPI_Request request;
MPI_Status status;
- HALO_TRACE(fprintf(stderr, "%d: halo_comm %p (num=%d)\n", this_node, hc, hc->num)) ;
+ HALO_TRACE(fprintf(stderr, "%d: halo_comm base=%p num=%d\n", this_node, base, hc->num)) ;
for (n = 0; n < hc->num; n++) {
HALO_TRACE(fprintf(stderr, "%d: halo_comm round %d\n", this_node, n)) ;
comm_type = hc->halo_info[n].type ;
- s_buffer = hc->halo_info[n].send_buffer ;
- r_buffer = hc->halo_info[n].recv_buffer ;
+ s_buffer = (char *)base + hc->halo_info[n].s_offset;
+ r_buffer = (char *)base + hc->halo_info[n].r_offset;
switch (comm_type) {
case HALO_LOCL:
fieldtype = hc->halo_info[n].fieldtype;
- halo_dtcopy(r_buffer,s_buffer,fieldtype);
+ halo_dtcopy(r_buffer,s_buffer,1,fieldtype);
break ;
case HALO_SENDRECV:
@@ -340,6 +400,7 @@ void halo_communication(HaloCommunicator *hc) {
HALO_TRACE(fprintf(stderr,"%d: halo_comm open boundaries\n",this_node));
+ /* \todo this does not work for the n_i - <n_i> */
halo_dtset(r_buffer,0,fieldtype);
break;
@@ -347,8 +408,6 @@ void halo_communication(HaloCommunicator *hc) {
}
- HALO_TRACE(fprintf(stderr, "%d: halo_comm %p finished\n", this_node, hc));
-
}
#endif /* LATTICE */
View
33 halo.h
@@ -7,7 +7,7 @@
* You should have received a copy of that license along with this program;
* if not, refer to http://www.espresso.mpg.de/license.html where its current version can be found, or
* write to Max-Planck-Institute for Polymer Research, Theory Group, PO Box 3148, 55021 Mainz, Germany.
- * Copyright (c) 2002-2009; all rights reserved unless otherwise stated.
+ * Copyright (c) 2002-2006; all rights reserved unless otherwise stated.
*/
/** \file halo.h
@@ -57,7 +57,8 @@
* is similar to MPI datatypes but a bit more compact. See \ref
* halo_create_fieldtype, \ref halo_create_field_vector and \ref
* halo_dtcopy to understand how it works. */
-typedef struct _Fieldtype {
+typedef struct _Fieldtype *Fieldtype;
+struct _Fieldtype {
int count; /**< number of subtypes in fieldtype */
int *disps; /**< displacements of the subtypes */
int *lengths; /**< lengths of the subtypes */
@@ -65,21 +66,30 @@ typedef struct _Fieldtype {
int vblocks; /**< number of blocks in field vectors */
int vstride; /**< size of strides in field vectors */
int vskip; /**< displacement between strides in field vectors */
-} *Fieldtype;
+ int vflag;
+ Fieldtype subtype;
+};
+
+/** Predefined fieldtypes */
+extern struct _Fieldtype fieldtype_double;
+#define FIELDTYPE_DOUBLE (&fieldtype_double)
/** Structure describing a Halo region */
typedef struct {
- int type; /**< type of halo communication */
+ int type; /**< type of halo communication */
+
+ int source_node; /**< index of processor which sends halo data */
+ int dest_node; /**< index of processor receiving halo data */
- int source_node; /**< index of processor which sends halo data */
- int dest_node; /**< index of processor receiving halo data */
+ //void *send_buffer; /**< pointer to data being sent */
+ //void *recv_buffer; /**< pointer to data being received */
- void *send_buffer; /**< pointer to data being sent */
- void *recv_buffer; /**< pointer to data being received */
+ unsigned long s_offset; /**< offset for send buffer */
+ unsigned long r_offset; /**< offset for receive buffer */
- Fieldtype fieldtype; /**< type layout of the data beeing exchanged */
- MPI_Datatype datatype; /**< MPI datatype of data beeing communicated */
+ Fieldtype fieldtype; /**< type layout of the data beeing exchanged */
+ MPI_Datatype datatype; /**< MPI datatype of data beeing communicated */
} HaloInfo ;
@@ -110,6 +120,7 @@ void halo_create_fieldtype(int count, int *lens, int *disps, int extent, Fieldty
* @param newtype newly created fieldtype (Input/Output)
*/
void halo_create_field_vector(int vblocks, int vstride, int vskip, Fieldtype oldtype, Fieldtype *newtype);
+void halo_create_field_hvector(int vblocks, int vstride, int vskip, Fieldtype oldtype, Fieldtype *newtype);
/** Frees a fieldtype
* @param ftype pointer to the type to be freed (Input)
@@ -134,7 +145,7 @@ void release_halo_communication(HaloCommunicator *hc);
* described by the halo communicator
* @param hc halo communicator describing the parallelization scheme
*/
-void halo_communication(HaloCommunicator *hc);
+void halo_communication(HaloCommunicator *hc, void *base);
#endif /* LATTICE */
View
38 initialize.c
@@ -93,6 +93,11 @@ int on_program_start(Tcl_Interp *interp)
fft_pre_init();
#endif
+#ifdef LB
+ lb_pre_init();
+#endif
+
+
/*
call all initializations to do only on the master node here.
*/
@@ -170,19 +175,19 @@ void on_integration_start()
#ifdef LB
if(lattice_switch & LATTICE_LB) {
- if (lbpar.agrid < 0.0) {
+ if (lbpar.agrid <= 0.0) {
errtext = runtime_error(128);
ERROR_SPRINTF(errtext,"{098 Lattice Boltzmann agrid not set} ");
}
- if (lbpar.tau < 0.0) {
+ if (lbpar.tau <= 0.0) {
errtext = runtime_error(128);
ERROR_SPRINTF(errtext,"{099 Lattice Boltzmann time step not set} ");
}
- if (lbpar.rho < 0.0) {
+ if (lbpar.rho <= 0.0) {
errtext = runtime_error(128);
ERROR_SPRINTF(errtext,"{100 Lattice Boltzmann fluid density not set} ");
}
- if (lbpar.viscosity < 0.0) {
+ if (lbpar.viscosity <= 0.0) {
errtext = runtime_error(128);
ERROR_SPRINTF(errtext,"{101 Lattice Boltzmann fluid viscosity not set} ");
}
@@ -358,12 +363,26 @@ void on_constraint_change()
EVENT_TRACE(fprintf(stderr, "%d: on_constraint_change\n", this_node));
invalidate_obs();
-#ifdef LB
-#ifdef CONSTRAINTS
+#ifdef LB_BOUNDARIES
if(lattice_switch & LATTICE_LB) {
- lb_init_constraints();
+ lb_init_boundaries();
}
#endif
+
+ recalc_forces = 1;
+}
+
+void on_lb_boundary_change()
+{
+ EVENT_TRACE(fprintf(stderr, "%d: on_lb_boundary_change\n", this_node));
+ invalidate_obs();
+
+#ifdef LB_BOUNDARIES
+ //printf("executing on_lb_boundary_change on node %d\n", this_node);
+
+ if(lattice_switch & LATTICE_LB) {
+ lb_init_boundaries();
+ }
#endif
recalc_forces = 1;
@@ -696,6 +715,10 @@ static void init_tcl(Tcl_Interp *interp)
REGISTER_COMMAND("blockfile", blockfile);
/* in constraint.c */
REGISTER_COMMAND("constraint", constraint);
+ /* in lb-boundaries.c */
+#ifdef LB_BOUNDARIES
+ REGISTER_COMMAND("lb_boundary", lb_boundary);
+#endif
/* in uwerr.c */
REGISTER_COMMAND("uwerr", uwerr);
/* in nemd.c */
@@ -706,6 +729,7 @@ static void init_tcl(Tcl_Interp *interp)
REGISTER_COMMAND("bin", bin);
/* in lb.c */
REGISTER_COMMAND("lbfluid", lbfluid_cmd);
+ REGISTER_COMMAND("lbnode", lbnode_cmd);
/* in utils.h */
REGISTER_COMMAND("replacestdchannel", replacestdchannel);
/* in iccp3m.h */
View
2  integrate.c
@@ -519,7 +519,7 @@ void integrate_vv(int n_steps)
rescale_forces_propagate_vel();
#ifdef LB
- if (lattice_switch & LATTICE_LB) lb_propagate();
+ if (lattice_switch & LATTICE_LB) lattice_boltzmann_update();
if (check_runtime_errors()) break;
#endif
View
10 lattice.c
@@ -7,7 +7,7 @@
* You should have received a copy of that license along with this program;
* if not, refer to http://www.espresso.mpg.de/license.html where its current version can be found, or
* write to Max-Planck-Institute for Polymer Research, Theory Group, PO Box 3148, 55021 Mainz, Germany.
- * Copyright (c) 2002-2009; all rights reserved unless otherwise stated.
+ * Copyright (c) 2002-2006; all rights reserved unless otherwise stated.
*/
/** \file lattice.c
@@ -20,6 +20,8 @@
#include "grid.h"
#include "lattice.h"
+#ifdef LATTICE
+
/** Switch determining the type of lattice dynamics. A value of zero
* means that there is no lattice dynamics. Different types can be
* combined by or'ing the respective flags.
@@ -27,8 +29,6 @@
*/
int lattice_switch = LATTICE_OFF ;
-#ifdef LATTICE
-
/** Initialize lattice.
*
* This function initializes the variables describing the lattice
@@ -51,8 +51,8 @@ void init_lattice(Lattice *lattice, double agrid, double tau) {
for (dir=0;dir<3;dir++) {
/* check if local_box_l is compatible with lattice spacing */
if (fabs(local_box_l[dir]-lattice->grid[dir]*agrid) > ROUND_ERROR_PREC) {
- char *errtxt = runtime_error(128 + 4*TCL_DOUBLE_SPACE + 4*TCL_INTEGER_SPACE);
- ERROR_SPRINTF(errtxt, "{097 Lattice spacing agrid=%g is incompatible with local_box_l[%d]=%g (box_l[%d]=%g node_grid[%d]=%d) %g} ",agrid,dir,local_box_l[dir],dir,box_l[dir],dir,node_grid[dir],local_box_l[dir]-lattice->grid[dir]*agrid);
+ char *errtxt = runtime_error(128);
+ ERROR_SPRINTF(errtxt, "{097 Lattice spacing agrid=%f is incompatible with local_box_l[%d]=%f (box_l[%d]=%f node_grid[%d]=%d) %f} ",agrid,dir,local_box_l[dir],dir,box_l[dir],dir,node_grid[dir],local_box_l[dir]-lattice->grid[dir]*agrid);
return;
}
}
View
58 lattice.h
@@ -7,7 +7,7 @@
* You should have received a copy of that license along with this program;
* if not, refer to http://www.espresso.mpg.de/license.html where its current version can be found, or
* write to Max-Planck-Institute for Polymer Research, Theory Group, PO Box 3148, 55021 Mainz, Germany.
- * Copyright (c) 2002-2009; all rights reserved unless otherwise stated.
+ * Copyright (c) 2002-2006; all rights reserved unless otherwise stated.
*/
/** \file lattice.h
@@ -24,6 +24,12 @@
#include "grid.h"
#include "particle_data.h"
+extern int lattice_switch;
+
+#ifdef LATTICE
+
+#define index_t long
+
/** Lattice off */
#define LATTICE_OFF 0
@@ -35,9 +41,6 @@
* combined by or'ing the respective flags.
* So far, only \ref LATTICE_OFF and \ref LATTICE_LB exist.
*/
-extern int lattice_switch ;
-
-#ifdef LATTICE
/** Data structure describing a lattice.
* Contains the lattice layout and pointers to the data fields.
@@ -52,14 +55,14 @@ typedef struct _Lattice {
int halo_grid[3] ;
/** total number (volume) of local lattice sites (excluding halo) */
- int grid_volume ;
+ index_t grid_volume;
/** total number (volume) of lattice sites (including halo) */
- int halo_grid_volume ;
+ index_t halo_grid_volume;
/** number of lattice sites in the halo region */
- int halo_grid_surface ;
+ index_t halo_grid_surface;
/** offset for number of halo sites stored in front of the local
lattice sites */
- int halo_offset ;
+ index_t halo_offset;
/** lattice constant */
double agrid ;
@@ -79,8 +82,9 @@ typedef struct _Lattice {
* Actually used are only the identity and the type. */
Particle part_rep;
- /** MPI datatype that describes the data layout of the lattice. */
+ /** datatypes that describe the data layout of the lattice. */
MPI_Datatype datatype;
+ struct _Fieldtype *fieldtype;
} Lattice;
@@ -103,7 +107,7 @@ void init_lattice(Lattice *lattice, double agrid, double tau);
*
* \param lattice pointer to the lattice
* \param ind global coordinates of the lattice site (Input)
- * \param ind local coordinates of the lattice site (Output)
+ * \param grid local coordinates of the lattice site (Output)
* \return index of the node for the lattice site
*/
MDINLINE int map_lattice_to_node(Lattice *lattice, int *ind, int *grid) {
@@ -115,15 +119,31 @@ MDINLINE int map_lattice_to_node(Lattice *lattice, int *ind, int *grid) {
//fprintf(stderr,"%d: (%d,%d,%d)\n",this_node,grid[0],grid[1],grid[2]);
- /* change from global to local lattice coordinates */
- ind[0] -= grid[0]*lattice->grid[0];
- ind[1] -= grid[1]*lattice->grid[1];
- ind[2] -= grid[2]*lattice->grid[2];
+ /* change from global to local lattice coordinates (+1 is for halo) */
+ ind[0] = ind[0] - grid[0]*lattice->grid[0] + 1;
+ ind[1] = ind[1] - grid[1]*lattice->grid[1] + 1;
+ ind[2] = ind[2] - grid[2]*lattice->grid[2] + 1;
/* return linear index into node array */
return map_array_node(grid);
}
+/** Map a local lattice site to the global position.
+ *
+ * This function determines the processor responsible for
+ * the specified lattice site. The coordinates of the site are
+ * taken as global coordinates andcoordinates of the site are
+ * taken as global coordinates and are returned as local coordinates.
+ *
+ * \param lattice pointer to the lattice
+ * \param ind global coordinates of the lattice site (Input)
+ * \param ind local coordinates of the lattice site (Output)
+ * \return index of the node for the lattice site
+ */
+MDINLINE int map_lattice_to_position(Lattice *lattice, int *ind, int *grid) {
+ return 0;
+}
+
/** Map a spatial position to the surrounding lattice sites.
*
* This function takes a global spatial position and determines the
@@ -140,7 +160,7 @@ MDINLINE int map_lattice_to_node(Lattice *lattice, int *ind, int *grid) {
* \param delta distance fraction of pos from the surrounding
* elementary cell, 6 directions (Output)
*/
-MDINLINE void map_position_to_lattice(Lattice *lattice, const double pos[3], int node_index[8], double delta[6]) {
+MDINLINE void map_position_to_lattice(Lattice *lattice, const double pos[3], index_t node_index[8], double delta[6]) {
int dir,ind[3] ;
double lpos, rel;
@@ -155,13 +175,17 @@ MDINLINE void map_position_to_lattice(Lattice *lattice, const double pos[3], int
/* surrounding elementary cell is not completely inside this box,
adjust if this is due to round off errors */
if (ind[dir] < 0) {
- fprintf(stderr,"%d: map_position_to_lattice: position (%f,%f,%f) not inside a local plaquette.\n",this_node,pos[0],pos[1],pos[2]);
+ if (abs(rel) < ROUND_ERROR_PREC) {
+ ind[dir] = 0; // TODO
+ } else {
+ fprintf(stderr,"%d: map_position_to_lattice: position (%f,%f,%f) not inside a local plaquette in dir %d ind[dir]=%d rel=%f lpos=%f.\n",this_node,pos[0],pos[1],pos[2],dir,ind[dir],rel,lpos);
+ }
}
else if (ind[dir] > lattice->grid[dir]) {
if (lpos - local_box_l[dir] < ROUND_ERROR_PREC)
ind[dir] = lattice->grid[dir];
else
- fprintf(stderr,"%d: map_position_to_lattice: position (%f,%f,%f) not inside a local plaquette.\n",this_node,pos[0],pos[1],pos[2]);
+ fprintf(stderr,"%d: map_position_to_lattice: position (%f,%f,%f) not inside a local plaquette in dir %d ind[dir]=%d rel=%f lpos=%f.\n",this_node,pos[0],pos[1],pos[2],dir,ind[dir],rel,lpos);
}
delta[3+dir] = rel - ind[dir]; // delta_x/a
View
637 lb-boundaries.c
@@ -11,7 +11,7 @@
* where its current version can be found, or write to
* Max-Planck-Institute for Polymer Research, Theory Group,
* PO Box 3148, 55021 Mainz, Germany.
- * Copyright (c) 2002-2009; all rights reserved unless otherwise stated.
+ * Copyright (c) 2002-2007; all rights reserved unless otherwise stated.
*/
/** \file lb-boundaries.c
@@ -20,65 +20,648 @@
* Header file for \ref lb-boundaries.h.
*
*/
-
#include "utils.h"
#include "constraint.h"
#include "lb-boundaries.h"
+#include "lb.h"
+#include "interaction_data.h"
+
+#ifdef LB_BOUNDARIES
+
+int n_lb_boundaries = 0;
+LB_Boundary *lb_boundaries = NULL;
+
+#ifdef LB_BOUNDARIES
-#ifdef LB
-#ifdef CONSTRAINTS
+// LB_Boundary lb_boundary_par = { LB_BOUNDARY_NONE, 0.0 }; //do we need this?
/** Initialize a planar boundary specified by a wall constraint.
* @param plane The \ref Constraint_wall struct describing the boundary.
*/
-static void lb_init_constraint_wall(Constraint_wall* plane) {
- int x, y, z;
- double pos[3], dist;
+int printLbBoundaryToResult(Tcl_Interp *interp, int i)
+{
+ LB_Boundary *lbb = &lb_boundaries[i];
+ char buffer[TCL_DOUBLE_SPACE + TCL_INTEGER_SPACE];
+ sprintf(buffer, "%d ", i);
+ Tcl_AppendResult(interp, buffer, (char *)NULL);
+
+ switch (lbb->type) {
+ case LB_BOUNDARY_WAL:
+ Tcl_PrintDouble(interp, lbb->c.wal.n[0], buffer);
+ Tcl_AppendResult(interp, "wall normal ", buffer, " ", (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.wal.n[1], buffer);
+ Tcl_AppendResult(interp, buffer, " ", (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.wal.n[2], buffer);
+ Tcl_AppendResult(interp, buffer, (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.wal.d, buffer);
+ Tcl_AppendResult(interp, " dist ", buffer, (char *) NULL);
+ break;
+ case LB_BOUNDARY_SPH:
+ Tcl_PrintDouble(interp, lbb->c.sph.pos[0], buffer);
+ Tcl_AppendResult(interp, "sphere center ", buffer, " ", (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.sph.pos[1], buffer);
+ Tcl_AppendResult(interp, buffer, " ", (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.sph.pos[2], buffer);
+ Tcl_AppendResult(interp, buffer, (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.sph.rad, buffer);
+ Tcl_AppendResult(interp, " radius ", buffer, (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.sph.direction, buffer);
+ Tcl_AppendResult(interp, " direction ", buffer, (char *) NULL);
+ break;
+ case LB_BOUNDARY_CYL:
+ Tcl_PrintDouble(interp, lbb->c.cyl.pos[0], buffer);
+ Tcl_AppendResult(interp, "cylinder center ", buffer, " ", (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.cyl.pos[1], buffer);
+ Tcl_AppendResult(interp, buffer, " ", (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.cyl.pos[2], buffer);
+ Tcl_AppendResult(interp, buffer, (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.cyl.axis[0], buffer);
+ Tcl_AppendResult(interp, " axis ", buffer, " ", (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.cyl.axis[1], buffer);
+ Tcl_AppendResult(interp, buffer, " ", (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.cyl.axis[2], buffer);
+ Tcl_AppendResult(interp, buffer, (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.cyl.rad, buffer);
+ Tcl_AppendResult(interp, " radius ", buffer, (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.cyl.length, buffer);
+ Tcl_AppendResult(interp, " length ", buffer, (char *) NULL);
+ Tcl_PrintDouble(interp, lbb->c.cyl.direction, buffer);
+ Tcl_AppendResult(interp, " direction ", buffer, (char *) NULL);
+ break;
+
+ default:
+ sprintf(buffer, "%d", lbb->type);
+ Tcl_AppendResult(interp, "unknown lb_boundary type ", buffer, ".", (char *) NULL);
+ return (TCL_OK);
+ }
+
+ return (TCL_OK);
+}
+
+int lb_boundary_print_all(Tcl_Interp *interp)
+{
+ int i;
+ if(n_lb_boundaries>0) Tcl_AppendResult(interp, "{", (char *)NULL);
+ for (i = 0; i < n_lb_boundaries; i++) {
+ if(i>0) Tcl_AppendResult(interp, " {", (char *)NULL);
+ printLbBoundaryToResult(interp, i);
+ Tcl_AppendResult(interp, "}", (char *)NULL);
+ }
+ return (TCL_OK);
+}
+
+/*void printLbBoundaryForceToResult(Tcl_Interp *interp, int con)
+{
+ double f[3];
+ char buffer[TCL_DOUBLE_SPACE];
+
+ mpi_get_lb_boundary_force(lbb, f);
+
+ Tcl_PrintDouble(interp, f[0], buffer);
+ Tcl_AppendResult(interp, buffer, " ", (char *) NULL);
+ Tcl_PrintDouble(interp, f[1], buffer);
+ Tcl_AppendResult(interp, buffer, " ", (char *) NULL);
+ Tcl_PrintDouble(interp, f[2], buffer);
+ Tcl_AppendResult(interp, buffer, (char *) NULL);
+}*/
+
+LB_Boundary *generate_lb_boundary()
+{
+ n_lb_boundaries++;
+ lb_boundaries = realloc(lb_boundaries,n_lb_boundaries*sizeof(LB_Boundary));
+ lb_boundaries[n_lb_boundaries-1].type = LB_BOUNDARY_NONE;
+
+ return &lb_boundaries[n_lb_boundaries-1];
+}
+
+int lb_boundary_wall(LB_Boundary *lbb, Tcl_Interp *interp,
+ int argc, char **argv)
+{
+ int i;
+ double norm;
+ lbb->type = LB_BOUNDARY_WAL;
+ /* invalid entries to start of */
+ lbb->c.wal.n[0] =
+ lbb->c.wal.n[1] =
+ lbb->c.wal.n[2] = 0;
+ lbb->c.wal.d = 0;
+ while (argc > 0) {
+ if(!strncmp(argv[0], "normal", strlen(argv[0]))) {
+ if(argc < 4) {
+ Tcl_AppendResult(interp, "lb_boundary wall normal <nx> <ny> <nz> expected", (char *) NULL);
+ return (TCL_ERROR);
+ }
+ if (Tcl_GetDouble(interp, argv[1], &(lbb->c.wal.n[0])) == TCL_ERROR ||
+ Tcl_GetDouble(interp, argv[2], &(lbb->c.wal.n[1])) == TCL_ERROR ||
+ Tcl_GetDouble(interp, argv[3], &(lbb->c.wal.n[2])) == TCL_ERROR)
+ return (TCL_ERROR);
+ argc -= 4; argv += 4;
+ }
+ else if(!strncmp(argv[0], "dist", strlen(argv[0]))) {
+ if (argc < 1) {
+ Tcl_AppendResult(interp, "lb_boundary wall dist <d> expected", (char *) NULL);
+ return (TCL_ERROR);
+ }
+ if (Tcl_GetDouble(interp, argv[1], &(lbb->c.wal.d)) == TCL_ERROR)
+ return (TCL_ERROR);
+ argc -= 2; argv += 2;
+ }
+ else if(!strncmp(argv[0], "type", strlen(argv[0]))) {
+ if (argc < 1) {
+ Tcl_AppendResult(interp, "lb_boundary wall type <t> expected", (char *) NULL);
+ return (TCL_ERROR);
+ }
+ argc -= 2; argv += 2;
+ }
+ else
+ break;
+ }
+ /* length of the normal vector */
+ norm = SQR(lbb->c.wal.n[0])+SQR(lbb->c.wal.n[1])+SQR(lbb->c.wal.n[2]);
+ if (norm < 1e-10) {
+ Tcl_AppendResult(interp, "usage: lb_boundary wall normal <nx> <ny> <nz> dist <d> type <t>",
+ (char *) NULL);
+ return (TCL_ERROR);
+ }
+ /* normalize the normal vector */
+ for (i=0;i<3;i++) lbb->c.wal.n[i] /= sqrt(norm);
+
+ return (TCL_OK);
+}
+
+int lb_boundary_sphere(LB_Boundary *lbb, Tcl_Interp *interp,
+ int argc, char **argv)
+{
+ lbb->type = LB_BOUNDARY_SPH;
+
+ /* invalid entries to start of */
+ lbb->c.sph.pos[0] =
+ lbb->c.sph.pos[1] =
+ lbb->c.sph.pos[2] = 0;
+ lbb->c.sph.rad = 0;
+ lbb->c.sph.direction = -1;
+
+ while (argc > 0) {
+ if(!strncmp(argv[0], "center", strlen(argv[0]))) {
+ if(argc < 4) {
+ Tcl_AppendResult(interp, "lb_boundary sphere center <x> <y> <z> expected", (char *) NULL);
+ return (TCL_ERROR);
+ }
+ if (Tcl_GetDouble(interp, argv[1], &(lbb->c.sph.pos[0])) == TCL_ERROR ||
+ Tcl_GetDouble(interp, argv[2], &(lbb->c.sph.pos[1])) == TCL_ERROR ||
+ Tcl_GetDouble(interp, argv[3], &(lbb->c.sph.pos[2])) == TCL_ERROR)
+ return (TCL_ERROR);
+ argc -= 4; argv += 4;
+ }
+ else if(!strncmp(argv[0], "radius", strlen(argv[0]))) {
+ if (argc < 1) {
+ Tcl_AppendResult(interp, "lb_boundary sphere radius <r> expected", (char *) NULL);
+ return (TCL_ERROR);
+ }
+ if (Tcl_GetDouble(interp, argv[1], &(lbb->c.sph.rad)) == TCL_ERROR)
+ return (TCL_ERROR);
+ argc -= 2; argv += 2;
+ }
+ else if(!strncmp(argv[0], "direction", strlen(argv[0]))) {
+ if (argc < 1) {
+ Tcl_AppendResult(interp, "-1/1 or inside/outside is expected", (char *) NULL);
+ return (TCL_ERROR);
+ }
+ if (!strncmp(argv[1], "inside", strlen(argv[1])))
+ lbb->c.sph.direction = -1;
+ else if (!strncmp(argv[1], "outside", strlen(argv[1])))
+ lbb->c.sph.direction = 1;
+ else if (Tcl_GetDouble(interp, argv[1], &(lbb->c.sph.direction)) == TCL_ERROR)
+ return (TCL_ERROR);
+ argc -= 2; argv += 2;
+ }
+ else if(!strncmp(argv[0], "type", strlen(argv[0]))) {
+ if (argc < 1) {
+ Tcl_AppendResult(interp, "lb_boundary sphere type <t> expected", (char *) NULL);
+ return (TCL_ERROR);
+ }
+ argc -= 2; argv += 2;
+ }
+ else
+ break;
+ }
+
+ if (lbb->c.sph.rad < 0.) {
+ Tcl_AppendResult(interp, "usage: lb_boundary sphere center <x> <y> <z> radius <d> direction <direction> type <t>",
+ (char *) NULL);
+ return (TCL_ERROR);
+ }
+
+ return (TCL_OK);
+}
+
+int lb_boundary_cylinder(LB_Boundary *lbb, Tcl_Interp *interp,
+ int argc, char **argv)
+{
+ double axis_len;
+ int i;
+
+ lbb->type = LB_BOUNDARY_CYL;
+ /* invalid entries to start of */
+ lbb->c.cyl.pos[0] =
+ lbb->c.cyl.pos[1] =
+ lbb->c.cyl.pos[2] = 0;
+ lbb->c.cyl.axis[0] =
+ lbb->c.cyl.axis[1] =
+ lbb->c.cyl.axis[2] = 0;
+ lbb->c.cyl.rad = 0;
+ lbb->c.cyl.length = 0;
+ lbb->c.cyl.direction = 0;
+ while (argc > 0) {
+ if(!strncmp(argv[0], "center", strlen(argv[0]))) {
+ if(argc < 4) {
+ Tcl_AppendResult(interp, "lb_boundary cylinder center <x> <y> <z> expected", (char *) NULL);
+ return (TCL_ERROR);
+ }
+ if (Tcl_GetDouble(interp, argv[1], &(lbb->c.cyl.pos[0])) == TCL_ERROR ||
+ Tcl_GetDouble(interp, argv[2], &(lbb->c.cyl.pos[1])) == TCL_ERROR ||
+ Tcl_GetDouble(interp, argv[3], &(lbb->c.cyl.pos[2])) == TCL_ERROR)
+ return (TCL_ERROR);
+ argc -= 4; argv += 4;
+ }
+ else if(!strncmp(argv[0], "axis", strlen(argv[0]))) {
+ if(argc < 4) {
+ Tcl_AppendResult(interp, "lb_boundary cylinder axis <rx> <ry> <rz> expected", (char *) NULL);
+ return (TCL_ERROR);
+ }
+ if (Tcl_GetDouble(interp, argv[1], &(lbb->c.cyl.axis[0])) == TCL_ERROR ||
+ Tcl_GetDouble(interp, argv[2], &(lbb->c.cyl.axis[1])) == TCL_ERROR ||
+ Tcl_GetDouble(interp, argv[3], &(lbb->c.cyl.axis[2])) == TCL_ERROR)
+ return (TCL_ERROR);
+
+ argc -= 4; argv += 4;
+ }
+ else if(!strncmp(argv[0], "radius", strlen(argv[0]))) {
+ if (argc < 1) {
+ Tcl_AppendResult(interp, "lb_boundary cylinder radius <rad> expected", (char *) NULL);
+ return (TCL_ERROR);
+ }
+ if (Tcl_GetDouble(interp, argv[1], &(lbb->c.cyl.rad)) == TCL_ERROR)
+ return (TCL_ERROR);
+ argc -= 2; argv += 2;
+ }
+ else if(!strncmp(argv[0], "length", strlen(argv[0]))) {
+ if (argc < 1) {
+ Tcl_AppendResult(interp, "lb_boundary cylinder length <len> expected", (char *) NULL);
+ return (TCL_ERROR);
+ }
+ if (Tcl_GetDouble(interp, argv[1], &(lbb->c.cyl.length)) == TCL_ERROR)
+ return (TCL_ERROR);
+ argc -= 2; argv += 2;
+ }
+ else if(!strncmp(argv[0], "direction", strlen(argv[0]))) {
+ if (argc < 1) {
+ Tcl_AppendResult(interp, "lb_boundary cylinder direction <dir> expected", (char *) NULL);
+ return (TCL_ERROR);
+ }
+ if (!strncmp(argv[1], "inside", strlen(argv[1])))
+ lbb->c.cyl.direction = -1;
+ else if (!strncmp(argv[1], "outside", strlen(argv[1])))
+ lbb->c.cyl.direction = 1;
+ else if (Tcl_GetDouble(interp, argv[1], &(lbb->c.cyl.direction)) == TCL_ERROR)
+ return (TCL_ERROR);
+ argc -= 2; argv += 2;
+ }
+ else if(!strncmp(argv[0], "type", strlen(argv[0]))) {
+ if (argc < 1) {
+ Tcl_AppendResult(interp, "lb_boundary cylinder type <t> expected", (char *) NULL);
+ return (TCL_ERROR);
+ }
+ argc -= 2; argv += 2;
+ }
+ else
+ break;
+ }
+
+ axis_len=0.;
+ for (i=0;i<3;i++)
+ axis_len += SQR(lbb->c.cyl.axis[i]);
+
+ if (lbb->c.cyl.rad < 0. || axis_len < 1e-30 ||
+ lbb->c.cyl.direction == 0 || lbb->c.cyl.length <= 0) {
+ Tcl_AppendResult(interp, "usage: lb_boundary cylinder center <x> <y> <z> axis <rx> <ry> <rz> radius <rad> length <length> direction <direction> type <t>",
+ (char *) NULL);
+ return (TCL_ERROR);
+ }
- for (x=0;x<lblattice.halo_grid[0];x++) {
- for (y=0;y<lblattice.halo_grid[1];y++) {
- for (z=0;z<lblattice.halo_grid[2];z++) {
+ /*normalize the axis vector */
+ axis_len = sqrt (axis_len);
+ for (i=0;i<3;i++) {
+ lbb->c.cyl.axis[i] /= axis_len;
+ }
+
+ return (TCL_OK);
+}
- pos[0] = my_left[0] + (x-1)*lblattice.agrid;
- pos[1] = my_left[1] + (y-1)*lblattice.agrid;
- pos[2] = my_left[2] + (z-1)*lblattice.agrid;
+#endif /* LB_BOUNDARIES */
- dist = scalar(pos,plane->n) - plane->d;
+int lb_boundary(ClientData _data, Tcl_Interp *interp,
+ int argc, char **argv)
+{
+#ifdef LB_BOUNDARIES
+ int status, c_num;
- if (fabs(dist) < lblattice.agrid) {
- //printf("%d %d %d\n",x,y,z);
- lbfluid[get_linear_index(x,y,z,lblattice.halo_grid)].boundary = 1;
- }
+ if (argc < 2) return lb_boundary_print_all(interp);
+
+ if(!strncmp(argv[1], "wall", strlen(argv[1]))) {
+ status = lb_boundary_wall(generate_lb_boundary(),interp, argc - 2, argv + 2);
+ mpi_bcast_lb_boundary(-1);
+ }
+ else if(!strncmp(argv[1], "sphere", strlen(argv[1]))) {
+ status = lb_boundary_sphere(generate_lb_boundary(),interp, argc - 2, argv + 2);
+ mpi_bcast_lb_boundary(-1);
+ }
+ else if(!strncmp(argv[1], "cylinder", strlen(argv[1]))) {
+ status = lb_boundary_cylinder(generate_lb_boundary(),interp, argc - 2, argv + 2);
+ mpi_bcast_lb_boundary(-1);
+ }
+ else if(!strncmp(argv[1], "delete", strlen(argv[1]))) {
+ if(argc < 3) {
+ /* delete all */
+ mpi_bcast_lb_boundary(-2);
+ status = TCL_OK;
+ }
+ else {
+ if(Tcl_GetInt(interp, argv[2], &(c_num)) == TCL_ERROR) return (TCL_ERROR);
+ if(c_num < 0 || c_num >= n_lb_boundaries) {
+ Tcl_AppendResult(interp, "Can not delete non existing lb_boundary",(char *) NULL);
+ return (TCL_ERROR);
}
+ mpi_bcast_lb_boundary(c_num);
+ status = TCL_OK;
}
}
+ else if (argc == 2 && Tcl_GetInt(interp, argv[1], &c_num) == TCL_OK) {
+ printLbBoundaryToResult(interp, c_num);
+ status = TCL_OK;
+ }
+ else {
+ Tcl_AppendResult(interp, "possible lb_boundaries: wall sphere cylinder lb_boundary delete {c} to delete lb_boundary or lb_boundaries",(char *) NULL);
+ return (TCL_ERROR);
+ }
+
+ return mpi_gather_runtime_errors(interp, status);
+#else /* !defined(LB_BOUNDARIES) */
+ Tcl_AppendResult(interp, "LB_BOUNDARIES not compiled in!" ,(char *) NULL);
+ return (TCL_ERROR);
+#endif /* LB_BOUNDARIES */
+}
+
+#ifdef LB_BOUNDARIES
+
+static void lb_init_boundary_wall(Constraint_wall* wall) {
+ int x, y, z, index, node_domain_position[3], offset[3];
+ double pos[3], dist, dist_vec[3];
+
+ map_node_array(this_node, node_domain_position); //constraint contains MD coordinates, therefore the conversion
+
+ offset[0] = node_domain_position[0]*lblattice.grid[0];
+ offset[1] = node_domain_position[1]*lblattice.grid[1];
+ offset[2] = node_domain_position[2]*lblattice.grid[2];
+
+ for (z=1; z<=lblattice.grid[2]; z++) {
+ for (y=1; y<=lblattice.grid[1]; y++) {
+ for (x=1; x<=lblattice.grid[0]; x++) {
+
+ pos[0] = (offset[0]+(x-1))*lblattice.agrid;
+ pos[1] = (offset[1]+(y-1))*lblattice.agrid;
+ pos[2] = (offset[2]+(z-1))*lblattice.agrid;
+
+ calculate_wall_dist((Particle*) NULL, pos, (Particle*) NULL, wall, &dist, dist_vec);
+
+ if (dist <= 0) {
+ lbfields[get_linear_index(x,y,z,lblattice.halo_grid)].boundary = 1;
+ }
+ }
+ }
+ }
+}
+
+static void lb_init_boundary_sphere(Constraint_sphere* sphere) {
+ int x, y, z, index, node_domain_position[3], offset[3];
+ double pos[3], dist, dist_vec[3];
+
+ map_node_array(this_node, node_domain_position); //constraint contains MD coordinates, therefore the conversion
+
+ offset[0] = node_domain_position[0]*lblattice.grid[0];
+ offset[1] = node_domain_position[1]*lblattice.grid[1];
+ offset[2] = node_domain_position[2]*lblattice.grid[2];
+
+ for (z=1; z<=lblattice.grid[2]; z++) {
+ for (y=1; y<=lblattice.grid[1]; y++) {
+ for (x=1; x<=lblattice.grid[0]; x++) {
+ pos[0] = (offset[0]+(x-1))*lblattice.agrid;
+ pos[1] = (offset[1]+(y-1))*lblattice.agrid;
+ pos[2] = (offset[2]+(z-1))*lblattice.agrid;
+
+ calculate_sphere_dist((Particle*) NULL, pos, (Particle*) NULL, sphere, &dist, dist_vec);
+
+ if (dist <= 0)
+ lbfields[get_linear_index(x,y,z,lblattice.halo_grid)].boundary = 1;
+ }
+ }
+ }
+}
+
+static void lb_init_boundary_cylinder(Constraint_cylinder* cylinder) {
+ int x, y, z, index, node_domain_position[3], offset[3];
+ double pos[3], dist, dist_vec[3];
+
+ map_node_array(this_node, node_domain_position);
+
+ offset[0] = node_domain_position[0]*lblattice.grid[0];
+ offset[1] = node_domain_position[1]*lblattice.grid[1];
+ offset[2] = node_domain_position[2]*lblattice.grid[2];
+
+ for (z=1; z<=lblattice.grid[2]; z++) {
+ for (y=1; y<=lblattice.grid[1]; y++) {
+ for (x=1; x<=lblattice.grid[0]; x++) {
+ pos[0] = (offset[0]+(x-1))*lblattice.agrid;
+ pos[1] = (offset[1]+(y-1))*lblattice.agrid;
+ pos[2] = (offset[2]+(z-1))*lblattice.agrid;
+
+ calculate_cylinder_dist((Particle*) NULL, pos, (Particle*) NULL, cylinder, &dist, dist_vec);
+
+ if (dist <= 0) {
+ lbfields[get_linear_index(x,y,z,lblattice.halo_grid)].boundary = 1;
+ }
+ }
+ }
+ }
}
/** Initialize boundary conditions for all constraints in the system. */
-void lb_init_constraints() {
+void lb_init_boundaries() {
+ int n, x, y, z, index, node_domain_position[3], offset[3];
+ char *errtxt;
+ double pos[3], dist, dist_tmp, dist_vec[3];
+
+ map_node_array(this_node, node_domain_position);
+
+ offset[0] = node_domain_position[0]*lblattice.grid[0];
+ offset[1] = node_domain_position[1]*lblattice.grid[1];
+ offset[2] = node_domain_position[2]*lblattice.grid[2];
+
+ for (n=0;n<lblattice.halo_grid_volume;n++) {
+ lbfields[n].boundary = 0;
+ }
+
+ for (z=1; z<=lblattice.grid[2]; z++) {
+ for (y=1; y<=lblattice.grid[1]; y++) {
+ for (x=1; x<=lblattice.grid[0]; x++) {
+ pos[0] = (offset[0]+(x-1))*lblattice.agrid;
+ pos[1] = (offset[1]+(y-1))*lblattice.agrid;
+ pos[2] = (offset[2]+(z-1))*lblattice.agrid;
+
+ dist = 0.;
+
+ for (n=0;n<n_lb_boundaries;n++) {
+ switch (lb_boundaries[n].type) {
+ case LB_BOUNDARY_WAL:
+ calculate_wall_dist((Particle*) NULL, pos, (Particle*) NULL, &lb_boundaries[n].c.wal, &dist_tmp, dist_vec);
+ break;
+ case LB_BOUNDARY_SPH:
+ calculate_sphere_dist((Particle*) NULL, pos, (Particle*) NULL, &lb_boundaries[n].c.sph, &dist_tmp, dist_vec);
+ break;
+ case LB_BOUNDARY_CYL:
+ calculate_cylinder_dist((Particle*) NULL, pos, (Particle*) NULL, &lb_boundaries[n].c.cyl, &dist_tmp, dist_vec);
+ break;
+ default:
+ errtxt = runtime_error(128);
+ ERROR_SPRINTF(errtxt, "{109 lb_boundary type %d not implemented in lb_init_boundaries()\n", lb_boundaries[n].type);
+ }
+
+ if (abs(dist) > abs(dist_tmp) || n == 0) {
+ dist = dist_tmp;
+ }
+ }
+
+ if (dist <= 0 && n_lb_boundaries > 0) {
+ lbfields[get_linear_index(x,y,z,lblattice.halo_grid)].boundary = 1;
+ }
+ }
+ }
+ }
+}
+/** Initialize boundary conditions for all constraints in the system. */
+/*void lb_init_boundaries() {
int n;
char *errtxt;
+
+ //printf("executing lb_init_boundaries on node %d\n", this_node);
for (n=0;n<lblattice.halo_grid_volume;n++) {
- lbfluid[n].boundary = 0;
+ lbfields[n].boundary = 0;
}
- for (n=0;n<n_constraints;n++) {
- switch (constraints[n].type) {
- case CONSTRAINT_WAL:
- lb_init_constraint_wall(&constraints[n].c.wal);
+ for (n=0;n<n_lb_boundaries;n++) {
+ switch (lb_boundaries[n].type) {
+ case LB_BOUNDARY_WAL:
+ lb_init_boundary_wall(&lb_boundaries[n].c.wal);
+ break;
+ case LB_BOUNDARY_SPH:
+ lb_init_boundary_sphere(&lb_boundaries[n].c.sph);
+ break;
+ case LB_BOUNDARY_CYL:
+ lb_init_boundary_cylinder(&lb_boundaries[n].c.cyl);
break;
default:
errtxt = runtime_error(128);
- ERROR_SPRINTF(errtxt, "{109 constraint type %d not implemented in lb_init_constraints()\n",constraints[n].type);
+ ERROR_SPRINTF(errtxt, "{109 lb_boundary type %d not implemented in lb_init_boundaries()\n",lb_boundaries[n].type);
}
}
+}*/
+
+
+static int lbboundaries_parse_slip_reflection(Tcl_Interp *interp, int argc, char **argv) {
+#if 0 //problems with slip_pref (georg, 03.08.10)
+ if (argc <1) {
+ Tcl_AppendResult(interp, "lbboundaries slip_reflection requires 1 argument", (char *)NULL);
+ return TCL_ERROR;
+ }
+ if (!ARG0_IS_D(lb_boundary_par.slip_pref)) {
+ Tcl_AppendResult(interp, "wrong argument for lbboundaries slip_reflection", (char *)NULL);
+ return TCL_ERROR;
+ }
+
+ lb_boundary_par.type = LB_BOUNDARY_SLIP_REFLECTION;
+
+// mpi_bcast_lb_params(LBPAR_BOUNDARY);
+
+ return TCL_OK;
+#endif //if 0
+}
+
+static int lbboundaries_parse_partial_slip(Tcl_Interp *interp, int argc, char **argv) {
+#if 0 //problems with slip_pref (georg, 03.08.10)
+ if (argc <1) {
+ Tcl_AppendResult(interp, "lbboundaries slip requires 1 argument", (char *)NULL);
+ return TCL_ERROR;
+ }
+ if (!ARG0_IS_D(lb_boundary_par.slip_pref)) {
+ Tcl_AppendResult(interp, "wrong argument for lbboundaries slip", (char *)NULL);
+ return TCL_ERROR;
+ }
+
+ lb_boundary_par.type = LB_BOUNDARY_PARTIAL_SLIP;
+
+// mpi_bcast_lb_params(LBPAR_BOUNDARY);
+
+ return TCL_OK;
+#endif //if 0
+}
+
+#endif /* LB_BOUNDARIES */
+
+/** Parser for the \ref lbfluid command. */
+int lbboundaries_cmd(ClientData data, Tcl_Interp *interp, int argc, char **argv) {
+#if 0 //problems with slip_pref (georg, 03.08.10)
+#ifdef LB_BOUNDARIES
+ int err = TCL_ERROR;
+
+ if (argc < 2) {
+ Tcl_AppendResult(interp, "too few arguments to \"lbboundaries\"", (char *)NULL);
+ err = TCL_ERROR;
+ }
+ else if (ARG1_IS_S("off")) {
+ err = TCL_ERROR;
+ }
+ else if (ARG1_IS_S("bounce_back")) {
+ lb_boundary_par.type = LB_BOUNDARY_BOUNCE_BACK;
+ err = TCL_OK;
+ }
+ else if (ARG1_IS_S("specular_reflections")) {
+ lb_boundary_par.type = LB_BOUNDARY_SPECULAR_REFLECTION;
+ err = TCL_OK;
+ }
+ else if (ARG1_IS_S("slip_reflection")) {
+ err = lbboundaries_parse_slip_reflection(interp, argc-2, argv+2);
+ }
+ else if (ARG1_IS_S("partial_slip")) {
+ err = lbboundaries_parse_partial_slip(interp, argc-2, argv+2);
+ }
+ else {
+ Tcl_AppendResult(interp, "unkown boundary condition \"", argv[1], (char *)NULL);
+ err = TCL_ERROR;
+ }
+
+ //mpi_bcast_lb_boundaries();
+ return err;
+#else /* !defined LB_BOUNDARIES */
+ Tcl_AppendResult(interp, "LB_BOUNDARIES not compiled in!", NULL);
+ return TCL_ERROR;
+#endif /* LB_BOUNDARIES */
+#endif //if 0
}
+#endif /* LB_BOUNDARIES */
-#endif /* CONSTRAINTS */
-#endif /* LB */
View
3,751 lb-boundaries.h
3,667 additions, 84 deletions not shown
View
27 lb-d3q19.h
@@ -74,7 +74,32 @@ static double d3q19_w[19] = { 1./3.,
1./36., 1./36., 1./36., 1./36.,
1./36., 1./36., 1./36., 1./36. };
-LB_Model d3q19_model = { 19, d3q19_lattice, d3q19_coefficients, d3q19_w, 1./3. };
+/** Basis of the mode space as described in [Duenweg, Schiller, Ladd] */
+static double d3q19_modebase[20][19] = {
+ { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 },
+ { 0.0, 1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 0.0 },
+ { 0.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, 1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, -1.0, 1.0, -1.0 },
+ { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0 },
+ { -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 },
+ { 0.0, 1.0, 1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0 },
+ { -0.0, 1.0, 1.0, 1.0, 1.0, -2.0, -2.0, 2.0, 2.0, 2.0, 2.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0 },
+ { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
+ { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 0.0 },
+ { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, -1.0, -1.0 },
+ { 0.0, -2.0, 2.0, 0.0, 0.0, 0.0, 0.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 0.0 },
+ { 0.0, 0.0, 0.0, -2.0, 2.0, 0.0, 0.0, 1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, -1.0, 1.0, -1.0 },
+ { 0.0, 0.0, 0.0, 0.0, 0.0, -2.0, 2.0, 0.0, 0.0, 0.0, 0.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0 },
+ { 0.0, -0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
+ { 0.0, 0.0, 0.0, -0.0, 0.0, 0.0, 0.0, 1.0, -1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 1.0, -1.0, 1.0 },
+ { 0.0, 0.0, 0.0, 0.0, 0.0, -0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0 },
+ { 1.0, -2.0, -2.0, -2.0, -2.0, -2.0, -2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 },
+ { 0.0, -1.0, -1.0, 1.0, 1.0, -0.0, -0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0 },
+ { 0.0, -1.0, -1.0, -1.0, -1.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0 },
+ /* the following values are the (weighted) lengths of the vectors */
+ { 1.0, 1./3., 1./3., 1./3., 2./3., 4./9., 4./3., 1./9., 1./9., 1./9., 2./3., 2./3., 2./3., 2./9., 2./9., 2./9., 2.0, 4./9., 4./3. }
+};
+
+LB_Model d3q19_model = { 19, d3q19_lattice, d3q19_coefficients, d3q19_w, NULL, 1./3. };
#endif /* LB */
View
2,702 lb.c
1,891 additions, 811 deletions not shown
View
398 lb.h
@@ -47,8 +47,26 @@
#define LBPAR_FRICTION 4 /**< friction coefficient for viscous coupling between particles and fluid */
#define LBPAR_EXTFORCE 5 /**< external force acting on the fluid */
#define LBPAR_BULKVISC 6 /**< fluid bulk viscosity */
+
/*@}*/
+ /** Some general remarks:
+ * This file implements the LB D3Q19 method to Espresso. The LB_Model
+ * construction is preserved for historical reasons and might be removed
+ * soon. It is constructed as a multi-relaxation time LB, thus all populations
+ * are converted to modes, then collision is performed and transfered back
+ * to population space, where the streaming is performed.
+ *
+ * For performance reasons it is clever to do streaming and collision at the same time
+ * because every fluid node has to be read and written only once. This increases
+ * mainly cache efficiency.
+ * Two alternatives are implemented: stream_collide and collide_stream.
+ *
+ * The hydrodynamic fields, corresponding to density, velocity and stress, are
+ * stored in LB_FluidNodes in the array lbfields, the populations in lbfluid
+ * which is constructed as 2 x (Nx x Ny x Nz) x 19 array.
+ */
+
/** Description of the LB Model in terms of the unit vectors of the
* velocity sub-lattice and the corresponding coefficients
* of the pseudo-equilibrium distribution */
@@ -58,7 +76,7 @@ typedef const struct {
int n_veloc ;
/** unit vectors of the velocity sublattice */
- double (*c)[3] ;
+ double (*c)[3];
/** coefficients in the pseudo-equilibrium distribution */
double (*coeff)[4];
@@ -66,6 +84,9 @@ typedef const struct {
/** weights in the functional for the equilibrium distribution */
double (*w);
+ /** basis of moment space */
+ double **e;
+
/** speed of sound squared */
double c_sound_sq;
@@ -74,6 +95,9 @@ typedef const struct {
/** Data structure for fluid on a local lattice site */
typedef struct {
+ /** flag indicating whether fields have to be recomputed */
+ int recalc_fields;
+
/** local density */
double rho[1];
@@ -83,19 +107,21 @@ typedef struct {
/** local stress tensor */
double pi[6];
- /** local populations of the velocity directions */
- double *n;
-#ifndef D3Q19
- /** temporary storage for populations */
- double *n_tmp;
-#endif
+ /** local populations of the velocity directions
+ * are stored seperately to achieve higher performance */
+
+ /** flag indicating whether a force is acting on the node */
+ int has_force;
-#ifdef CONSTRAINTS
+ /** local force density TODO: FORCE DENSITY or FORCE?*/
+ double force[3];
+
+#ifdef LB_BOUNDARIES
/** flag indicating whether this site belongs to a boundary */
int boundary;
/** normal vector of the boundary surface */
- double *nvec;
+ double *nvec; //doesn't work like that any more, I think (georg, 17.08.10)
#endif
} LB_FluidNode;
@@ -124,8 +150,9 @@ typedef struct {
* lead to numerical artifacts with low order integrators */
double friction;
- /** external force applied to the fluid at each lattice site (LJ units) */
- double ext_force[3];
+ /** external force applied to the fluid at each lattice site (MD units) */
+ double ext_force[3]; /* TODO: local force or global force? */
+ double rho_lb_units;
} LB_Parameters;
@@ -138,23 +165,38 @@ extern LB_Parameters lbpar;
/** The underlying lattice */
extern Lattice lblattice;
-/** Pointer to the fluid nodes */
-extern LB_FluidNode *lbfluid;
+/** Pointer to the velocity populations of the fluid.
+ * lbfluid[0] contains pre-collision populations, lbfluid[1]
+ * contains post-collision populations*/
+extern double **lbfluid[2];
+
+/** Pointer to the hydrodynamic fields of the fluid */
+extern LB_FluidNode *lbfields;
/** Switch indicating momentum exchange between particles and fluid */
extern int transfer_momentum;
/** Eigenvalue of collision operator corresponding to shear viscosity. */
-//extern double lblambda;
+extern double lblambda;
/** Eigenvalue of collision operator corresponding to bulk viscosity. */
-//extern double lblambda_bulk;
+extern double lblambda_bulk;
/************************************************************/
/** \name Exported Functions */
/************************************************************/
/*@{*/
+/** Updates the Lattice Boltzmann system for one time step.
+ * This function performs the collision step and the streaming step.
+ * If external forces are present, they are applied prior to the collisions.
+ * If boundaries are present, it also applies the boundary conditions.
+ */
+void lattice_boltzmann_update();
+
+/** (Pre-)initializes data structures. */
+void lb_pre_init();
+
/** Performs a full initialization of
* the Lattice Boltzmann system. All derived parameters
* and the fluid are reset to their default values. */
@@ -168,20 +210,27 @@ void lb_reinit_parameters();
/** (Re-)initializes the fluid. */
void lb_reinit_fluid();
+/** Resets the forces on the fluid nodes */
+void lb_reinit_forces();
+
/** Sets the density and momentum on a local lattice site.
- * @param index The index of the lattice site within the local domain (Input)
+ * @param node Pointer to the Node of the lattice site within the local domain (Input)
* @param rho Local density of the fluid (Input)
* @param j Local momentum of the fluid (Input)
*/
-void lb_set_local_fields(int index, const double rho, const double *v, const double *pi);
+void lb_set_local_fields(LB_FluidNode *node, const double rho, const double *v, const double *pi);
/** Returns the mass, momentum and stress of a local lattice site.
- * @param index The index of the lattice site within the local domain (Input)
+ * @param node The index of the lattice site within the local domain (Input)
* @param rho Local density of the fluid (Output)
* @param j Local momentum of the fluid (Output)
* @param pi Local stress tensor of the fluid (Output)
*/
-void lb_get_local_fields(int index, double *rho, double *j, double *pi);
+void lb_get_local_fields(LB_FluidNode *node, double *rho, double *j, double *pi);
+
+/** Calculates the equilibrium distributions TODO: What does it mean?*/
+void lb_calc_n_equilibrium(const index_t index, const double rho, const double *j, double *pi);
+
/** Propagates the Lattice Boltzmann system for one time step.
* This function performs the collision step and the streaming step.
@@ -202,26 +251,27 @@ void calc_particle_lattice_ia();
* The calculation is implemented explicitly for the special case of D3Q19.
* @param local_node The local lattice site (Input).
*/
-MDINLINE void lb_calc_local_rho(LB_FluidNode *local_node) {
-
- double *local_n = local_node->n;
- double *local_rho = local_node->rho;
- double avg_rho = lbpar.rho*(lbpar.agrid*lbpar.agrid*lbpar.agrid);
+MDINLINE void lb_calc_local_rho(index_t index, double *rho) {
+ // unit conversion: mass density
+ double avg_rho = lbpar.rho*lbpar.agrid*lbpar.agrid*lbpar.agrid;
#ifdef D3Q19
- *local_rho = avg_rho
- + local_n[0]
- + local_n[1] + local_n[2]
- + local_n[3] + local_n[4]
- + local_n[5] + local_n[6]
- + local_n[7] + local_n[8] + local_n[9] + local_n[10]
- + local_n[11] + local_n[12] + local_n[13] + local_n[14]
- + local_n[15] + local_n[16] + local_n[17] + local_n[18];
+ *rho = avg_rho
+ + lbfluid[0][0][index]
+ + lbfluid[0][1][index] + lbfluid[0][2][index]
+ + lbfluid[0][3][index] + lbfluid[0][4][index]
+ + lbfluid[0][5][index] + lbfluid[0][6][index]
+ + lbfluid[0][7][index] + lbfluid[0][8][index]
+ + lbfluid[0][9][index] + lbfluid[0][10][index]
+ + lbfluid[0][11][index] + lbfluid[0][12][index]
+ + lbfluid[0][13][index] + lbfluid[0][14][index]
+ + lbfluid[0][15][index] + lbfluid[0][16][index]
+ + lbfluid[0][17][index] + lbfluid[0][18][index];
#else
int i;
- *local_rho = 0.0;
+ *rho = avg_rho;
for (i=0;i<lbmodel.n_veloc;i++) {
- *local_rho += local_n[i] + lbmodel.coeff[i][0]*avg_rho;
+ *rho += lbfluid[0][i][index];// + lbmodel.coeff[i][0]*avg_rho;
}
#endif
@@ -231,82 +281,99 @@ MDINLINE void lb_calc_local_rho(LB_FluidNode *local_node) {
* The calculation is implemented explicitly for the special case of D3Q19.
* @param local_node The local lattice site (Input).
*/
-MDINLINE void lb_calc_local_j(LB_FluidNode *local_node) {
-
- double *local_n = local_node->n;
- double *local_j = local_node->j;
+MDINLINE void lb_calc_local_j(index_t index, double *j) {
#ifdef D3Q19
- local_j[0] = local_n[1] - local_n[2]
- + local_n[7] - local_n[8] + local_n[9] - local_n[10]
- + local_n[11] - local_n[12] + local_n[13] - local_n[14];
- local_j[1] = local_n[3] - local_n[4]
- + local_n[7] - local_n[8] - local_n[9] + local_n[10]
- + local_n[15] - local_n[16] + local_n[17] - local_n[18];
- local_j[2] = local_n[5] - local_n[6]
- + local_n[11] - local_n[12] - local_n[13] + local_n[14]
- + local_n[15] - local_n[16] - local_n[17] + local_n[18];
+ j[0] = lbfluid[0][1][index] - lbfluid[0][2][index]
+ + lbfluid[0][7][index] - lbfluid[0][8][index]
+ + lbfluid[0][9][index] - lbfluid[0][10][index]
+ + lbfluid[0][11][index] - lbfluid[0][12][index]
+ + lbfluid[0][13][index] - lbfluid[0][14][index];
+ j[1] = lbfluid[0][3][index] - lbfluid[0][4][index]
+ + lbfluid[0][7][index] - lbfluid[0][8][index]
+ - lbfluid[0][9][index] + lbfluid[0][10][index]
+ + lbfluid[0][15][index] - lbfluid[0][16][index]
+ + lbfluid[0][17][index] - lbfluid[0][18][index];
+ j[2] = lbfluid[0][5][index] - lbfluid[0][6][index]
+ + lbfluid[0][11][index] - lbfluid[0][12][index]
+ - lbfluid[0][13][index] + lbfluid[0][14][index]
+ + lbfluid[0][15][index] - lbfluid[0][16][index]
+ - lbfluid[0][17][index] + lbfluid[0][18][index];
#else
int i;
double tmp;
- double avg_rho = lbpar.rho*(lbpar.agrid*lbpar.agrid*lbpar.agrid);
- local_j[0] = 0.0;
- local_j[1] = 0.0;
- local_j[2] = 0.0;
+ //double avg_rho = lbpar.rho/(lbpar.agrid*lbpar.agrid*lbpar.agrid);
+ j[0] = 0.0;
+ j[1] = 0.0;
+ j[2] = 0.0;
for (i=0;i<lbmodel.n_veloc;i++) {
- tmp = local_n[i] + lbmodel.coeff[i][0]*avg_rho;
- local_j[0] += lbmodel.c[i][0] * tmp;
- local_j[1] += lbmodel.c[i][1] * tmp;
- local_j[2] += lbmodel.c[i][2] * tmp;
+ tmp = lbfluid[0][i][index];// + lbmodel.coeff[i][0]*avg_rho;
+ j[0] += lbmodel.c[i][0] * tmp;
+ j[1] += lbmodel.c[i][1] * tmp;
+ j[2] += lbmodel.c[i][2] * tmp;
}
#endif
+#ifdef EXTERNAL_FORCES
+ /* the coupling forces are not yet included self-consistently */
+ j[0] += 0.5*lbpar.ext_force[0];
+ j[1] += 0.5*lbpar.ext_force[1];
+ j[2] += 0.5*lbpar.ext_force[2];
+#endif
+
}
/** Calculate the local fluid stress.
* The calculation is implemented explicitly for the special case of D3Q19.
* @param local_node The local lattice site (Input).
*/
-MDINLINE void lb_calc_local_pi(LB_FluidNode *local_node) {
+MDINLINE void lb_calc_local_pi(index_t index, double *pi) {
- double *local_n = local_node->n;
- double *local_pi = local_node->pi;
- double avg_rho = lbpar.rho*(lbpar.agrid*lbpar.agrid*lbpar.agrid);
+ double avg_rho = lbpar.rho*lbpar.agrid*lbpar.agrid*lbpar.agrid;
#ifdef D3Q19
- local_pi[0] = avg_rho/3.0
- + local_n[1] + local_n[2]
- + local_n[7] + local_n[8] + local_n[9] + local_n[10]
- + local_n[11] + local_n[12] + local_n[13] + local_n[14];
- local_pi[2] = avg_rho/3.0
- + local_n[3] + local_n[4]
- + local_n[7] + local_n[8] + local_n[9] + local_n[10]
- + local_n[15] + local_n[16] + local_n[17] + local_n[18];
- local_pi[5] = avg_rho/3.0
- + local_n[5] + local_n[6]
- + local_n[11] + local_n[12] + local_n[13] + local_n[14]
- + local_n[15] + local_n[16] + local_n[17] + local_n[18];
- local_pi[1] = local_n[7] + local_n[8] - local_n[9] - local_n[10];
- local_pi[3] = local_n[11] + local_n[12] - local_n[13] - local_n[14];
- local_pi[4] = local_n[15] + local_n[16] - local_n[17] - local_n[18];
+ pi[0] = avg_rho/3.0
+ + lbfluid[0][1][index] + lbfluid[0][2][index]
+ + lbfluid[0][7][index] + lbfluid[0][8][index]
+ + lbfluid[0][9][index] + lbfluid[0][10][index]
+ + lbfluid[0][11][index] + lbfluid[0][12][index]
+ + lbfluid[0][13][index] + lbfluid[0][14][index];
+ pi[2] = avg_rho/3.0
+ + lbfluid[0][3][index] + lbfluid[0][4][index]
+ + lbfluid[0][7][index] + lbfluid[0][8][index]
+ + lbfluid[0][9][index] + lbfluid[0][10][index]
+ + lbfluid[0][15][index] + lbfluid[0][16][index]
+ + lbfluid[0][17][index] + lbfluid[0][18][index];
+ pi[5] = avg_rho/3.0
+ + lbfluid[0][5][index] + lbfluid[0][6][index]
+ + lbfluid[0][11][index] + lbfluid[0][12][index]
+ + lbfluid[0][13][index] + lbfluid[0][14][index]
+ + lbfluid[0][15][index] + lbfluid[0][16][index]
+ + lbfluid[0][17][index] + lbfluid[0][18][index];
+ pi[1] = lbfluid[0][7][index] + lbfluid[0][8][index]
+ - lbfluid[0][9][index] - lbfluid[0][10][index];
+ pi[3] = lbfluid[0][11][index] + lbfluid[0][12][index]
+ - lbfluid[0][13][index] - lbfluid[0][14][index];
+ pi[4] = lbfluid[0][15][index] + lbfluid[0][16][index]
+ - lbfluid[0][17][index] - lbfluid[0][18][index];
#else
int i;
double tmp;
double (*c)[3] = lbmodel.c;
- local_pi[0] = 0.0;
- local_pi[1] = 0.0;
- local_pi[2] = 0.0;
- local_pi[3] = 0.0;
- local_pi[4] = 0.0;
- local_pi[5] = 0.0;
+ pi[0] = 0.0;
+ pi[1] = 0.0;
+ pi[2] = 0.0;
+ pi[3] = 0.0;
+ pi[4] = 0.0;
+ pi[5] = 0.0;
for (i=0;i<lbmodel.n_veloc;i++) {
- tmp = local_n[i] + lbmodel.coeff[i][0]*avg_rho;
- local_pi[0] += c[i][0] * c[i][0] * tmp;
- local_pi[1] += c[i][0] * c[i][1] * tmp;
- local_pi[2] += c[i][1] * c[i][1] * tmp;
- local_pi[3] += c[i][0] * c[i][2] * tmp;
- local_pi[4] += c[i][1] * c[i][2] * tmp;
- local_pi[5] += c[i][2] * c[i][2] * tmp;
+ tmp = lbfluid[0][i][index] + lbmodel.coeff[i][0]*avg_rho;
+ pi[0] += c[i][0] * c[i][0] * tmp;
+ pi[1] += c[i][0] * c[i][1] * tmp;
+ pi[2] += c[i][1] * c[i][1] * tmp;
+ pi[3] += c[i][0] * c[i][2] * tmp;
+ pi[4] += c[i][1] * c[i][2] * tmp;
+ pi[5] += c[i][2] * c[i][2] * tmp;
}
#endif
@@ -321,93 +388,126 @@ MDINLINE void lb_calc_local_pi(LB_FluidNode *local_node) {
* @param calc_pi_flag Flag indicating whether stress tensor should be
* computed.
*/
-MDINLINE void lb_calc_local_fields(LB_FluidNode *local_node,int calc_pi_flag) {
+MDINLINE void lb_calc_local_fields(index_t index, double *rho, double *j, double *pi) {
- double *local_n = local_node->n;
- double *local_rho = local_node->rho;
- double *local_j = local_node->j;
- double *local_pi = local_node->pi;
- double avg_rho = lbpar.rho*(lbpar.agrid*lbpar.agrid*lbpar.agrid);
+ double avg_rho = lbpar.rho*lbpar.agrid*lbpar.agrid*lbpar.agrid;
#ifdef D3Q19
- *local_rho = avg_rho
- + local_n[0]
- + local_n[1] + local_n[2]
- + local_n[3] + local_n[4]
- + local_n[5] + local_n[6]
- + local_n[7] + local_n[8] + local_n[9] + local_n[10]
- + local_n[11] + local_n[12] + local_n[13] + local_n[14]
- + local_n[15] + local_n[16] + local_n[17] + local_n[18];
-
- local_j[0] = local_n[1] - local_n[2]
- + local_n[7] - local_n[8] + local_n[9] - local_n[10]
- + local_n[11] - local_n[12] + local_n[13] - local_n[14];
- local_j[1] = local_n[3] - local_n[4]
- + local_n[7] - local_n[8] - local_n[9] + local_n[10]
- + local_n[15] - local_n[16] + local_n[17] - local_n[18];
- local_j[2] = local_n[5] - local_n[6]
- + local_n[11] - local_n[12] - local_n[13] + local_n[14]
- + local_n[15] - local_n[16] - local_n[17] + local_n[18];
+#ifdef LB_BOUNDARIES
+ if ( lbfields[index].boundary ) {
+ *rho = avg_rho;
+ j[0] = 0.; j[1] = 0.; j[2] = 0.;
+ pi[0] = 0.; pi[1] = 0.; pi[2] = 0.; pi[3] = 0.; pi[4] = 0.; pi[5] = 0.;
+ return;
+ }
+#endif
+ *rho = avg_rho
+ + lbfluid[0][0][index]
+ + lbfluid[0][1][index] + lbfluid[0][2][index]
+ + lbfluid[0][3][index] + lbfluid[0][4][index]
+ + lbfluid[0][5][index] + lbfluid[0][6][index]
+ + lbfluid[0][7][index] + lbfluid[0][8][index]
+ + lbfluid[0][9][index] + lbfluid[0][10][index]
+ + lbfluid[0][11][index] + lbfluid[0][12][index]
+ + lbfluid[0][13][index] + lbfluid[0][14][index]
+ + lbfluid[0][15][index] + lbfluid[0][16][index]
+ + lbfluid[0][17][index] + lbfluid[0][18][index];
+
+ j[0] = lbfluid[0][1][index] - lbfluid[0][2][index]
+ + lbfluid[0][7][index] - lbfluid[0][8][index]
+ + lbfluid[0][9][index] - lbfluid[0][10][index]
+ + lbfluid[0][11][index] - lbfluid[0][12][index]
+ + lbfluid[0][13][index] - lbfluid[0][14][index];
+ j[1] = lbfluid[0][3][index] - lbfluid[0][4][index]
+ + lbfluid[0][7][index] - lbfluid[0][8][index]
+ - lbfluid[0][9][index] + lbfluid[0][10][index]
+ + lbfluid[0][15][index] - lbfluid[0][16][index]
+ + lbfluid[0][17][index] - lbfluid[0][18][index];
+ j[2] = lbfluid[0][5][index] - lbfluid[0][6][index]
+ + lbfluid[0][11][index] - lbfluid[0][12][index]
+ - lbfluid[0][13][index] + lbfluid[0][14][index]
+ + lbfluid[0][15][index] - lbfluid[0][16][index]
+ - lbfluid[0][17][index] + lbfluid[0][18][index];
- if (calc_pi_flag) {
- local_pi[0] = avg_rho/3.0
- + local_n[1] + local_n[2]
- + local_n[7] + local_n[8] + local_n[9] + local_n[10]
- + local_n[11] + local_n[12] + local_n[13] + local_n[14];
- local_pi[2] = avg_rho/3.0
- + local_n[3] + local_n[4]
- + local_n[7] + local_n[8] + local_n[9] + local_n[10]
- + local_n[15] + local_n[16] + local_n[17] + local_n[18];
- local_pi[5] = avg_rho/3.0
- + local_n[5] + local_n[6]
- + local_n[11] + local_n[12] + local_n[13] + local_n[14]
- + local_n[15] + local_n[16] + local_n[17] + local_n[18];
- local_pi[1] = local_n[7] - local_n[9] + local_n[8] - local_n[10];
- local_pi[3] = local_n[11] + local_n[12] - local_n[13] - local_n[14];
- local_pi[4] = local_n[15] + local_n[16] - local_n[17] - local_n[18];
+ if (pi) {
+ pi[0] = avg_rho/3.0
+ + lbfluid[0][1][index] + lbfluid[0][2][index]
+ + lbfluid[0][7][index] + lbfluid[0][8][index]
+ + lbfluid[0][9][index] + lbfluid[0][10][index]
+ + lbfluid[0][11][index] + lbfluid[0][12][index]
+ + lbfluid[0][13][index] + lbfluid[0][14][index];
+ pi[2] = avg_rho/3.0
+ + lbfluid[0][3][index] + lbfluid[0][4][index]
+ + lbfluid[0][7][index] + lbfluid[0][8][index]
+ + lbfluid[0][9][index] + lbfluid[0][10][index]
+ + lbfluid[0][15][index] + lbfluid[0][16][index]
+ + lbfluid[0][17][index] + lbfluid[0][18][index];
+ pi[5] = avg_rho/3.0
+ + lbfluid[0][5][index] + lbfluid[0][6][index]
+ + lbfluid[0][11][index] + lbfluid[0][12][index]
+ + lbfluid[0][13][index] + lbfluid[0][14][index]
+ + lbfluid[0][15][index] + lbfluid[0][16][index]
+ + lbfluid[0][17][index] + lbfluid[0][18][index];
+ pi[1] = lbfluid[0][7][index] - lbfluid[0][9][index]
+ + lbfluid[0][8][index] - lbfluid[0][10][index];
+ pi[3] = lbfluid[0][11][index] + lbfluid[0][12][index]
+ - lbfluid[0][13][index] - lbfluid[0][14][index];
+ pi[4] = lbfluid[0][15][index] + lbfluid[0][16][index]
+ - lbfluid[0][17][index] - lbfluid[0][18][index];
}
-#else
+#else /* if not D3Q19 */
int i;
double tmp;
double (*c)[3] = lbmodel.c;
- *local_rho = 0.0;
+ *rho = 0.0;
- local_j[0] = 0.0;
- local_j[1] = 0.0;
- local_j[2] = 0.0;
+ j[0] = 0.0;
+ j[1] = 0.0;
+ j[2] = 0.0;
if (calc_pi_flag) {
- local_pi[0] = 0.0;
- local_pi[1] = 0.0;