Wall tests

solvers/DPStokes/mobility.h

@@ -51,23 +51,52 @@ class DPStokes: public libmobility::Mobility{
     this->lanczosTolerance = ipar.tolerance;
     this->dppar.mode = this->wallmode;
     this->dppar.hydrodynamicRadius = ipar.hydrodynamicRadius[0];
-    this->dppar.w = 6;
-    this->dppar.beta = 1.714*this->dppar.w;
-    real h = this->dppar.hydrodynamicRadius/1.554;
+    this->dppar.w = 4;
+    this->dppar.beta = 1.785*this->dppar.w;
+    real h = this->dppar.hydrodynamicRadius/1.205;
     this->dppar.alpha = this->dppar.w/2.0;
     this->dppar.tolerance = 1e-6;
-    this->dppar.nx = int(this->dppar.Lx/h + 0.5);
-    this->dppar.ny = int(this->dppar.Ly/h + 0.5);
-    // Add a buffer of w*h/2 when there is an open boundary
+
+    // adjust box size to be a multiple of h
+    real N_in = this->dppar.Lx/h;
+    int N_up = ceil(N_in);
+    int N_down = floor(N_in);
+    int N;
+    // either N_up or N_down will be a multiple of 2. pick the even one for a more FFT friendly grid.
+    if(N_up % 2 == 0){
+      N = N_up;
+    }else{
+      N = N_down;
+    }
+
+    // note: only set up for square boxes
+    this->dppar.Lx = N*h;
+    this->dppar.Ly = N*h;
+    this->dppar.nx = N;
+    this->dppar.ny = N;
+
+    // Add a buffer of 1.5*w*h/2 when there is an open boundary
     if(this->wallmode == "nowall"){
-      this->dppar.zmax += this->dppar.w*h/2;
-      this->dppar.zmin -= this->dppar.w*h/2;
+      this->dppar.zmax += 1.5*this->dppar.w*h/2;
+      this->dppar.zmin -= 1.5*this->dppar.w*h/2;
     }
     if(this->wallmode == "bottom"){
-      this->dppar.zmin -= this->dppar.w*h/2;
+      this->dppar.zmax += 1.5*this->dppar.w*h/2;
     }
     real Lz = this->dppar.zmax - this->dppar.zmin;
-    this->dppar.nz = M_PI*Lz/(h);
+    real H = Lz/2;
+    // sets chebyshev node spacing at its coarsest (in the middle) to be h
+    real nz_actual = M_PI/(asin(h/H)) + 1;
+
+    // pick nearby N such that 2(Nz-1) is FFT friendly
+    N_up = ceil(nz_actual);
+    N_down = floor(nz_actual);
+    if(N_up % 2 == 1){
+      this->dppar.nz = N_up;
+    } else {
+      this->dppar.nz = N_down;
+    }
+
     dpstokes->initialize(dppar, this->numberParticles);
     Mobility::initialize(ipar);
   }

solvers/NBody/extra/hydrodynamicKernels.cuh

@@ -70,21 +70,21 @@ namespace nbody_rpy{
 	real invZi = real(1.0) / hj;
 	real invZi3 = invZi * invZi * invZi;
 	real invZi5 = invZi3 * invZi * invZi;
-	correction.x += -vj.x*(real(9.0)*invZi - real(2.0)*invZi3 + invZi5 ) / real(12.0);
-	correction.y += -vj.y*(real(9.0)*invZi - real(2.0)*invZi3 + invZi5 ) / real(12.0);
-	correction.z += -vj.z*(real(9.0)*invZi - real(4.0)*invZi3 + invZi5 ) / real(6.0);
+	correction.x += -vj.x*(real(9.0)*invZi - real(2.0)*invZi3 + invZi5 ) / real(16.0);
+	correction.y += -vj.y*(real(9.0)*invZi - real(2.0)*invZi3 + invZi5 ) / real(16.0);
+	correction.z += -vj.z*(real(9.0)*invZi - real(4.0)*invZi3 + invZi5 ) / real(8.0);
       }
       else{
 	real h_hat = hj / rij.z;
 	real invR = rsqrt(dot(rij, rij));
 	real3 e = rij*invR;
 	real invR3 = invR * invR * invR;
 	real invR5 = invR3 * invR * invR;
-	real fact1 = -(real(3.0)*(real(1.0)+real(2.0)*h_hat*(real(1.0)-h_hat)*e.z*e.z) * invR + real(2.0)*(real(1.0)-real(3.0)*e.z*e.z) * invR3 - real(2.0)*(real(1.0)-real(5.0)*e.z*e.z) * invR5)  / real(3.0);
-	real fact2 = -(real(3.0)*(real(1.0)-real(6.0)*h_hat*(real(1.0)-h_hat)*e.z*e.z) * invR - real(6.0)*(real(1.0)-real(5.0)*e.z*e.z) * invR3 + real(10.0)*(real(1.0)-real(7.0)*e.z*e.z) * invR5) / real(3.0);
-	real fact3 =  e.z * (real(3.0)*h_hat*(real(1.0)-real(6.0)*(real(1.0)-h_hat)*e.z*e.z) * invR - real(6.0)*(real(1.0)-real(5.0)*e.z*e.z) * invR3 + real(10.0)*(real(2.0)-real(7.0)*e.z*e.z) * invR5) * real(2.0) / real(3.0);
-	real fact4 =  e.z * (real(3.0)*h_hat*invR - real(10.0)*invR5) * real(2.0) / real(3.0);
-	real fact5 = -(real(3.0)*h_hat*h_hat*e.z*e.z*invR + real(3.0)*e.z*e.z*invR3 + (real(2.0)-real(15.0)*e.z*e.z)*invR5) * real(4.0) / real(3.0);
+	real fact1 = -(real(3.0)*(real(1.0)+real(2.0)*h_hat*(real(1.0)-h_hat)*e.z*e.z) * invR + real(2.0)*(real(1.0)-real(3.0)*e.z*e.z) * invR3 - real(2.0)*(real(1.0)-real(5.0)*e.z*e.z) * invR5)  / real(4.0);
+	real fact2 = -(real(3.0)*(real(1.0)-real(6.0)*h_hat*(real(1.0)-h_hat)*e.z*e.z) * invR - real(6.0)*(real(1.0)-real(5.0)*e.z*e.z) * invR3 + real(10.0)*(real(1.0)-real(7.0)*e.z*e.z) * invR5) / real(4.0);
+	real fact3 =  e.z * (real(3.0)*h_hat*(real(1.0)-real(6.0)*(real(1.0)-h_hat)*e.z*e.z) * invR - real(6.0)*(real(1.0)-real(5.0)*e.z*e.z) * invR3 + real(10.0)*(real(2.0)-real(7.0)*e.z*e.z) * invR5) / real(2.0);
+	real fact4 =  e.z * (real(3.0)*h_hat*invR - real(10.0)*invR5) / real(2.0);
+	real fact5 = -(real(3.0)*h_hat*h_hat*e.z*e.z*invR + real(3.0)*e.z*e.z*invR3 + (real(2.0)-real(15.0)*e.z*e.z)*invR5);
 	correction.x += (fact1 + fact2 * e.x*e.x)*vj.x;
 	correction.x += (fact2 * e.x*e.y)*vj.y;
 	correction.x += (fact2 * e.x*e.z + fact3 * e.x)*vj.z;

tests/ref/readme.txt

@@ -0,0 +1,17 @@
+Some of this data is generated from the original DoublyPeriodicStokes repository as a correctness check. Some of this data was generated from this repository as a consistency check.
+
+The reason we've regenerated data is that default parameter selection has been changed to fix the grid parameters to the values that best preserve the hydrodynamic radius. When data has been regenerated, consistency was checked by first matching the original DoublyPeriodicStokes repository by using their same parameters, then switching to the new parameter selection and regenerating.
+
+Abbreviations: 
+bw- bottom wall
+sc- slit channel
+
+-------------- Files --------------
+SELF:
+self_mobility_bw_ref: from the DPStokes repository. Periodized NBody using kernels with wall corrections in double precision.
+self_mobility_bw_ref_noimg.mat: from the DPStokes repository. NBody kernels with wall corrections in double precision, but without using periodized images.
+self_mobility_bw_w4, self_mobility_sc_w4: generated from this repository in double precision with default parameters. 
+
+PAIR:
+pair_mobility_bw_w4, pair_mobility_sc_w4: generated from this repository in double precision with default parameters. 
+pair_mobility_bw_ref_noimg: from the DPStokes repository. NBody kernels with wall corrections in double precision, but without using periodized images.

tests/test_fluctuation_dissipation.py

@@ -62,7 +62,7 @@ def fluctuation_dissipation_KS(M, fluctuation_method):
     ],
 )
 @pytest.mark.parametrize("hydrodynamicRadius", [0.95, 1.12])
-@pytest.mark.parametrize("numberParticles", [1, 2, 3, 10])
+@pytest.mark.parametrize("numberParticles", [10])
 def test_fluctuation_dissipation(
     Solver, periodicity, hydrodynamicRadius, numberParticles
 ):

tests/test_mobility_matrix.py

@@ -35,8 +35,8 @@ def test_mobility_matrix(Solver, periodicity, hydrodynamicRadius, numberParticle
     assert M.dtype == precision
     sym = M - M.T
     assert np.allclose(
-        sym, 0.0, rtol=0, atol=1e-7
-    ), f"Mobility matrix is not symmetric within 1e-6, max diff: {np.max(np.abs(sym))}"
+        sym, 0.0, rtol=0, atol=5e-5
+    ), f"Mobility matrix is not symmetric within 1e-5, max diff: {np.max(np.abs(sym))}"
 
 
 def test_self_mobility_selfmobility():
@@ -62,7 +62,7 @@ def test_self_mobility_selfmobility():
 
 
 def test_self_mobility_nbody():
-    # Mobility should be just 1/(6\pi\eta R)*(1 - 2.83729748/L) for a single particle.
+    # Mobility should be just 1/(6\pi\eta R)
     Solver = NBody
     precision = np.float32 if Solver.precision == "float" else np.float64
     solver = Solver("open", "open", "open")

tests/test_wall_mobility.py

@@ -0,0 +1,164 @@
+import pytest
+import numpy as np
+import scipy.interpolate
+import scipy.io
+
+from libMobility import DPStokes, NBody
+from utils import compute_M
+
+self_mobility_params = {
+    "DPStokes": {"dt": 1, "Lx": 76.8, "Ly": 76.8, "zmin": 0, "zmax": 19.2},
+    "NBody": {"algorithm": "advise"},
+}
+
+@pytest.mark.parametrize(
+    ("Solver", "periodicity", "tol", "start_height", "ref_file"),
+    [
+        (DPStokes, ("periodic", "periodic", "single_wall"), 5e-3, 4, "self_mobility_bw_ref.mat"), # correctness check
+        (DPStokes, ("periodic", "periodic", "single_wall"), 1e-6, 0, "self_mobility_bw_w4.mat"), # consistency check
+        (DPStokes, ("periodic", "periodic", "two_walls"), 1e-6, 0, "self_mobility_sc_w4.mat"),
+        (NBody, ("open", "open", "single_wall"), 1e-6, 1, "self_mobility_bw_ref_noimg.mat")
+    ],
+)
+def test_self_mobility(Solver, periodicity, tol, start_height, ref_file):
+    zmax = 19.2
+    xymax = 76.8
+    params = self_mobility_params[Solver.__name__]
+
+    ref_dir = "./ref/"
+    ref = scipy.io.loadmat(ref_dir + ref_file)
+    refM = ref['M']
+    refHeights = ref['heights'][0]
+
+    hydrodynamicRadius = 1.0
+    eta = 1/4/np.sqrt(np.pi)
+
+    precision = np.float32 if Solver.precision == "float" else np.float64
+
+    solver = Solver(*periodicity)
+    solver.setParameters(**params)
+    numberParticles = 1
+    solver.initialize(
+        temperature=0,
+        viscosity=eta,
+        hydrodynamicRadius=hydrodynamicRadius,
+        numberParticles=numberParticles,
+    )
+
+    start_ind = np.where(refHeights >= start_height)[0][0]
+    refHeights = refHeights[start_ind:]
+    nHeights = len(refHeights)
+
+    refM = refM[start_ind:]
+
+    normMat = np.ones((3*numberParticles, 3*numberParticles), dtype=precision)
+    diag_ind = np.diag_indices_from(normMat)
+    normMat[diag_ind] = 1/(6*np.pi*eta*hydrodynamicRadius) # only for diag elements as this is for self mobility
+
+    allM = np.zeros((nHeights, 3*numberParticles, 3*numberParticles), dtype=precision)
+    for i in range(0,nHeights):
+        positions = np.array([[xymax/2, xymax/2, refHeights[i]]], dtype=precision)
+        solver.setPositions(positions)
+
+        M = compute_M(solver, numberParticles)
+        M /= normMat
+        allM[i] = M
+
+    # uncomment to save datafile for test plots
+    # scipy.io.savemat('./temp/test_data/test_' + ref_file, {'M': allM, 'heights': refHeights})
+
+    diags = [np.diag(matrix) for matrix in allM]
+    ref_diags = [np.diag(matrix)[0:3] for matrix in refM] # only take diagonal elements from forces
+
+    diff = np.abs([(diag - ref_diag) for diag, ref_diag in zip(diags, ref_diags)])
+
+    avgErr = np.mean(diff)
+
+    assert avgErr < tol, "Self mobility does not match reference"
+
+@pytest.mark.parametrize(
+    ("Solver", "periodicity", "tol", "ref_file"),
+    [
+        (DPStokes, ("periodic", "periodic", "single_wall"), 1e-6, "pair_mobility_bw_w4.mat"),
+        (DPStokes, ("periodic", "periodic", "two_walls"), 1e-6, "pair_mobility_sc_w4.mat"),
+        (NBody, ("open", "open", "single_wall"), 1e-4, "pair_mobility_bw_ref_noimg.mat"),
+    ],
+)
+def test_pair_mobility(Solver, periodicity, ref_file, tol):
+    zmax = 19.2
+    xymax = 76.8
+    params = self_mobility_params[Solver.__name__]
+
+    ref_dir = "./ref/"
+    ref = scipy.io.loadmat(ref_dir + ref_file)
+    refM = ref['M']
+    refHeights = ref['heights'].flatten()
+    nHeights = len(refHeights)
+
+    radH = 1.0 # hydrodynamic radius
+    eta = 1/4/np.sqrt(np.pi)
+
+    precision = np.float32 if Solver.precision == "float" else np.float64
+
+    solver = Solver(*periodicity)
+    solver.setParameters(**params)
+    numberParticles = 2
+    solver.initialize(
+        temperature=0,
+        viscosity=eta,
+        hydrodynamicRadius=radH,
+        numberParticles=numberParticles,
+    )
+
+    normMat = (1/(6*np.pi*eta))*np.ones((3*numberParticles, 3*numberParticles), dtype=precision)
+
+    seps = np.array([3 * radH, 4 * radH, 8 * radH])
+    nSeps = len(seps)
+
+    allM = np.zeros((nSeps, nHeights, 3*numberParticles, 3*numberParticles), dtype=precision)
+    for i in range(0,nSeps):
+        for k in range(0, nHeights):
+            xpos = xymax/2
+            positions = np.array([[xpos+seps[i]/2, xpos, refHeights[k]],
+                                  [xpos-seps[i]/2, xpos, refHeights[k]]], dtype=precision)
+            solver.setPositions(positions)
+
+            M = compute_M(solver, numberParticles)
+            M /= normMat
+            allM[i][k] = M
+
+    # uncomment to save datafile for test plots
+    # scipy.io.savemat('./temp/test_data/test_' + ref_file, {'M': allM, 'heights': refHeights})
+
+
+    indx, indy = 4, 1 ## xx
+    checkComponent(indx, indy, allM, refM, nSeps, tol)
+
+    indx, indy = 5, 2 # yy
+    checkComponent(indx, indy, allM, refM, nSeps, tol)
+
+    indx, indy = 6, 3 # zz
+    checkComponent(indx, indy, allM, refM, nSeps, tol)
+
+    indx, indy = 5, 1 # yx
+    checkComponent(indx, indy, allM, refM, nSeps, tol)
+
+    indx, indy = 3, 4 # zx
+    checkComponent(indx, indy, allM, refM, nSeps, tol)
+
+    indx, indy = 3, 5 # zy
+    checkComponent(indx, indy, allM, refM, nSeps, tol)
+
+def checkComponent(indx, indy, allM, refM, nSeps, tol):
+
+    indx -= 1 # shift from matlab to python indexing
+    indy -= 1
+    for i in range(0, nSeps):
+
+        xx = allM[i, :, indx, indy]
+        xx_ref = refM[i, :, indx, indy]
+
+        relDiff = np.abs([np.linalg.norm(xx - xx_ref)/np.linalg.norm(xx_ref + 1e-6) for xx, xx_ref in zip(xx, xx_ref)])
+        avgErr = np.mean(relDiff)
+
+        assert avgErr < tol, f"Pair mobility does not match reference for component {indx+1}, {indy+1}"