diff --git a/CMakeLists.txt b/CMakeLists.txt index c5e6c42c..c0e44224 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -46,7 +46,7 @@ endif() FetchContent_Declare( uammd GIT_REPOSITORY https://github.com/RaulPPelaez/uammd/ - GIT_TAG v2.11.0 + GIT_TAG v3.0.0 EXCLUDE_FROM_ALL ) FetchContent_MakeAvailable(uammd) diff --git a/solvers/DPStokes/extra/poly_fits.h b/solvers/DPStokes/extra/poly_fits.h index 0ea11a55..c4e18b33 100644 --- a/solvers/DPStokes/extra/poly_fits.h +++ b/solvers/DPStokes/extra/poly_fits.h @@ -37,9 +37,17 @@ double polyEval(std::vector polyCoeffs, double x) { } // Coefficients for the polynomial fit of c^{-1} from [1] -std::vector cbetam_inv = { +// the magic numbers come from the original DPStokes repo +// https://github.com/stochasticHydroTools/DoublyPeriodicStokes/blob/main/source/cpu/python/GridAndKernelConfig.py +std::vector cbeta_monopole_inv = { 4131643418.193291, -10471683395.26777, 11833009228.6429, -7851132955.882548, 3388121732.651829, -994285251.2185925, 201183449.7086889, -27776767.88241613, 2515647.646492857, -136305.2970161326, 3445.959503226691}; + +std::vector cbeta_dipole_inv = { + 39129605.95992675, -141095001.9380851, 226461386.9399433, + -213225975.7007014, 130589579.4920038, -54470021.75408537, + 15723274.21737542, -3119398.842468358, 411336.473729197, + -33239.91568304499, 1311.988715157221}; } // namespace dpstokes_polys \ No newline at end of file diff --git a/solvers/DPStokes/extra/uammd_interface.h b/solvers/DPStokes/extra/uammd_interface.h index cc502fa7..c7e00570 100644 --- a/solvers/DPStokes/extra/uammd_interface.h +++ b/solvers/DPStokes/extra/uammd_interface.h @@ -8,8 +8,10 @@ namespace uammd_dpstokes { // Instead of using uammd::real I have to re define real here. #ifndef DOUBLE_PRECISION using real = float; +using real3 = float3; #else using real = double; +using real3 = double3; #endif // This function returns either 'single' or 'double' according to the UAMMD's @@ -23,7 +25,6 @@ struct PyParameters { int nx = -1; int ny = -1; int nz = -1; - real dt = 0; real viscosity; real Lx; real Ly; @@ -32,8 +33,8 @@ struct PyParameters { real tolerance = 1e-5; real w, w_d; real hydrodynamicRadius = -1; - real beta = -1; - real beta_d = -1; + real3 beta = {-1.0, -1.0, -1.0}; + real3 beta_d = {-1.0, -1.0, -1.0}; real alpha = -1; real alpha_d = -1; // Can be either none, bottom, slit or periodic diff --git a/solvers/DPStokes/extra/uammd_wrapper.cu b/solvers/DPStokes/extra/uammd_wrapper.cu index e68b8717..b18cf075 100644 --- a/solvers/DPStokes/extra/uammd_wrapper.cu +++ b/solvers/DPStokes/extra/uammd_wrapper.cu @@ -53,10 +53,13 @@ auto createFCMParameters(PyParameters pypar) { par.tolerance = pypar.tolerance; par.box = uammd::Box({pypar.Lx, pypar.Ly, pypar.zmax - pypar.zmin}); par.cells = {pypar.nx, pypar.ny, pypar.nz}; - par.kernel = std::make_shared(pypar.w, pypar.alpha, pypar.beta, - pypar.Lx / pypar.nx); + par.kernel = + std::make_shared(pypar.w, pypar.alpha, + pypar.beta.x, // TODO beta parameter may need to + // be adjusted for non-square? + pypar.Lx / pypar.nx); par.kernelTorque = std::make_shared( - pypar.w_d, pypar.alpha_d, pypar.beta_d, pypar.Lx / pypar.nx); + pypar.w_d, pypar.alpha_d, pypar.beta_d.x, pypar.Lx / pypar.nx); return par; } @@ -76,14 +79,12 @@ auto createDPStokesParameters(PyParameters pypar) { par.nx = pypar.nx; par.ny = pypar.ny; par.nz = pypar.nz; - par.dt = pypar.dt; par.viscosity = pypar.viscosity; par.Lx = pypar.Lx; par.Ly = pypar.Ly; par.H = pypar.zmax - pypar.zmin; par.w = pypar.w; par.w_d = pypar.w_d; - par.hydrodynamicRadius = pypar.hydrodynamicRadius; par.beta = pypar.beta; par.beta_d = pypar.beta_d; par.alpha = pypar.alpha; @@ -98,8 +99,7 @@ private: auto computeHydrodynamicDisplacements(const auto *d_pos, const uammd::real4 *d_force, const uammd::real4 *d_torques, - int numberParticles, real dt, - real dtTorque, cudaStream_t st) { + int numberParticles, cudaStream_t st) { if (fcm) { return fcm->computeHydrodynamicDisplacements( (uammd::real4 *)(d_pos), (uammd::real4 *)(d_force), @@ -192,8 +192,7 @@ public: Real3ToReal4SubstractOriginZ(zOrigin)); auto mob = this->computeHydrodynamicDisplacements( stored_positions.data().get(), force4.data().get(), - (includeAngular) ? torque4.data().get() : nullptr, numberParticles, 0.0, - 0.0, st); + (includeAngular) ? torque4.data().get() : nullptr, numberParticles, st); // always copy linear translation but only angular if needed thrust::copy(thrust::cuda::par.on(st), mob.first.begin(), mob.first.end(), (uammd::real3 *)h_MF); diff --git a/solvers/DPStokes/mobility.h b/solvers/DPStokes/mobility.h index d98b191f..3384b236 100644 --- a/solvers/DPStokes/mobility.h +++ b/solvers/DPStokes/mobility.h @@ -50,10 +50,6 @@ class DPStokes : public libmobility::Mobility { void setParametersDPStokes(DPStokesParameters i_dppar) { this->dppar = i_dppar; - if (this->dppar.Lx != this->dppar.Ly) - throw std::runtime_error("[DPStokes] Only square periodic boxes (Lx = " - "Ly) are currently supported.\n"); - dpstokes = std::make_shared(); } @@ -71,33 +67,76 @@ class DPStokes : public libmobility::Mobility { if (ipar.includeAngular) { this->dppar.w = 6; this->dppar.w_d = 6; - this->dppar.beta = 1.327 * this->dppar.w; - this->dppar.beta_d = 2.217 * this->dppar.w; + this->dppar.beta = {real(1.327) * this->dppar.w, + real(1.327) * this->dppar.w, -1.0}; + this->dppar.beta_d = {real(2.217) * this->dppar.w_d, + real(2.217) * this->dppar.w_d, -1.0}; h = this->dppar.hydrodynamicRadius / 1.731; this->dppar.alpha_d = this->dppar.w_d * 0.5; } else { + // w=4 this->dppar.w = 4; - this->dppar.beta = 1.785 * this->dppar.w; + this->dppar.beta = {real(1.785) * this->dppar.w, + real(1.785) * this->dppar.w, -1.0}; h = this->dppar.hydrodynamicRadius / 1.205; + + // w=6 + // this->dppar.w = 6; + // this->dppar.beta = {real(1.714) * this->dppar.w, real(1.714) * + // this->dppar.w, -1.0}; + // h = this->dppar.hydrodynamicRadius / 1.554; } this->dppar.alpha = this->dppar.w * 0.5; this->dppar.tolerance = 1e-6; - int N = floor(this->dppar.Lx / h); - N += N % 2; + int Nx = floor(this->dppar.Lx / h); + Nx += Nx % 2; + int Ny = floor(this->dppar.Ly / h); + Ny += Ny % 2; - this->dppar.nx = N; - this->dppar.ny = N; + this->dppar.nx = Nx; + this->dppar.ny = Ny; - // note: this part is only configured for square boxes if (this->dppar.allowChangingBoxSize) { // adjust box size to suit h - this->dppar.Lx = N * h; - this->dppar.Ly = N * h; + this->dppar.Lx = Nx * h; + this->dppar.Ly = Ny * h; } else { // adjust h so that L/h is an integer - h = this->dppar.Lx / N; - double arg = this->dppar.hydrodynamicRadius / (this->dppar.w * h); - this->dppar.beta = - dpstokes_polys::polyEval(dpstokes_polys::cbetam_inv, arg); + real h_x = this->dppar.Lx / Nx; + real h_y = this->dppar.Ly / Ny; + h = min(h_x, h_y); + double arg = this->dppar.hydrodynamicRadius / (this->dppar.w * h_x); + real beta_x = + dpstokes_polys::polyEval(dpstokes_polys::cbeta_monopole_inv, arg); + arg = this->dppar.hydrodynamicRadius / (this->dppar.w * h_y); + real beta_y = + dpstokes_polys::polyEval(dpstokes_polys::cbeta_monopole_inv, arg); + + if (beta_x < 4.0 || beta_x > 18.0 || beta_y < 4.0 || beta_y > 18.0) { + throw std::runtime_error( + "[DPStokes] Could not find (h,beta) within interp range. This " + "means the particle radius and grid spacing are incompatible- try " + "a square domain."); + } + + this->dppar.beta = {beta_x, beta_y, min(beta_x, beta_y)}; + + if (ipar.includeAngular) { + arg = this->dppar.hydrodynamicRadius / (this->dppar.w_d * h_x); + real beta_xd = + dpstokes_polys::polyEval(dpstokes_polys::cbeta_dipole_inv, arg); + arg = this->dppar.hydrodynamicRadius / (this->dppar.w_d * h_y); + real beta_yd = + dpstokes_polys::polyEval(dpstokes_polys::cbeta_dipole_inv, arg); + if (beta_xd < 4.0 || beta_xd > 18.0 || beta_yd < 4.0 || + beta_yd > 18.0) { + throw std::runtime_error( + "[DPStokes] Could not find (h,beta) within interp range. This " + "means the particle radius and grid spacing are incompatible- " + "try " + "a square domain."); + } + this->dppar.beta_d = {beta_xd, beta_yd, min(beta_xd, beta_yd)}; + } } // Add a buffer of 1.5*w*h/2 when there is an open boundary @@ -113,7 +152,8 @@ class DPStokes : public libmobility::Mobility { // 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) has two factors of 2 and is FFT friendly + // pick nearby N such that 2(Nz-1) has two factors of 2 and is FFT + // friendly this->dppar.nz = floor(nz_actual); this->dppar.nz += (int)ceil(nz_actual) % 2; diff --git a/tests/ref/pair_mobility_bw_w6_offset_x.npz b/tests/ref/pair_mobility_bw_w6_offset_x.npz index 896be620..6b54215f 100644 Binary files a/tests/ref/pair_mobility_bw_w6_offset_x.npz and b/tests/ref/pair_mobility_bw_w6_offset_x.npz differ diff --git a/tests/ref/pair_mobility_bw_w6_offset_y.npz b/tests/ref/pair_mobility_bw_w6_offset_y.npz index 96bd5607..ce459dad 100644 Binary files a/tests/ref/pair_mobility_bw_w6_offset_y.npz and b/tests/ref/pair_mobility_bw_w6_offset_y.npz differ diff --git a/tests/ref/pair_mobility_sc_w6_offset_x.npz b/tests/ref/pair_mobility_sc_w6_offset_x.npz index 43c2c4ec..a1b61cc6 100644 Binary files a/tests/ref/pair_mobility_sc_w6_offset_x.npz and b/tests/ref/pair_mobility_sc_w6_offset_x.npz differ diff --git a/tests/ref/pair_mobility_sc_w6_offset_y.npz b/tests/ref/pair_mobility_sc_w6_offset_y.npz index cdf51b07..7f485c49 100644 Binary files a/tests/ref/pair_mobility_sc_w6_offset_y.npz and b/tests/ref/pair_mobility_sc_w6_offset_y.npz differ diff --git a/tests/ref/readme.txt b/tests/ref/readme.txt index 1edcfb64..c84e8dd8 100644 --- a/tests/ref/readme.txt +++ b/tests/ref/readme.txt @@ -2,6 +2,9 @@ Some of this data is generated from the original DoublyPeriodicStokes repository 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. +8/21/2025: Regenerated DPStokes data for angular tests since the procedure for selecting beta for torques was updated. + Agreed with previous data at tol=5e-2. + Abbreviations: bw- bottom wall sc- slit channel diff --git a/tests/ref/self_mobility_bw_w6.npz b/tests/ref/self_mobility_bw_w6.npz index 5f1d1552..6148a37d 100644 Binary files a/tests/ref/self_mobility_bw_w6.npz and b/tests/ref/self_mobility_bw_w6.npz differ diff --git a/tests/ref/self_mobility_sc_w6.npz b/tests/ref/self_mobility_sc_w6.npz index 9d6abffd..ccdc940c 100644 Binary files a/tests/ref/self_mobility_sc_w6.npz and b/tests/ref/self_mobility_sc_w6.npz differ diff --git a/tests/test_dpstokes.py b/tests/test_dpstokes.py new file mode 100644 index 00000000..9961440a --- /dev/null +++ b/tests/test_dpstokes.py @@ -0,0 +1,128 @@ +import numpy as np +import libMobility as lm + + +def compute_with_dpstokes(pos, lx, ly, forces=None, torques=None): + params = {"viscosity": 1 / (6 * np.pi), "hydrodynamicRadius": 1.73} + dpstokes = lm.DPStokes("periodic", "periodic", "single_wall") + torques_on = True if torques is not None else False + dpstokes.setParameters(Lx=lx, Ly=ly, zmin=-8.0, zmax=8.0) + dpstokes.initialize(**params, includeAngular=torques_on) + dpstokes.setPositions(pos) + mf_dp, mt_dp = dpstokes.Mdot(forces=forces, torques=torques) + if torques_on: + return mt_dp + else: + return mf_dp + + +def test_non_square_box(): + nP = 10 + pos_or = np.random.uniform(-8, 8, (nP, 3)).astype(np.float32) + forces_or = np.random.uniform(-1, 1, (nP, 3)).astype(np.float32) + forces_or -= np.mean(forces_or, axis=0) # Center the forces around zero + + L_fact = 2.0 + L_short = 16.0 + L_long = L_short * L_fact + + mf_dp_cube = compute_with_dpstokes(pos_or, forces=forces_or, lx=L_short, ly=L_short) + mt_dp_cube = compute_with_dpstokes( + pos_or, torques=forces_or, lx=L_short, ly=L_short + ) + + # This must be identical than doubling the box size in x and repeating the particles + pos = np.vstack((pos_or, pos_or + np.array([0.5 * L_long, 0.0, 0.0]))) + forces = np.vstack((forces_or, forces_or)) + mf_dp_lx = compute_with_dpstokes(pos, forces=forces, lx=L_long, ly=L_short)[ + : len(pos_or), : + ] + mt_dp_lx = compute_with_dpstokes(pos, torques=forces, lx=L_long, ly=L_short)[ + : len(pos_or), : + ] + + # And the same for doubling the box size in y + pos = np.vstack((pos_or, pos_or + np.array([0.0, 0.5 * L_long, 0.0]))) + forces = np.vstack((forces_or, forces_or)) + mf_dp_ly = compute_with_dpstokes(pos, forces=forces, lx=L_short, ly=L_long)[ + : len(pos_or), : + ] + mt_dp_ly = compute_with_dpstokes(pos, torques=forces, lx=L_short, ly=L_long)[ + : len(pos_or), : + ] + + assert np.allclose(mf_dp_lx, mf_dp_cube, rtol=1e-3, atol=1e-2) + assert np.allclose(mf_dp_ly, mf_dp_cube, rtol=1e-3, atol=1e-2) + + assert np.allclose(mt_dp_lx, mt_dp_cube, rtol=1e-3, atol=1e-2) + assert np.allclose(mt_dp_ly, mt_dp_cube, rtol=1e-3, atol=1e-2) + + +def test_isotropy(): + L0 = 16.0 + L_fact = 2.0 + + f_x = np.array([1.0, 0.0, 0.0]) + f_y = np.array([0.0, 1.0, 0.0]) + + n_runs = 5 + for _ in range(n_runs): + pos = np.random.uniform(-8, 8, (3)).astype(np.float32) + + mf_x = compute_with_dpstokes(pos, forces=f_x, lx=L0 * L_fact, ly=L0) + mf_y = compute_with_dpstokes(pos, forces=f_y, lx=L0, ly=L0 * L_fact) + + mt_x = compute_with_dpstokes(pos, torques=f_x, lx=L0 * L_fact, ly=L0) + mt_y = compute_with_dpstokes(pos, torques=f_y, lx=L0, ly=L0 * L_fact) + + assert np.allclose(mf_x[0], mf_y[1], rtol=1e-3, atol=1e-2) + assert np.allclose(mf_x[1], mf_y[0], rtol=1e-3, atol=1e-2) + assert np.allclose(mf_x[2], mf_y[2], rtol=1e-3, atol=1e-2) + + assert np.allclose(mt_x[0], mt_y[1], rtol=1e-3, atol=1e-2) + assert np.allclose(mt_x[1], mt_y[0], rtol=1e-3, atol=1e-2) + assert np.allclose(mt_x[2], mt_y[2], rtol=1e-3, atol=1e-2) + + +def test_non_square_matching_rpy(): + + a = np.random.uniform(0.5, 2.0) + eta = 1 / (6 * np.pi * a) + L_min = 100 * a # need a fairly large domain to neglect periodic effects + Lx = np.random.uniform(L_min, 3.0 * L_min) + Ly = np.random.uniform(L_min, 3.0 * L_min) + max_height = 10.0 * a + min_height = 1.5 * a # need a buffer for regularized kernels to match + + solver = lm.DPStokes("periodic", "periodic", "single_wall") + solver.setParameters(Lx=Lx, Ly=Ly, zmin=0.0, zmax=max_height) + solver.initialize(hydrodynamicRadius=a, viscosity=eta, includeAngular=True) + + rpy = lm.NBody("open", "open", "single_wall") + rpy.setParameters(wallHeight=0.0) + rpy.initialize(hydrodynamicRadius=a, viscosity=eta, includeAngular=True) + + F = np.eye(6) + + heights = np.linspace(min_height, max_height, 20) + + results_dpstokes = np.zeros((2 * len(heights), 3)) + results_rpy = np.zeros((2 * len(heights), 3)) + + for i, h in enumerate(heights): + pos = np.array([0.0, 0.0, h]) + + rpy.setPositions(pos) + solver.setPositions(pos) + for j in range(6): + f_j = F[0:3, j].copy() + t_j = F[3:6, j].copy() + mf_j, mt_j = solver.Mdot(forces=f_j, torques=t_j) + results_dpstokes[i] += mf_j + results_dpstokes[i + 3] += mt_j + + mf_j, mt_j = rpy.Mdot(forces=f_j, torques=t_j) + results_rpy[i] += mf_j + results_rpy[i + 3] += mt_j + + assert np.allclose(results_dpstokes, results_rpy, rtol=1e-3, atol=1e-2) diff --git a/tests/test_initialization.py b/tests/test_initialization.py index 9700de07..a001ed67 100644 --- a/tests/test_initialization.py +++ b/tests/test_initialization.py @@ -32,14 +32,6 @@ def test_invalid_throws(Solver): solver = Solver("periodicasdas", "periodic", "open") -def test_dpstokes_invalid_box(): - with pytest.raises(RuntimeError): - Lx, Ly = 10, 20 - solver = DPStokes("periodic", "periodic", "single_wall") - params = {"Lx": Lx, "Ly": Ly, "zmin": 0, "zmax": 19.2} - solver.setParameters(**params) - - @pytest.mark.parametrize( ("NBatch", "NperBatch"), [ diff --git a/tests/test_wall_mobility.py b/tests/test_wall_mobility.py index 569cf05f..67786785 100644 --- a/tests/test_wall_mobility.py +++ b/tests/test_wall_mobility.py @@ -2,48 +2,26 @@ import numpy as np from libMobility import DPStokes, NBody -import libMobility from utils import compute_M, get_wall_params import os # NOTE: Some of the following tests will only pass if compiled with double precision. # This is because the reference data was generated in double precision. -precision = np.float32 if NBody.precision == "float" else np.float64 - +precision_str = "single" if NBody.precision == "float" else "double" ref_dir = os.path.dirname(os.path.abspath(__file__)) + "/ref/" @pytest.mark.parametrize( - ("Solver", "periodicity", "tol", "ref_file"), + ("Solver", "periodicity", "ref_file"), [ - ( - DPStokes, - ("periodic", "periodic", "single_wall"), - 1e-6, - "self_mobility_bw_w4.npz", - ), - ( - DPStokes, - ("periodic", "periodic", "two_walls"), - 1e-6, - "self_mobility_sc_w4.npz", - ), - ( - NBody, - ("open", "open", "single_wall"), - 1e-6, - "self_mobility_bw_ref_noimg.npz", - ), + (DPStokes, ("periodic", "periodic", "single_wall"), "self_mobility_bw_w4.npz"), + (DPStokes, ("periodic", "periodic", "two_walls"), "self_mobility_sc_w4.npz"), + (NBody, ("open", "open", "single_wall"), "self_mobility_bw_ref_noimg.npz"), ], ) @pytest.mark.parametrize("wallHeight", [0, 5.4, -10]) -def test_self_mobility_linear(Solver, periodicity, tol, ref_file, wallHeight): - if precision == np.float32 and Solver.__name__ == "DPStokes": - pytest.skip( - "The test is only valid for double precision due to how reference data was generated." - ) - +def test_self_mobility_linear(Solver, periodicity, ref_file, wallHeight): xymax = 76.8 params = get_wall_params(Solver.__name__, wallHeight) includeAngular = False @@ -54,6 +32,8 @@ def test_self_mobility_linear(Solver, periodicity, tol, ref_file, wallHeight): hydrodynamicRadius = 1.0 eta = 1 / 4 / np.sqrt(np.pi) + tol = 100 * np.finfo(precision_str).eps + solver = Solver(*periodicity) solver.setParameters(**params) numberParticles = 1 @@ -64,18 +44,18 @@ def test_self_mobility_linear(Solver, periodicity, tol, ref_file, wallHeight): nHeights = len(refHeights) - normMat = np.ones((3 * numberParticles, 3 * numberParticles), dtype=precision) + normMat = np.ones((3 * numberParticles, 3 * numberParticles), dtype=precision_str) 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 + (nHeights, 3 * numberParticles, 3 * numberParticles), dtype=precision_str ) for i in range(0, nHeights): positions = np.array( - [[xymax / 2, xymax / 2, refHeights[i] + wallHeight]], dtype=precision + [[xymax / 2, xymax / 2, refHeights[i] + wallHeight]], dtype=precision_str ) solver.setPositions(positions) @@ -88,42 +68,23 @@ def test_self_mobility_linear(Solver, periodicity, tol, ref_file, wallHeight): np.diag(matrix)[0:3] for matrix in refM ] # only take diagonal elements from forces - for diag, ref_diag in zip(diags, ref_diags): - assert np.allclose( - diag, ref_diag, rtol=tol, atol=tol - ), f"Self mobility does not match reference" + assert np.allclose(diags, ref_diags, rtol=tol, atol=tol) @pytest.mark.parametrize( - ("Solver", "periodicity", "tol", "ref_file"), + ("Solver", "periodicity", "ref_file"), [ - ( - DPStokes, - ("periodic", "periodic", "single_wall"), - 1e-6, - "pair_mobility_bw_w4.npz", - ), - ( - DPStokes, - ("periodic", "periodic", "two_walls"), - 1e-6, - "pair_mobility_sc_w4.npz", - ), + (DPStokes, ("periodic", "periodic", "single_wall"), "pair_mobility_bw_w4.npz"), + (DPStokes, ("periodic", "periodic", "two_walls"), "pair_mobility_sc_w4.npz"), ( NBody, ("open", "open", "single_wall"), - 1e-4, "pair_mobility_bw_ref_noimg_offset_x.npz", ), ], ) @pytest.mark.parametrize("wallHeight", [0, 5.4, -10]) -def test_pair_mobility_linear(Solver, periodicity, ref_file, tol, wallHeight): - if precision == np.float32 and Solver.__name__ == "DPStokes": - pytest.skip( - "The test is only valid for double precision due to how reference data was generated." - ) - +def test_pair_mobility_linear(Solver, periodicity, ref_file, wallHeight): xymax = 76.8 params = get_wall_params(Solver.__name__, wallHeight) includeAngular = False @@ -135,7 +96,7 @@ def test_pair_mobility_linear(Solver, periodicity, ref_file, tol, wallHeight): radH = 1.0 # hydrodynamic radius eta = 1 / 4 / np.sqrt(np.pi) - tol = 100 * np.finfo(precision).eps + tol = 100 * np.finfo(precision_str).eps solver = Solver(*periodicity) solver.setParameters(**params) @@ -146,14 +107,14 @@ def test_pair_mobility_linear(Solver, periodicity, ref_file, tol, wallHeight): ) normMat = (1 / (6 * np.pi * eta)) * np.ones( - (3 * numberParticles, 3 * numberParticles), dtype=precision + (3 * numberParticles, 3 * numberParticles), dtype=precision_str ) seps = np.array([3 * radH, 4 * radH, 8 * radH]) nSeps = len(seps) allM = np.zeros( - (nSeps, nHeights, 3 * numberParticles, 3 * numberParticles), dtype=precision + (nSeps, nHeights, 3 * numberParticles, 3 * numberParticles), dtype=precision_str ) for i in range(0, nSeps): for k in range(0, nHeights): @@ -163,7 +124,7 @@ def test_pair_mobility_linear(Solver, periodicity, ref_file, tol, wallHeight): [xpos + seps[i] / 2, xpos, refHeights[k] + wallHeight], [xpos - seps[i] / 2, xpos, refHeights[k] + wallHeight], ], - dtype=precision, + dtype=precision_str, ) solver.setPositions(positions) @@ -171,10 +132,7 @@ def test_pair_mobility_linear(Solver, periodicity, ref_file, tol, wallHeight): M /= normMat allM[i][k] = M - for i in range(0, nSeps): - for k in range(0, nHeights): - diff = abs(allM[i, k] - refM[i, k][0:6, 0:6]) - assert np.all(diff < tol) + assert np.allclose(allM, refM[:, :, 0:6, 0:6], atol=tol, rtol=tol) @pytest.mark.parametrize( @@ -187,7 +145,7 @@ def test_pair_mobility_linear(Solver, periodicity, ref_file, tol, wallHeight): ) @pytest.mark.parametrize("wallHeight", [0, 5.4, -10]) def test_self_mobility_angular(Solver, periodicity, ref_file, wallHeight): - if precision == np.float32 and Solver.__name__ == "DPStokes": + if precision_str == "single" and Solver.__name__ == "DPStokes": pytest.skip( "The test is only valid for double precision due to how reference data was generated." ) @@ -198,8 +156,9 @@ def test_self_mobility_angular(Solver, periodicity, ref_file, wallHeight): hydrodynamicRadius = 1.0 eta = 1 / 4 / np.sqrt(np.pi) + tol = 100 * np.finfo(precision_str).eps + includeAngular = True - tol = 1e-6 ref = np.load(ref_dir + ref_file) refM = ref["M"] @@ -216,18 +175,18 @@ def test_self_mobility_angular(Solver, periodicity, ref_file, wallHeight): nHeights = len(refHeights) - normMat = np.zeros((6 * numberParticles, 6 * numberParticles), dtype=precision) + normMat = np.zeros((6 * numberParticles, 6 * numberParticles), dtype=precision_str) normMat[0:3, 0:3] = 1 / (6 * np.pi * eta * hydrodynamicRadius) # tt normMat[3:, 3:] = 1 / (8 * np.pi * eta * hydrodynamicRadius**3) # rr normMat[3:, 0:3] = 1 / (6 * np.pi * eta * hydrodynamicRadius**2) # tr normMat[0:3, 3:] = 1 / (6 * np.pi * eta * hydrodynamicRadius**2) # rt allM = np.zeros( - (nHeights, 6 * numberParticles, 6 * numberParticles), dtype=precision + (nHeights, 6 * numberParticles, 6 * numberParticles), dtype=precision_str ) for i in range(0, nHeights): positions = np.array( - [[xymax / 2, xymax / 2, refHeights[i] + wallHeight]], dtype=precision + [[xymax / 2, xymax / 2, refHeights[i] + wallHeight]], dtype=precision_str ) solver.setPositions(positions) @@ -235,9 +194,7 @@ def test_self_mobility_angular(Solver, periodicity, ref_file, wallHeight): M /= normMat allM[i] = M - for i in range(0, nHeights): - diff = abs(allM[i] - refM[i]) - assert np.all(diff < tol) + assert np.allclose(allM, refM, atol=tol, rtol=tol) @pytest.mark.parametrize( @@ -251,7 +208,7 @@ def test_self_mobility_angular(Solver, periodicity, ref_file, wallHeight): @pytest.mark.parametrize("offset", ["x", "y"]) @pytest.mark.parametrize("wallHeight", [0, 5.4, -10]) def test_pair_mobility_angular(Solver, periodicity, ref_file, offset, wallHeight): - if precision == np.float32 and Solver.__name__ == "DPStokes": + if precision_str == "single" and Solver.__name__ == "DPStokes": pytest.skip( "The test is only valid for double precision due to how reference data was generated." ) @@ -262,7 +219,7 @@ def test_pair_mobility_angular(Solver, periodicity, ref_file, offset, wallHeight eta = 1 / 4 / np.sqrt(np.pi) includeAngular = True - tol = 1e-6 + tol = 100 * np.finfo(precision_str).eps ref_file += "_offset_" + offset + ".npz" ref = np.load(ref_dir + ref_file) @@ -285,7 +242,7 @@ def test_pair_mobility_angular(Solver, periodicity, ref_file, offset, wallHeight ) nSeps = len(seps) - normMat = np.zeros((6 * nP, 6 * nP), dtype=precision) + normMat = np.zeros((6 * nP, 6 * nP), dtype=precision_str) normMat[0 : 3 * nP, 0 : 3 * nP] = 1 / (6 * np.pi * eta * hydrodynamicRadius) # tt normMat[3 * nP :, 3 * nP :] = 1 / (8 * np.pi * eta * hydrodynamicRadius**3) # rr normMat[3 * nP :, 0 : 3 * nP] = 1 / (6 * np.pi * eta * hydrodynamicRadius**2) # tr @@ -299,7 +256,7 @@ def test_pair_mobility_angular(Solver, periodicity, ref_file, offset, wallHeight raise ValueError("Test for offset in {} not implemented".format(offset)) xpos = xymax / 2 - allM = np.zeros((nSeps, nHeights, 6 * nP, 6 * nP), dtype=precision) + allM = np.zeros((nSeps, nHeights, 6 * nP, 6 * nP), dtype=precision_str) for i in range(0, nSeps): seps_vec = (seps[i] * offset_vec) / 2 for k in range(0, nHeights): @@ -308,7 +265,7 @@ def test_pair_mobility_angular(Solver, periodicity, ref_file, offset, wallHeight [xpos, xpos, refHeights[k] + wallHeight] + seps_vec, [xpos, xpos, refHeights[k] + wallHeight] - seps_vec, ], - dtype=precision, + dtype=precision_str, ) solver.setPositions(positions) @@ -316,7 +273,48 @@ def test_pair_mobility_angular(Solver, periodicity, ref_file, offset, wallHeight M /= normMat allM[i, k] = M - for i in range(0, nSeps): - for k in range(0, nHeights): - diff = abs(allM[i, k] - refM[i, k]) - assert np.all(diff < tol) + assert np.allclose(allM, refM, atol=tol, rtol=tol) + + +def test_dpstokes_matching_rpy(): + + a = np.random.uniform(0.5, 2.0) + eta = 1 / (6 * np.pi * a) + L_min = 100 * a # need a fairly large domain to neglect periodic effects + L_fact = np.random.uniform(1.0, 3.0) + L = L_min * L_fact + max_height = 10.0 * a + min_height = 2.0 * a # need a buffer for regularized kernels to match + + solver = DPStokes("periodic", "periodic", "single_wall") + solver.setParameters(Lx=L, Ly=L, zmin=0.0, zmax=max_height) + solver.initialize(hydrodynamicRadius=a, viscosity=eta, includeAngular=True) + + rpy = NBody("open", "open", "single_wall") + rpy.setParameters(wallHeight=0.0) + rpy.initialize(hydrodynamicRadius=a, viscosity=eta, includeAngular=True) + + nP = 10 + # generate particles in interior of domain to avoid periodic artifacts + pos = np.random.uniform(L / 4, 3 * L / 4, (nP, 3)).astype(np.float32) + pos_z = np.random.uniform(min_height, max_height, nP).astype(np.float32) + pos[:, 2] = pos_z + F = np.eye(6 * nP) + + results_dpstokes = np.zeros((2 * 6 * nP, 3 * nP)) + results_rpy = np.zeros((2 * 6 * nP, 3 * nP)) + + rpy.setPositions(pos) + solver.setPositions(pos) + for j in range(6 * nP): + f_j = F[: 3 * nP, j].copy() + t_j = F[3 * nP :, j].copy() + mf_j, mt_j = solver.Mdot(forces=f_j, torques=t_j) + results_dpstokes[j] += mf_j + results_dpstokes[j + 6 * nP] += mt_j + + mf_j, mt_j = rpy.Mdot(forces=f_j, torques=t_j) + results_rpy[j] += mf_j + results_rpy[j + 6 * nP] += mt_j + + assert np.allclose(results_dpstokes, results_rpy, rtol=1e-3, atol=1e-2)