From 4c642a0abfe9ff525f153b1e0c75ad3f11320517 Mon Sep 17 00:00:00 2001 From: Joakim Stenhammar Date: Mon, 14 Nov 2016 15:59:10 +0100 Subject: [PATCH 1/2] Fixed some bugs: 1. Corrected forcing term (d3q15.c) 2. Fixed inconsistecy between forces and force densities (d3q15.c) 3. Included forcing in initialization (Lattice.i) 4. Fixed indexing error in PBCs (bc_pbc.c) --- libd3q15/Makefile | 4 +-- libd3q15/bc_pbc.c | 2 +- libd3q15/d3q15.c | 55 +++++++++++++++++++++--------------- python/d3q15/Lattice.i | 6 ++-- python/dqTools/Controller.py | 1 - 5 files changed, 38 insertions(+), 30 deletions(-) diff --git a/libd3q15/Makefile b/libd3q15/Makefile index c500b08..e317809 100644 --- a/libd3q15/Makefile +++ b/libd3q15/Makefile @@ -1,8 +1,8 @@ LIBDIR = $(HOME)/lib INCDIR = $(HOME)/include -CC = gcc -CFLAGS = -Wall -std=c99 -g -I$(HOME)/include +CC = icc +CFLAGS = -Wall -O2 -std=c99 -g -I$(HOME)/include AR = ar rc diff --git a/libd3q15/bc_pbc.c b/libd3q15/bc_pbc.c index 3468250..762efad 100644 --- a/libd3q15/bc_pbc.c +++ b/libd3q15/bc_pbc.c @@ -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; diff --git a/libd3q15/d3q15.c b/libd3q15/d3q15.c index 1a6b30f..3537046 100644 --- a/libd3q15/d3q15.c +++ b/libd3q15/d3q15.c @@ -133,9 +133,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) @@ -158,21 +158,21 @@ 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); } } } @@ -184,8 +184,9 @@ void calc_hydro_site(Site *site, Lattice* lat) { int i, a; double rho = 0; double mom[DQ_d]; + for (a=0; aforce[a] / 2.0; + mom[a] = 0.0; } for (i=0; irho[0] = rho; + for (a=0; au[a] = mom[a] / rho; + site->u[a] = mom[a] / rho + site->force[a]/2.0; } } + 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; @@ -214,19 +217,20 @@ 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); // rho & u evaluated at t + calc_hydro_site(&site, lat); calc_equil(lat, site.rho[0], site.u, fEq); double u_F = 0.0; - for (a=0; axi[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] = site.rho[0]*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/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 */ - } /************************************************************/ @@ -271,7 +280,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++) { diff --git a/python/d3q15/Lattice.i b/python/d3q15/Lattice.i index 77b4085..3121b8a 100644 --- a/python/d3q15/Lattice.i +++ b/python/d3q15/Lattice.i @@ -322,9 +322,9 @@ 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 Date: Thu, 24 Nov 2016 11:04:00 +0100 Subject: [PATCH 2/2] Reverted to force densities instead of forces --- libd3q15/Lattice.h | 1 - libd3q15/Makefile | 2 +- libd3q15/bc_pbc.c | 2 +- libd3q15/d3q15.c | 38 +++++++++++++++++++++++--------------- libd3q15/d3q15.h | 5 +++-- python/d3q15/Lattice.i | 4 +++- 6 files changed, 31 insertions(+), 21 deletions(-) diff --git a/libd3q15/Lattice.h b/libd3q15/Lattice.h index bd53ff1..253c56f 100644 --- a/libd3q15/Lattice.h +++ b/libd3q15/Lattice.h @@ -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 diff --git a/libd3q15/Makefile b/libd3q15/Makefile index e317809..5c24566 100644 --- a/libd3q15/Makefile +++ b/libd3q15/Makefile @@ -2,7 +2,7 @@ LIBDIR = $(HOME)/lib INCDIR = $(HOME)/include CC = icc -CFLAGS = -Wall -O2 -std=c99 -g -I$(HOME)/include +CFLAGS = -Wall -O3 -std=c99 -g -I$(HOME)/include AR = ar rc diff --git a/libd3q15/bc_pbc.c b/libd3q15/bc_pbc.c index 762efad..e17f82a 100644 --- a/libd3q15/bc_pbc.c +++ b/libd3q15/bc_pbc.c @@ -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]; diff --git a/libd3q15/d3q15.c b/libd3q15/d3q15.c index 3537046..88f20e2 100644 --- a/libd3q15/d3q15.c +++ b/libd3q15/d3q15.c @@ -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; @@ -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; tforce_func)(lat); /* Collide */ collide(lat); - + /* Update boundary */ (*lat->bc_func)(lat); /* Propagate */ propagate(lat); - + lat->time_step++; } } @@ -172,7 +178,8 @@ void propagate (Lattice *lat) { //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); //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); + } } } @@ -180,15 +187,15 @@ void propagate (Lattice *lat) { /************************************************************/ -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; aforce[a]*0.5; } - + for (i=0; if[i]; for (a=0; arho[0] = rho; - for (a=0; au[a] = mom[a] / rho + site->force[a]/2.0; + site->u[a] = mom[a] / rho; } } - void collide (Lattice *lat) { /* loop over cells */ int i,j,k; @@ -222,11 +227,12 @@ void collide (Lattice *lat) { 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); // rho & u evaluated at t - calc_hydro_site(&site, lat); - 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; axi[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)); */ // Old version - phi[p] = site.rho[0]*lat->w[p] * ((F_ep - u_F) / lat->cs2 + + 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) */ diff --git a/libd3q15/d3q15.h b/libd3q15/d3q15.h index 6c40ca5..6be19ce 100644 --- a/libd3q15/d3q15.h +++ b/libd3q15/d3q15.h @@ -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 @@ -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[]); diff --git a/python/d3q15/Lattice.i b/python/d3q15/Lattice.i index 3121b8a..2fc9c4e 100644 --- a/python/d3q15/Lattice.i +++ b/python/d3q15/Lattice.i @@ -323,7 +323,9 @@ EXC_CHECK(force_set) for (k=1; k<=$self->nz; k++) { double rho = DQ_rho_get($self,i,j,k); for (a=0; a