Skip to content

Commit 60abc4b

Browse files
author
Davoud Saffar
committed
Merge branch 'master' into kaiser-all
2 parents 7a4d340 + 7de2be7 commit 60abc4b

File tree

95 files changed

+3553
-7
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

95 files changed

+3553
-7
lines changed

‎CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,7 @@ add_subdirectory("${PROJECT_SOURCE_DIR}/SE_Stresslet/mex")
7171
add_subdirectory("${PROJECT_SOURCE_DIR}/SE_Rotlet/mex")
7272
add_subdirectory("${PROJECT_SOURCE_DIR}/SE2P_Stokes/mex")
7373
add_subdirectory("${PROJECT_SOURCE_DIR}/SE1P/mex")
74+
add_subdirectory("${PROJECT_SOURCE_DIR}/SE0P/mex")
7475

7576
# Enable running Matlab test suite through ctest
7677
enable_testing()

‎SE0P/find_params.m

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
function [M, P, rc, sl, nl, R] = find_params(L, Q, N, xi, tol)
2+
3+
[M, P, rc, ~] = ewald_k_inf2(L, Q, N, xi, tol*10);
4+
5+
for nl = 1:max(M/10,5);
6+
v = erfc(nl/2/xi);
7+
if v<tol;
8+
break;
9+
end
10+
end
11+
12+
for sl = 2:3
13+
if exp(-2*pi*nl*sl)<tol
14+
break
15+
end
16+
end
17+
18+
R = L*sqrt(2)*(1+P/M/2);
19+
20+
for sl = 1:4
21+
v = 100/sqrt(M)*exp(-50/sqrt(M)*sl);
22+
if v<tol
23+
break
24+
end
25+
end
26+
27+
for nl = 1:max(M/10,5);
28+
v = 2/sqrt(M)*exp(-160/M*nl);
29+
if v<tol;
30+
break;
31+
end
32+
end
33+
sl = sl -1 ;
34+
35+
36+
nl = -log10(tol);
37+
if(tol>=1e-3)
38+
sl = 1;
39+
elseif tol>=1e-6
40+
sl = 2;
41+
elseif tol>=1e-12
42+
sl = 3;
43+
else
44+
sl = 4;
45+
end
46+
47+
nl = 1:20;
48+
f = exp(-nl*M/(xi*L));
49+
idx = find(f<tol);
50+
nl = min(idx);

‎SE0P/init.m

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
base = '';
2+
addpath(base)
3+
addpath([base 'tests']);
4+
addpath([base 'src']);
5+
addpath([base 'mex']);
6+
addpath([base '../SE_fast_gridding']);
7+
addpath([base '../bin'])
8+
addpath([base '../util'])

‎SE0P/laplace_demo.m

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
clear
2+
rng(1);
3+
4+
N = 10000;
5+
L = 3;
6+
box = [L L L];
7+
[x, f] = laplace_system(N, box);
8+
9+
M0 = 28; % Set M0, the rest is auto
10+
11+
opt.M = M0*box;
12+
opt.xi = pi*M0 / 12;
13+
opt.P = 24;
14+
opt.oversampling = 3.2;
15+
opt.rc = 6 / opt.xi;
16+
opt.box = box
17+
18+
% Precomp
19+
disp('Precomputing...')
20+
tic
21+
pre = laplace_precomp(opt);
22+
toc
23+
24+
25+
% Ewald
26+
disp('Free-space Ewald...')
27+
tic
28+
uf = laplace_fourier_space(x, f, opt, pre);
29+
ur = laplace_real_space(x, f, opt);
30+
us = -f*opt.xi *2/sqrt(pi);
31+
ue = uf+ur+us;
32+
toc
33+
% Direct
34+
disp('Direct sum...')
35+
tic
36+
u = laplace_direct(x, f, box);
37+
toc
38+
39+
error = norm(u-ue, inf) / norm(u,inf)

‎SE0P/laplace_demo_force.m

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
clear
2+
rng(1);
3+
4+
N = 10000;
5+
L = 3;
6+
box = [L L L];
7+
[x, f] = laplace_system(N, box);
8+
9+
M0 = 28; % Set M0, the rest is auto
10+
11+
opt.M = M0*box;
12+
opt.xi = pi*M0 / 12;
13+
opt.P = 24;
14+
opt.oversampling = 3;
15+
opt.rc = 6 / opt.xi;
16+
opt.box = box
17+
18+
% Precomp
19+
disp('Precomputing...')
20+
tic
21+
pre = laplace_precomp(opt);
22+
toc
23+
24+
25+
% Ewald
26+
disp('Free-space Ewald...')
27+
tic
28+
uf = laplace_fourier_space(x, f, opt, pre,'force');
29+
ur = laplace_real_space_force(x, f, opt);
30+
us = -f*opt.xi *2/sqrt(pi)*0;
31+
ue = uf+ur+us;
32+
toc
33+
% Direct
34+
disp('Direct sum...')
35+
tic
36+
u = laplace_direct_force(x, f, box);
37+
toc
38+
39+
error = norm(u-ue, inf) / norm(u,inf)

‎SE0P/mex/CMakeLists.txt

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
## MEX functions
2+
3+
matlab_add_mex(
4+
NAME SE0P_stokeslet_fast_fs_k_scaling
5+
SRC stokeslet_fast_fs_k_scaling.c
6+
)
7+
8+
matlab_add_mex(
9+
NAME SE0P_stokeslet_direct_mex
10+
SRC stokeslet_direct.c
11+
)
12+
13+
matlab_add_mex(
14+
NAME SE0P_stokeslet_rsrc_mex
15+
SRC stokeslet_rsrc.c
16+
)
17+
18+
matlab_add_mex(
19+
NAME SE0P_stresslet_fast_fs_k_scaling
20+
SRC stresslet_fast_fs_k_scaling.c
21+
)
22+
23+
matlab_add_mex(
24+
NAME SE0P_stresslet_direct_mex
25+
SRC stresslet_direct.c
26+
)
27+
28+
matlab_add_mex(
29+
NAME SE0P_stresslet_rsrc_mex
30+
SRC stresslet_rsrc.c
31+
)
32+
33+
matlab_add_mex(
34+
NAME SE0P_rotlet_direct_mex
35+
SRC rotlet_direct.c
36+
)
37+
38+
matlab_add_mex(
39+
NAME SE0P_rotlet_rsrc_mex
40+
SRC rotlet_rsrc.c
41+
)
42+
43+
if(MKL_FOUND)
44+
45+
# Build stokeslet RS using cell list and MKL
46+
matlab_add_mex(
47+
NAME SE0P_stokeslet_rsrc_cell_mkl_mex
48+
SRC stokeslet_rsrc_cell.c cell_list.c
49+
)
50+
51+
# Linking twice, this might have to do with circular dependencies
52+
target_link_libraries(SE0P_stokeslet_rsrc_cell_mkl_mex ${MKL_STATIC})
53+
target_link_libraries(SE0P_stokeslet_rsrc_cell_mkl_mex ${MKL_STATIC})
54+
55+
target_compile_definitions(SE0P_stokeslet_rsrc_cell_mkl_mex PUBLIC INTEL_MKL)
56+
57+
endif()
58+
59+
# no-MKL version
60+
matlab_add_mex(
61+
NAME SE0P_stokeslet_rsrc_cell_mex
62+
SRC stokeslet_rsrc_cell.c cell_list.c
63+
)

‎SE0P/mex/cell_list.c

Lines changed: 88 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,88 @@
1+
#include "cell_list.h"
2+
#include "math.h"
3+
#include "mex_compat.h"
4+
5+
// =============================================================================
6+
// ==== BUILD CELL LIST (NEW) ==================================================
7+
//
8+
// Alternative implementation that returns a straight list
9+
// Reads faster, but requires two sweeps for creation
10+
//
11+
// TODO: Add some assertions to make sure rc not too big,
12+
// and that box can be divided into square cells.
13+
void build_cell_list(
14+
// Input
15+
const double* restrict x,
16+
const int N,
17+
const double* restrict box,
18+
const double rc,
19+
// Output
20+
double* rn_p,
21+
int ncell[3],
22+
int* restrict *cell_list_p,
23+
int* restrict *cell_idx_p
24+
)
25+
{
26+
int i,j;
27+
int ncell_tot;
28+
int icell[3];
29+
double boxmin, rn;
30+
// Outputs
31+
int* restrict cell_list;
32+
int* restrict cell_idx;
33+
// Intermediates (could do this with fewer vars, but this is clear)
34+
int* restrict cell_count;
35+
int* restrict point_cell_map;
36+
37+
38+
// Setup cell partitioning
39+
boxmin = box[0];
40+
41+
if(box[1]<boxmin)
42+
boxmin = box[1];
43+
if (box[2]<boxmin)
44+
boxmin = box[2];
45+
rn = boxmin / floor(boxmin/rc);
46+
for(i=0;i<3;i++)
47+
ncell[i] = round( box[i]/rn ) + 1;
48+
ncell_tot = ncell[0]*ncell[1]*ncell[2];
49+
50+
// Prepare arrays
51+
cell_list = __MALLOC(N*sizeof(int));
52+
cell_idx = __CALLOC(ncell_tot+1, sizeof(int));
53+
point_cell_map = __MALLOC(N*sizeof(int));
54+
55+
// Build list in two sweeps,
56+
// First sweep, count number of points in each cell
57+
for(i=0; i<N; i++)
58+
{
59+
for(j=0; j<3; j++)
60+
icell[j] = x[i*3+j]/rn;
61+
int icell_idx =
62+
icell[0] +
63+
icell[1]*ncell[0] +
64+
icell[2]*ncell[1]*ncell[0];
65+
//ASSERT(icell_idx < ncell_tot, 'cell index out of bounds, point x=box?');
66+
cell_idx[icell_idx + 1]++;
67+
point_cell_map[i] = icell_idx;
68+
}
69+
// Generate adressing
70+
cell_idx[0]=0;
71+
for (int i=0; i<ncell_tot; i++)
72+
cell_idx[i+1] += cell_idx[i];
73+
// Setup new vector
74+
cell_count = __CALLOC(ncell_tot, sizeof(int));
75+
// Second sweep, build list
76+
for(i=0; i<N; i++)
77+
{
78+
int icell_idx = point_cell_map[i];
79+
int adr = cell_idx[icell_idx] + cell_count[icell_idx];
80+
cell_list[adr] = i;
81+
cell_count[icell_idx]++;
82+
}
83+
__FREE(cell_count);
84+
__FREE(point_cell_map);
85+
*rn_p = rn;
86+
*cell_list_p = cell_list;
87+
*cell_idx_p = cell_idx;
88+
}

‎SE0P/mex/cell_list.h

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
#ifndef __CELL_LIST_H__
2+
#define __CELL_LIST_H__
3+
4+
void build_cell_list(
5+
// Input
6+
const double* restrict x,
7+
const int N,
8+
const double* restrict box,
9+
const double rc,
10+
// Output
11+
double* rn_p,
12+
int ncell[3],
13+
int* restrict *cell_list_p,
14+
int* restrict *cell_idx_p
15+
);
16+
17+
#endif

‎SE0P/mex/mex_compat.h

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
#ifndef __MEX_COMPAT_H__
2+
#define __MEX_COMPAT_H__
3+
4+
#define GiB *((size_t) 1024*1024*1024)
5+
6+
/** Maximum total amount of memory to malloc in applications,
7+
not strictly followed, but tries to avoid some swap-deaths. */
8+
#define MALLOC_MAX 6 GiB
9+
10+
11+
#ifdef MATLAB_MEX_FILE
12+
#include "mex.h"
13+
#define __MALLOC mxMalloc
14+
#define __CALLOC mxCalloc
15+
#define __REALLOC mxRealloc
16+
#define __FREE mxFree
17+
#define __PRINTF mexPrintf
18+
#define __ERROR(msg) { \
19+
char ebuffer[1024]; \
20+
snprintf(ebuffer, 1024, "%s at %s, line %d.\n", msg, __FILE__, __LINE__); \
21+
mexErrMsgTxt(ebuffer); \
22+
}
23+
#define __WARNING(msg) { \
24+
char wbuffer[1024]; \
25+
snprintf(wbuffer, 1024, "%s at %s, line %d.\n", msg, __FILE__, __LINE__); \
26+
mexWarnMsgTxt(wbuffer); \
27+
}
28+
#define ASSERT(expr, err) \
29+
if(! (expr)) \
30+
__ERROR("Assertion failed: (" # expr ") " # err);
31+
#define __FLUSH() mexEvalString("drawnow;")
32+
#else
33+
#include <stdlib.h>
34+
#include <stdio.h>
35+
#include "assert.h"
36+
#define __MALLOC malloc
37+
#define __CALLOC calloc
38+
#define __REALLOC realloc
39+
#define __FREE free
40+
#define __PRINTF printf
41+
#define __ERROR(msg) assert(0 && msg)
42+
#define __WARNING(msg) printf("[WARNING] %s at %s, line %d.\n", msg, __FILE__, __LINE__);
43+
#define ASSERT(expr, err) \
44+
assert( (expr) && "Assert (" # expr ") failed: " # err "\n");
45+
#define __FLUSH()
46+
#endif
47+
48+
#endif

0 commit comments

Comments
 (0)