diff --git a/include/MobilityInterface/random_finite_differences.h b/include/MobilityInterface/random_finite_differences.h index 22b580a..4722661 100644 --- a/include/MobilityInterface/random_finite_differences.h +++ b/include/MobilityInterface/random_finite_differences.h @@ -31,8 +31,7 @@ template void random_finite_differences(Mdot mdot, device_span positions, device_span ilinear, device_span iangular, uint seed, - real prefactor = 1) { - constexpr real delta = 1e-3; + const real delta, real prefactor = 1) { const uint N = ilinear.size() / 3; using device_vector = thrust::device_vector>; diff --git a/solvers/DPStokes/extra/uammd_interface.h b/solvers/DPStokes/extra/uammd_interface.h index c7e0057..1139320 100644 --- a/solvers/DPStokes/extra/uammd_interface.h +++ b/solvers/DPStokes/extra/uammd_interface.h @@ -31,6 +31,7 @@ struct PyParameters { real zmin, zmax; // Tolerance will be ignored in DP mode, TP will use only tolerance and nxy/nz real tolerance = 1e-5; + real delta = 1e-3; // RFD step size real w, w_d; real hydrodynamicRadius = -1; real3 beta = {-1.0, -1.0, -1.0}; diff --git a/solvers/DPStokes/mobility.h b/solvers/DPStokes/mobility.h index 7da91fc..876179e 100644 --- a/solvers/DPStokes/mobility.h +++ b/solvers/DPStokes/mobility.h @@ -222,8 +222,9 @@ class DPStokes : public libmobility::Mobility { uint seed = rng(); device_vector thermal_drift_m(ilinear.size(), 0); device_vector thermal_drift_d(iangular.size(), 0); - libmobility::random_finite_differences(mdot, original_pos, thermal_drift_m, - thermal_drift_d, seed, prefactor); + libmobility::random_finite_differences( + mdot, original_pos, thermal_drift_m, thermal_drift_d, seed, + this->dppar.delta * this->dppar.hydrodynamicRadius, prefactor); this->setPositions(original_pos); thrust::transform(thrust::cuda::par, thermal_drift_m.begin(), thermal_drift_m.end(), linear.begin(), linear.begin(), diff --git a/solvers/DPStokes/python_wrapper.cu b/solvers/DPStokes/python_wrapper.cu index 4cf4e59..402e118 100644 --- a/solvers/DPStokes/python_wrapper.cu +++ b/solvers/DPStokes/python_wrapper.cu @@ -22,6 +22,8 @@ zmax : float The maximum value of the z coordinate. This is the position of the top wall if the Z periodicity is `two_walls`. allowChangingBoxSize : bool Whether the periodic extents Lx & Ly can be modified during parameter selection. Default: false. +delta : float + The finite difference step size for random finite differences. Specified in units of hydrodynamicRadius. Default is 1e-3. )pbdoc"; static const char *docstring = R"pbdoc( @@ -35,18 +37,20 @@ The periodicity must be set to `periodic` in the X and Y directions. The Z perio )pbdoc"; MOBILITY_PYTHONIFY_WITH_EXTRA_CODE( - DPStokes, solver.def( - "setParameters", - [](DPStokes &self, real Lx, real Ly, real zmin, real zmax, - bool allowChangingBoxSize) { - DPStokesParameters params; - params.Lx = Lx; - params.Ly = Ly; - params.zmin = zmin; - params.zmax = zmax; - params.allowChangingBoxSize = allowChangingBoxSize; - self.setParametersDPStokes(params); - }, - "Lx"_a, "Ly"_a, "zmin"_a, "zmax"_a, - "allowChangingBoxSize"_a = false, setparameters_docstring); + DPStokes, + solver.def( + "setParameters", + [](DPStokes &self, real Lx, real Ly, real zmin, real zmax, + bool allowChangingBoxSize, real delta) { + DPStokesParameters params; + params.Lx = Lx; + params.Ly = Ly; + params.zmin = zmin; + params.zmax = zmax; + params.allowChangingBoxSize = allowChangingBoxSize; + params.delta = delta; + self.setParametersDPStokes(params); + }, + "Lx"_a, "Ly"_a, "zmin"_a, "zmax"_a, "allowChangingBoxSize"_a = false, + "delta"_a = 1e-3, setparameters_docstring); , docstring); diff --git a/solvers/NBody/mobility.h b/solvers/NBody/mobility.h index 6ead265..13cdeb2 100644 --- a/solvers/NBody/mobility.h +++ b/solvers/NBody/mobility.h @@ -28,6 +28,7 @@ class NBody : public libmobility::Mobility { nbody_rpy::algorithm algorithm = nbody_rpy::algorithm::advise; real wallHeight; // location of the wall in z + real delta; // finite difference step size for random finite differences // Batched functionality configuration int Nbatch = -1; @@ -53,6 +54,7 @@ class NBody : public libmobility::Mobility { int Nbatch = -1; int NperBatch = -1; std::optional wallHeight = std::nullopt; + real delta = 1e-3; // in units of particle radius }; /** * @brief Sets the parameters for the N-body computation @@ -94,6 +96,7 @@ class NBody : public libmobility::Mobility { "you want to use a wall, set periodicityZ to single_wall in the " "configuration."); } + this->delta = par.delta; } void initialize(Parameters ipar) override { @@ -195,8 +198,9 @@ class NBody : public libmobility::Mobility { device_vector thermal_drift_m(ilinear.size(), 0); device_vector thermal_drift_d(iangular.size(), 0); - libmobility::random_finite_differences(mdot, original_pos, thermal_drift_m, - thermal_drift_d, seed, prefactor); + libmobility::random_finite_differences( + mdot, original_pos, thermal_drift_m, thermal_drift_d, seed, + this->delta * this->hydrodynamicRadius, prefactor); device_adapter linear(ilinear, device::cuda); this->setPositions(original_pos); thrust::transform(thrust::cuda::par, thermal_drift_m.begin(), diff --git a/solvers/NBody/python_wrapper.cu b/solvers/NBody/python_wrapper.cu index add597d..7d6c31f 100644 --- a/solvers/NBody/python_wrapper.cu +++ b/solvers/NBody/python_wrapper.cu @@ -5,19 +5,21 @@ #include static const char *docstringSetParameters = R"pbdoc( - Set the parameters for the NBody solver. +Set the parameters for the NBody solver. - Parameters - ---------- - algorithm : str - The algorithm to use. Options are "naive", "fast", "block" and "advise". Default is "advise". - NBatch : int - The number of batches to use. If -1 (default), the number of batches is automatically determined. - NperBatch : int - The number of particles per batch. If -1 (default), the number of particles per batch is automatically determined. - wallHeight : float - The height of the wall. Only valid if periodicityZ is single_wall. - )pbdoc"; +Parameters +---------- +algorithm : str + The algorithm to use. Options are "naive", "fast", "block" and "advise". Default is "advise". +NBatch : int + The number of batches to use. If -1 (default), the number of batches is automatically determined. +NperBatch : int + The number of particles per batch. If -1 (default), the number of particles per batch is automatically determined. +wallHeight : float + The height of the wall. Only valid if periodicityZ is single_wall. +delta : float + The finite difference step size for random finite differences. Specified in units of hydrodynamicRadius. Default is 1e-3. +)pbdoc"; static const char *docstring = R"pbdoc( This module computes hydrodynamic interactions using an :math:`O(N^2)` algorithm. @@ -50,10 +52,10 @@ MOBILITY_PYTHONIFY_WITH_EXTRA_CODE( solver.def( "setParameters", [](NBody &myself, std::string algo, int NBatch, int NperBatch, - std::optional wallHeight) { + std::optional wallHeight, real delta) { myself.setParametersNBody({nbody_rpy::string2NBodyAlgorithm(algo), - NBatch, NperBatch, wallHeight}); + NBatch, NperBatch, wallHeight, delta}); }, docstringSetParameters, "algorithm"_a = "advise", "Nbatch"_a = -1, - "NperBatch"_a = -1, "wallHeight"_a = std::nullopt); + "NperBatch"_a = -1, "wallHeight"_a = std::nullopt, "delta"_a = 1e-3); , docstring); diff --git a/tests/test_thermal_drift.py b/tests/test_thermal_drift.py index 16dcdfa..bde3e27 100644 --- a/tests/test_thermal_drift.py +++ b/tests/test_thermal_drift.py @@ -1,45 +1,45 @@ import pytest import numpy as np -from libMobility import SelfMobility, PSE, NBody, DPStokes from utils import ( solver_configs_all, - solver_configs_torques, get_sane_params, generate_positions_in_box, ) -def average(function, num_averages): - result_m, result_d = function() - for i in range(num_averages - 1): - new_result_m, new_result_d = function() - result_m += new_result_m - if result_d is not None: - result_d += new_result_d - if result_d is not None: - result_d /= num_averages - return result_m / num_averages, result_d - - -def divM_rfd(solver, positions): - # RFD works by approximating \partial_q \dot M = 1/\delta \langle M(q+\delta/2 W)W - M(q-\delta/2 W)W \rangle +@pytest.mark.parametrize(("Solver", "periodicity"), solver_configs_all) +@pytest.mark.parametrize("includeAngular", [True, False]) +@pytest.mark.parametrize("a", [0.4, 1.0, 2.5]) +def test_deterministic_divM_matches_rfd(Solver, periodicity, a, includeAngular): + if Solver.__name__ == "PSE" and includeAngular: + pytest.skip("PSE does not support torques") + N_iter = 250001 if Solver.__name__ == "NBody" else 25001 + tol = 1e-2 if Solver.__name__ == "DPStokes" else 1e-3 + solver = Solver(*periodicity) + parameters = get_sane_params(Solver.__name__, periodicity[2]) delta = 1e-3 + if "delta" in parameters: + parameters["delta"] = delta + solver.setParameters(**parameters) + precision = np.float32 if solver.precision == "float" else np.float64 + solver.initialize( + viscosity=1.0, + hydrodynamicRadius=a, + includeAngular=includeAngular, + ) + numberParticles = 10 + positions = generate_positions_in_box(parameters, numberParticles).astype(precision) + solver.setPositions(positions) - def rfd_func(): - W = np.random.normal(size=positions.shape).astype(positions.dtype) - solver.setPositions(positions + delta / 2 * W) - _tdriftp_m, _tdriftp_d = solver.Mdot(W) - solver.setPositions(positions - delta / 2 * W) - _tdriftm_m, _tdriftm_d = solver.Mdot(W) - _tdrift_m = (_tdriftp_m - _tdriftm_m) / delta - _tdrift_d = ( - (_tdriftp_d - _tdriftm_d) / delta if _tdriftm_d is not None else None - ) - return _tdrift_m, _tdrift_d + # use delta*a in deterministic since libMobility solvers expect delta to be in units of a + det_div_m = deterministic_div_m(solver, positions, includeAngular, delta * a) + rfd_m, rfd_d = average_until_error_tolerance( + solver.divM, N_iter, det_div_m, numberParticles, includeAngular, tol + ) - solver.setPositions(positions) - tdrift_m, tdrift_d = average(lambda: average(rfd_func, 400), 100) - return tdrift_m, tdrift_d + assert np.allclose(det_div_m[0 : 3 * numberParticles], rfd_m, atol=tol, rtol=tol) + if includeAngular: + assert np.allclose(det_div_m[3 * numberParticles :], rfd_d, atol=tol, rtol=tol) @pytest.mark.parametrize(("Solver", "periodicity"), solver_configs_all) @@ -70,46 +70,6 @@ def test_divM_does_not_change_positions(Solver, periodicity): ), f"Mdot has changed: {np.max(np.abs(mf - mf2))}" -@pytest.mark.parametrize(("Solver", "periodicity"), solver_configs_all) -@pytest.mark.parametrize("hydrodynamicRadius", [1.0, 0.95, 1.12]) -@pytest.mark.parametrize("numberParticles", [1, 2, 3, 10]) -@pytest.mark.parametrize("includeAngular", [False, True]) -def test_divM_is_zero( - Solver, periodicity, hydrodynamicRadius, numberParticles, includeAngular -): - if not np.all(np.array(periodicity) == "open") and not np.all( - np.array(periodicity) == "periodic" - ): - pytest.skip("Only periodic and open boundary conditions have zero divergence") - if Solver.__name__ == "PSE" and includeAngular: - pytest.skip("PSE does not support torques") - precision = np.float32 if Solver.precision == "float" else np.float64 - solver = Solver(*periodicity) - parameters = get_sane_params(Solver.__name__, periodicity[2]) - solver.setParameters(**parameters) - solver.initialize( - viscosity=1.0, - hydrodynamicRadius=hydrodynamicRadius, - includeAngular=includeAngular, - ) - positions = generate_positions_in_box(parameters, numberParticles).astype(precision) - solver.setPositions(positions) - thermal_drift_m, thermal_drift_d = average(solver.divM, 3000) - assert np.allclose( - np.abs(thermal_drift_m), - 0.0, - atol=1e-5, - rtol=0, - ), f"Linear RFD drift is not zero: {np.max(np.abs(thermal_drift_m))}" - if thermal_drift_d is not None: - assert np.allclose( - np.abs(thermal_drift_d), - 0.0, - atol=1e-5, - rtol=0, - ), f"Dipolar RFD drift is not zero: {np.max(np.abs(thermal_drift_d))}" - - @pytest.mark.parametrize(("Solver", "periodicity"), solver_configs_all) @pytest.mark.parametrize("hydrodynamicRadius", [1.0, 0.95, 1.12]) @pytest.mark.parametrize("numberParticles", [1, 2, 3, 10]) @@ -142,42 +102,106 @@ def test_divM_returns_different_numbers( ), f"RFD is not different: {np.max(np.abs(rfd1 - rfd2))}" -@pytest.mark.parametrize(("Solver", "periodicity"), solver_configs_all) -@pytest.mark.parametrize("hydrodynamicRadius", [0.95]) -@pytest.mark.parametrize("numberParticles", [1, 10]) -@pytest.mark.parametrize("includeAngular", [False, True]) -def test_divM_matches_rfd( - Solver, periodicity, hydrodynamicRadius, numberParticles, includeAngular -): - if Solver.__name__ == "PSE" and includeAngular: - pytest.skip("PSE does not support torques") - precision = np.float32 if Solver.precision == "float" else np.float64 - solver = Solver(*periodicity) - parameters = get_sane_params(Solver.__name__, periodicity[2]) - solver.setParameters(**parameters) - solver.initialize( - viscosity=1.0, - hydrodynamicRadius=hydrodynamicRadius, - includeAngular=includeAngular, - ) - positions = np.asarray( - generate_positions_in_box(parameters, numberParticles).astype(precision) - ) - solver.setPositions(positions) - reference_m, reference_d = divM_rfd(solver, positions) +def average_until_error_tolerance(function, num_averages, soln, N, includeAngular, tol): + # track error of the running average of RFDs after each new sample. exit early if the error is below the tolerance + avg_err_m = np.zeros(num_averages) + rfd_m = np.zeros(3 * N) + avg_err_d = np.zeros(num_averages) + rfd_d = np.zeros(3 * N) - solver.setPositions(positions) - rfd_m, rfd_d = average(lambda: average(solver.divM, 400), 100) - assert np.allclose( - reference_m, - rfd_m, - atol=1e-3, - rtol=1e-3, - ), f"Linear RFD does not match: {np.max(np.abs(reference_m - rfd_m))}" - if reference_d is not None and rfd_d is not None: - assert np.allclose( - reference_d, - rfd_d, - atol=1e-3, - rtol=1e-3, - ), f"Dipole RFD does not match: {np.max(np.abs(reference_d - rfd_d))}" + soln_m = soln[0 : 3 * N] + soln_d = soln[3 * N :] if includeAngular else None + + err_check_iter = 0 + check_every = 5000 + exit_early = False + last_iter = -1 + for i in range(num_averages): + new_result_m, new_result_d = function() + rfd_m += new_result_m + running_m = rfd_m / float(i + 1) + avg_err_m[i] = np.linalg.norm(running_m - soln_m) + if new_result_d is not None: + rfd_d += new_result_d + running_d = rfd_d / float(i + 1) + avg_err_d[i] = np.linalg.norm((running_d - soln_d)) + + if err_check_iter == check_every: + err_m = np.allclose(running_m, soln_m, atol=tol, rtol=tol) + errs = [err_m] + if includeAngular: + err_d = np.allclose(running_d, soln_d, atol=tol, rtol=tol) + errs.append(err_d) + if all(errs): + last_iter = i + exit_early = True + break + err_check_iter = 0 + err_check_iter += 1 + + if exit_early: + num_averages = last_iter + 1 + + rfd_m /= float(num_averages) + rfd_d /= float(num_averages) + + avg_err_m = avg_err_m[:num_averages] + avg_err_d = avg_err_d[:num_averages] + + # plot_error(avg_err_m, avg_err_d, num_averages, includeAngular, delta, a) + + return rfd_m, rfd_d + + +def plot_error(avg_err_m, avg_err_d, num_averages, includeAngular): + import matplotlib.pyplot as plt + + plt.figure() + n = np.arange(1, num_averages + 1) + plt.loglog(n, avg_err_m, marker="o", linestyle="-") + plt.loglog(n, avg_err_d, marker="o", linestyle="-") + plt.loglog(n, 0.1 * 1 / np.sqrt(n), linestyle="--", color="black") + if includeAngular: + plt.legend( + ["Measured error (linear)", "Measured error (dipole)", "O(1/sqrt(N))"] + ) + else: + plt.legend(["Measured error", "O(1/sqrt(N))"]) + plt.xlabel("RFD count") + plt.ylabel("Error (||running_avg - sol||)") + plt.title("Error of running average of RFDs compared to deterministic divM") + plt.grid(True) + plt.tight_layout() + plt.savefig("rfd_error.png", dpi=150) + plt.close() + + +def deterministic_div_m(solver, positions, includeAngular, delta): + N = np.size(positions) // 3 + size = 6 * N if includeAngular else 3 * N + + positions = positions.flatten() + div_M = np.zeros(size, dtype=positions.dtype) + for j in range(3 * N): + pos_plus = positions.copy() + pos_minus = positions.copy() + pos_plus[j] += 0.5 * delta + pos_minus[j] -= 0.5 * delta + + e_j = np.zeros(size, dtype=positions.dtype) + e_j[j] = 1.0 + f = e_j[0 : 3 * N] if includeAngular else e_j + t = e_j[3 * N :] if includeAngular else None + + solver.setPositions(pos_plus) + Mf_plus, Mt_plus = solver.Mdot(forces=f, torques=t) + + solver.setPositions(pos_minus) + Mf_minus, Mt_minus = solver.Mdot(forces=f, torques=t) + + div_M[0 : 3 * N] += Mf_plus - Mf_minus + if includeAngular: + div_M[3 * N :] += Mt_plus - Mt_minus + + div_M /= delta + return div_M diff --git a/tests/utils.py b/tests/utils.py index 5e95380..5b414f0 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -4,14 +4,15 @@ sane_parameters = { "PSE": {"psi": 1.0, "Lx": 32, "Ly": 32, "Lz": 32, "shearStrain": 0.0}, - "NBody": {"algorithm": "advise"}, - "NBody_wall": {"algorithm": "advise", "wallHeight": 0.0}, + "NBody": {"algorithm": "advise", "delta": 1e-3}, + "NBody_wall": {"algorithm": "advise", "wallHeight": 0.0, "delta": 1e-3}, "DPStokes": { "Lx": 16, "Ly": 16, "zmin": -6, "zmax": 6, "allowChangingBoxSize": False, + "delta": 1e-3, }, "SelfMobility": {"parameter": 5.0}, }