Skip to content

Merge JS's changes in #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rupertnash wants to merge 2 commits into
base: bgk
Choose a base branch
from
Open

Merge JS's changes in #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4c642a0
Fixed some bugs:
joakimstenhammar Nov 14, 2016
357d1b6
Reverted to force densities instead of forces
joakimstenhammar Nov 24, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 0 additions & 1 deletion libd3q15/Lattice.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,5 +80,4 @@ typedef struct Site {
s.rho = L->rho_ptr + ijk; \
s.u = L->u_ptr + ijk*DQ_d; \
s.force = L->force_ptr + ijk*DQ_d; }

#endif
4 changes: 2 additions & 2 deletions libd3q15/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
LIBDIR = $(HOME)/lib
INCDIR = $(HOME)/include

CC = gcc
CFLAGS = -Wall -std=c99 -g -I$(HOME)/include
CC = icc
CFLAGS = -Wall -O3 -std=c99 -g -I$(HOME)/include

AR = ar rc

Expand Down
4 changes: 2 additions & 2 deletions libd3q15/bc_pbc.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ void bc_pbc_update(Lattice *lat) {

/* XY */
for (i[0] = 0; i[0] <= n[0]+1; i[0]++) {
for (i[1] = 0; i[1] <= n[0]+1; i[1]++) {
for (i[1] = 0; i[1] <= n[1]+1; i[1]++) {
i[2] = 0;
bc_pbc_do_site(lat, i, n);
i[2] = n[2] + 1;
Expand Down Expand Up @@ -72,7 +72,7 @@ void bc_pbc_do_site(Lattice *lat, int i[DQ_d], int n[DQ_d]) {

if (i[d]==0) {
src[d] = n[d];
} else if (i[d]==n[d]+1) {
} else if (i[d]==(n[d]+1)) {
src[d] = 1;
} else {
src[d] = i[d];
Expand Down
75 changes: 46 additions & 29 deletions libd3q15/d3q15.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,9 @@ Lattice *d3q15_init(int nx, int ny, int nz, double tau_s, double tau_b) {
lat->tau_s = tau_s;

lat->time_step = 0;


//printf("%5i %5i %5i %5i \n", lat->f_ptr, lat->rho_ptr, lat->u_ptr, lat->force_ptr);

/* Set our function pointers to NULL so we can know if they
* are initialised or not. */
lat->force_func = NULL;
Expand Down Expand Up @@ -78,19 +80,23 @@ void d3q15_destroy(Lattice *lat) {
void d3q15_iterate(Lattice *lat, int n_steps) {


double *totmass;
double *totmom[DQ_d];

for (int t=0; t<n_steps; t++) {

/* We have to work out the force first */
(*lat->force_func)(lat);

/* Collide */
collide(lat);

/* Update boundary */
(*lat->bc_func)(lat);

/* Propagate */
propagate(lat);

lat->time_step++;
}
}
Expand Down Expand Up @@ -133,9 +139,9 @@ void propagate (Lattice *lat) {
for (j=0; j<=ny+1; j++) {
for (k=0; k<=nz+1; k++) {
/* propagate */

/* [1,0,0] */
if (i<=nx)
if (i<=nx)
swap(DQ_f_get(lat, i,j,k, 1), DQ_f_get(lat, i+1,j,k, 2), tmp);
/* [0,1,0] */
if (j<=ny)
Expand All @@ -158,36 +164,38 @@ void propagate (Lattice *lat) {
swap(DQ_f_get(lat, i,j,k, 10),DQ_f_get(lat, i+1, j-1, k-1, 11), tmp);

/* reorder */
//if (i<=nx)
swap(DQ_f_get(lat, i,j,k, 1), DQ_f_get(lat, i,j,k, 2), tmp);
//if (i<=nx)
swap(DQ_f_get(lat, i,j,k, 1), DQ_f_get(lat, i,j,k, 2), tmp);
//if (j<=ny)
swap(DQ_f_get(lat, i,j,k, 3), DQ_f_get(lat, i,j,k, 4), tmp);
swap(DQ_f_get(lat, i,j,k, 3), DQ_f_get(lat, i,j,k, 4), tmp);
//if (k<=nz)
swap(DQ_f_get(lat, i,j,k, 5), DQ_f_get(lat, i,j,k, 6), tmp);
swap(DQ_f_get(lat, i,j,k, 5), DQ_f_get(lat, i,j,k, 6), tmp);

//if (i<=nz && j<=ny && k<=nz)
swap(DQ_f_get(lat, i,j,k, 7), DQ_f_get(lat, i,j,k, 14), tmp);
swap(DQ_f_get(lat, i,j,k, 7), DQ_f_get(lat, i,j,k, 14), tmp);
//if (i<=nz && j<=ny && k>0)
swap(DQ_f_get(lat, i,j,k, 8), DQ_f_get(lat, i,j,k, 13), tmp);
swap(DQ_f_get(lat, i,j,k, 8), DQ_f_get(lat, i,j,k, 13), tmp);
//if (i<=nz && j>0 && k<=nz)
swap(DQ_f_get(lat, i,j,k, 9), DQ_f_get(lat, i,j,k, 12), tmp);
swap(DQ_f_get(lat, i,j,k, 9), DQ_f_get(lat, i,j,k, 12), tmp);
//if (i<=nz && j>0 && k>0)
swap(DQ_f_get(lat, i,j,k, 10),DQ_f_get(lat, i,j,k, 11), tmp);
swap(DQ_f_get(lat, i,j,k, 10),DQ_f_get(lat, i,j,k, 11), tmp);

}
}
}
}

/************************************************************/

void calc_hydro_site(Site *site, Lattice* lat) {
void calc_hydro_site(Site *site, Lattice *lat) {
int i, a;
double rho = 0;
double rho = 0.0;
double mom[DQ_d];

for (a=0; a<DQ_d; a++) {
mom[a] = site->force[a] / 2.0;
mom[a] = site->force[a]*0.5;
}

for (i=0; i<DQ_q; i++) {
rho += site->f[i];
for (a=0; a<DQ_d; a++) {
Expand All @@ -204,7 +212,7 @@ void collide (Lattice *lat) {
/* loop over cells */
int i,j,k;
/* loop indices for dimension */
int a;
int a,b;
/* loop indices for velocities & modes */
int p;

Expand All @@ -214,39 +222,48 @@ void collide (Lattice *lat) {

const double tau_s = lat->tau_s;
const double omega_s = 1.0 / (tau_s + 0.5);

for (i=1; i<=lat->nx; i++) {
for (j=1; j<=lat->ny; j++) {
for (k=1; k<=lat->nz; k++) {
set_site(lat, site, i,j,k);
calc_hydro_site(&site, lat);

calc_hydro_site(&site, lat);
// rho & u evaluated at t
calc_equil(lat, site.rho[0], site.u, fEq);

calc_equil(lat, site.rho[0], site.u, fEq);

double u_F = 0.0;
for (a=0; a<DQ_d; a++)
for (a=0; a<DQ_d; a++) {
u_F += site.u[a]*site.force[a];
}

for (p=0; p<DQ_q; p++) {
double u_ep = 0.0;
double F_ep = 0.0;

for (a=0; a<DQ_d; a++) {
u_ep += site.u[a] * lat->xi[p][a];
F_ep += site.force[a] * lat->xi[p][a];
}

phi[p] = lat->w[p] * ((F_ep - u_ep) / lat->cs2 +
(u_ep * u_F) / (lat->cs2 * lat->cs2));


/* phi[p] = lat->w[p] * ((F_ep - u_ep) / lat->cs2 +
(u_ep * u_F) / (lat->cs2 * lat->cs2)); */ // Old version

phi[p] = lat->w[p] * ((F_ep - u_F) / lat->cs2 +
(u_ep * F_ep) / (lat->cs2 * lat->cs2)); // New version

// phi[p] = lat->w[p] * ((F_ep - u_F) / lat->cs2);

// phi[p] = lat->w[p] * (F_ep/lat->cs2); // Simplified version

/* Collide - Eq. (17) */
site.f[p] += omega_s * (tau_s * phi[p] - (site.f[p] - fEq[p]));
site.f[p] += omega_s * (tau_s * phi[p] - (site.f[p] - fEq[p]));

}
} /* k */
} /* j */
} /* i */

}

/************************************************************/
Expand All @@ -271,7 +288,7 @@ void calc_hydro(Lattice *lat) {
void calc_equil(Lattice *lat, double rho, double u[], double f_eq[]) {
int i, a, b;
const double cs2 = lat->cs2;

for (i = 0; i< DQ_q; i++) {
double feqi = 1.0;
for (a = 0; a < DQ_d; a++) {
Expand Down
5 changes: 3 additions & 2 deletions libd3q15/d3q15.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@
(k) * L->strides[DQ_Z]) * DQ_d + (m))
#define DQ_force_get(L, i,j,k, m) *(L->force_ptr + ((i) * L->strides[DQ_X] + \
(j) * L->strides[DQ_Y] + \
(k) * L->strides[DQ_Z])*DQ_d+(m))
(k) * L->strides[DQ_Z]) * DQ_d + (m))
#define DQ_access_macros
#include "eigenvectors.h"
#undef DQ_access_macros
Expand All @@ -42,7 +42,8 @@ void d3q15_destroy(Lattice *lat);
void propagate(Lattice *lat);
void collide(Lattice *lat);

void calc_equil(Lattice *lat, double rho, double *u, double *f_eq);
//void calc_equil(Lattice *lat, double rho, double *u, double *f_eq);
void calc_equil(Lattice *lat, double rho, double u[], double f_eq[]);
void calc_hydro(Lattice *lat);

void calc_phi(double force[], double rho, double u[], double ans[]);
Expand Down
8 changes: 5 additions & 3 deletions python/d3q15/Lattice.i
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -322,9 +322,11 @@ EXC_CHECK(force_set)
for (j=1; j<=$self->ny; j++)
for (k=1; k<=$self->nz; k++) {
double rho = DQ_rho_get($self,i,j,k);
for (a=0; a<DQ_d; a++)
u[a] = DQ_u_get($self, i,j,k, a) - DQ_force_get($self, i,j,k, a) / (2.0 * rho);

for (a=0; a<DQ_d; a++) {
/* u[a] = DQ_u_get($self, i,j,k, a) - DQ_force_get($self, i,j,k, a) / (2.0 * rho); */
u[a] = DQ_u_get($self, i,j,k, a);
/* if(i==2 && j==3 && a == 2) printf("Init %5i %5i %5i %15.10e %15.10e \n", i,j,k, u[a], DQ_force_get($self, i,j,k,a)); */
}
calc_equil($self, rho, u, &DQ_f_get($self, i,j,k,0));
}
}
Expand Down
1 change: 0 additions & 1 deletion python/dqTools/Controller.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,6 @@ def step(self):

"""
self.lat.step(1)

if self.isRecordStep():
self.log('Time step %d of %d' %
(self.lat.time_step, self.totalSteps))
Expand Down