Fix a bug in convolution. Also add a resample function to resample the grind...

Fix a bug in convolution. Also add a resample function to resample the grind data on different resolution lattice.

GridDxFormat/grid_reader.c

@@ -9,6 +9,7 @@
 #include <assert.h>
 #include <unistd.h>
 #include <ctype.h>
+#include <math.h>
 #include "grid_reader.h"
 typedef struct
 {
@@ -386,6 +387,130 @@ void set_grid_data(void* src, compute_type compute, void* kernel)
     } 
 }
 
+static void* create_new_grid(size_t nx, size_t ny, size_t nz, double* origin, double** basis)
+{
+    void* g = malloc(sizeof(grid));
+    ((grid*)g)->nx = nx;
+    ((grid*)g)->ny = ny;
+    ((grid*)g)->nz = nz;
+    ((grid*)g)->size = nx*ny*nz;
+    ((grid*)g)->data = (double*)malloc(sizeof(double)*(((grid*)g)->size)); 
+
+    ((grid*)g)->origin = (double*)malloc(sizeof(double)*3);
+    for(size_t i = 0; i < 3; ++i)
+        ((grid*)g)->origin[i] = origin[i];
+
+    ((grid*)g)->basis  = (double**)malloc(sizeof(double*)*3);
+    for(size_t i = 0; i < 3; ++i)
+        ((grid*)g)->basis[i] = (double*)malloc(sizeof(double)*3);
+
+    for(size_t i = 0; i < 3; ++i)
+        for(size_t j = 0; j < 3; ++j)
+            ((grid*)g)->basis[i][j] = basis[i][j];
+    return g;
+
+}
+
+static size_t Index(size_t nx, size_t ny, size_t nz, int x, int y, int z)
+{
+    x = (x < 0) ? 0 : x;
+    y = (y < 0) ? 0 : y;
+    z = (z < 0) ? 0 : z;
+    x = ((size_t)x >= nx) ? (int)nx-1 : x;
+    y = ((size_t)y >= ny) ? (int)ny-1 : y;
+    z = ((size_t)z >= nz) ? (int)nz-1 : z;
+    return z + nz * (y + ny * x);
+}
+/*
+static void inverse_3by3(double** m, double inv[3][3])
+{
+    double det = m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1])-
+                 m[0][1]*(m[1][0]*m[2][2]-m[1][2]*m[2][0])+ 
+                 m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0]);
+
+    inv[0][0] = (m[1][1]*m[2][2]-m[1][2]*m[2][1])/det;
+    inv[0][1] =-(m[0][1]*m[2][2]-m[0][2]*m[2][1])/det;
+    inv[0][2] = (m[0][1]*m[1][2]-m[0][2]*m[1][1])/det;
+
+    inv[1][0] =-(m[1][0]*m[2][2]-m[1][2]*m[2][0])/det;
+    inv[1][1] = (m[0][0]*m[2][2]-m[0][2]*m[2][0])/det;
+    inv[1][2] =-(m[0][0]*m[1][2]-m[0][2]*m[1][0])/det;
+
+    inv[2][0] = (m[1][0]*m[2][1]-m[1][1]*m[2][0])/det;
+    inv[2][1] =-(m[0][0]*m[2][1]-m[0][1]*m[2][0])/det;
+    inv[2][2] = (m[0][0]*m[1][1]-m[0][1]*m[1][0])/det;
+
+}*/
+//resample the grid data such that each dimensition is even number by using linear interpolation
+//I currently implement for diagonal basis set
+void* resample_grid_data(void* src)
+{
+    grid* old_g = (grid*)src;
+    grid* new_g;
+    size_t nx, ny, nz;
+    double** b = (double**)malloc(sizeof(double*)*3);
+    for(size_t i = 0; i < 3; ++i)
+    {
+        b[i] = (double*)malloc(sizeof(double)*3);
+        memset(b[i],0,sizeof(double)*3);
+    }
+
+    nx = (old_g->nx/2)*2;
+    ny = (old_g->ny/2)*2;
+    nz = (old_g->nz/2)*2;
+    
+    b[0][0] = old_g->basis[0][0]*(old_g->nx/(double)nx);
+    b[1][1] = old_g->basis[1][1]*(old_g->ny/(double)ny);
+    b[2][2] = old_g->basis[2][2]*(old_g->nz/(double)nz);
+
+    new_g = create_new_grid(nx,ny,nz,old_g->origin, b);
+
+    for(size_t i = 0; i < 3; ++i)
+        free(b[i]);
+    free(b);
+
+    for(size_t i = 0; i < new_g->size; ++i)
+    {
+        size_t ix = i / new_g->nz / new_g->ny;
+        size_t iy = i / new_g->nz % new_g->ny;
+        size_t iz = i % new_g->nz;
+
+        double posx = new_g->basis[0][0]*ix;
+        double posy = new_g->basis[1][1]*iy;
+        double posz = new_g->basis[2][2]*iz;
+
+        int x0 = floor(posx/old_g->basis[0][0]);
+        int y0 = floor(posy/old_g->basis[1][1]);
+        int z0 = floor(posz/old_g->basis[2][2]);
+ 
+        double posx0 = old_g->basis[0][0] * x0;
+        double posy0 = old_g->basis[1][1] * y0;
+        double posz0 = old_g->basis[2][2] * z0;
+
+        int x1 = x0 + 1;
+        int y1 = y0 + 1;
+        int z1 = z0 + 1;
+
+        double posx1 = old_g->basis[0][0] * x1;
+        double posy1 = old_g->basis[1][1] * y1;
+        double posz1 = old_g->basis[2][2] * z1;
+
+        double xd = (posx-posx0)/(posx1-posx0);
+        double yd = (posy-posy0)/(posy1-posy0);
+        double zd = (posz-posz0)/(posz1-posz0);
+        double c0 = ((1.-xd)* old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x0, y0, z0)] + 
+                         xd * old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x1, y0, z0)])*(1.-yd) +
+                    ((1.-xd)* old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x0, y1, z0)] + 
+                         xd * old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x1, y1, z0)])*yd;
+ 
+       double c1 = ((1.-xd)* old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x0, y0, z1)] +
+                        xd * old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x1, y0, z1)])*(1.-yd) +
+                   ((1.-xd)* old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x0, y1, z1)] + 
+                        xd * old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x1, y1, z1)])*yd;
+       new_g->data[i] = c0*(1.-zd) + c1*zd;
+    }
+    return new_g; 
+}
 double* get_grid_origin(void* src)
 {
     return ((grid*)src)->origin;

Kernel/lenard_jones_kernel.c

@@ -63,16 +63,16 @@ double lenard_jones_compute(size_t i, size_t j, size_t k, double* origin, double
     double pos[3];
     for(size_t dim = 0; dim < 3; ++dim)
         pos[dim] = basis[dim][0]*i + basis[dim][1]*j + basis[dim][2]*k + origin[dim];
-    double d = sqrt(pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]);
-    d = d*d*d*d*d*d;
+    double d = pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2];
+    d = d*d*d;
     if(d > 0.)
     {
         double e = -(kernel->c6)/d+(kernel->c12)/(d*d);
-        return (e>1e4) ? 1e4 : e;
+        return (e>20) ? 20 : e;
     }
     else
     {
-        return 1e4;
+        return 20;
     }
 }
 

inc/grid_reader.h

@@ -16,4 +16,5 @@ size_t get_size(void* src);
 void* copy_grid_format(void*);
 typedef double (*compute_type)(size_t,size_t,size_t,double*,double**, void*);
 void set_grid_data(void* src, compute_type compute, void* kernel);
+void* resample_grid_data(void* src);
 #endif

main.c

@@ -12,8 +12,11 @@ static void inner_product_grid(fftw_complex* dest, fftw_complex* src, size_t nx,
     {
         re = dest[i][0];
         im = dest[i][1];
-        dest[i][0] = (re*src[i][0]-im*src[i][1])/V;
-        dest[i][1] = (re*src[i][1]+im*src[i][0])/V;
+        size_t x = i / (nz/2+1) / ny; 
+        size_t y = i / (nz/2+1) % ny;
+        size_t z = i % (nz/2+1);
+        dest[i][0] = (re*src[i][0]-im*src[i][1])/V*((0.5-(x+y+z)%2)*2);
+        dest[i][1] = (re*src[i][1]+im*src[i][0])/V*((0.5-(x+y+z)%2)*2);
     }
 }
 static void shift_zero_frequency(size_t nx, size_t ny, size_t nz, double* data)
@@ -75,7 +78,11 @@ int main(int argc, char** argv)
     char* kernel_name   = argv[2];
     char* kernel_infile = argv[3];
 
-    void* density   = gridformat_reader(density_file);
+    void* density          = gridformat_reader(density_file);
+    void* resample_density = resample_grid_data(density);
+    grid_destroy(density);
+    density = resample_density;
+
     void* potential = copy_grid_format(density);
    
     create_type creator = find_list(create_map, kernel_name);