-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
commit ac2b405 Author: Ryker Fish <rfish@mines.edu> Date: Tue Jul 16 19:41:16 2024 -0600 regenerated data w new defaults after verifying against old repository commit 07b4794 Author: Ryker Fish <rfish@mines.edu> Date: Tue Jul 16 19:05:08 2024 -0600 added parameter to initialize that sets whether the box size can get altered or not commit 447e543 Author: Ryker Fish <rfish@mines.edu> Date: Tue Jul 16 18:22:27 2024 -0600 added mode that preserves L commit 74ec436 Author: Ryker Fish <rfish@mines.edu> Date: Tue Jul 16 17:23:24 2024 -0600 checkpoint commit- adding functionality to have parameter select that modifies the box size and another version that preserves it commit 1113390 Author: Ryker Fish <rfish@mines.edu> Date: Wed Jul 10 14:59:33 2024 -0600 the changes to DPStokes parameters seem to made the mobility matrix less symmetric commit 800f634 Author: Ryker Fish <rfish@mines.edu> Date: Wed Jul 10 14:56:08 2024 -0600 had to change the ref file to use the right positions. nbody test working. commit 830866e Author: Ryker Fish <rfish@mines.edu> Date: Mon Jul 8 18:46:51 2024 -0600 the two pair DPStokes tests are passing. currently having issues with the nbody test commit 9b3a4d3 Author: Ryker Fish <rfish@mines.edu> Date: Fri Jun 28 14:19:28 2024 -0600 self mobility tests done commit c4f2664 Author: Ryker Fish <rfish@mines.edu> Date: Wed Jun 12 15:31:55 2024 -0600 larger safety factor near open boundary for dpstokes commit 5da7450 Author: Ryker Fish <rfish@mines.edu> Date: Wed Jun 12 12:15:33 2024 -0600 modify how N gets rounded in all dimensions to be fft friendly commit 51db36a Author: Ryker Fish <rfish@mines.edu> Date: Wed Jun 12 11:58:39 2024 -0600 changed selection of grid parameters to not modify h. instead, modify the box size to make L/h an integer commit 5127394 Author: Ryker Fish <rfish@mines.edu> Date: Thu Jun 6 14:07:01 2024 -0600 changed test pass threshold commit 4fb82bf Author: Ryker Fish <rfish@mines.edu> Date: Thu Jun 6 13:42:52 2024 -0600 regenerated all ref data using libmobility implementations in single precision commit 4392903 Author: Ryker Fish <rfish@mines.edu> Date: Thu Jun 6 13:25:04 2024 -0600 regenerated refs using gpu implementations from dpstokes library. tests are now passing if the same kernel parameters are used as the ref implementation commit 1ce606a Author: Ryker Fish <rfish@mines.edu> Date: Thu May 30 16:06:32 2024 -0600 fixed bug in nbody wall correction that came from normalizing by a factor of 1/8pi in the original implementation instead of 1/6pi in the Swan paper. commit edf78f5 Author: Ryker Fish <rfish@mines.edu> Date: Thu May 30 16:05:21 2024 -0600 switched all tests to w=4 kernels, added nbody to pair mobility tests commit c575ef2 Author: Ryker Fish <rfish@mines.edu> Date: Wed May 29 15:11:48 2024 -0600 change dpstokes to use w=4 kernel by default for forces commit 070728b Author: Ryker Fish <rfish@mines.edu> Date: Tue May 28 16:57:09 2024 -0600 started pair mobility test commit ca4e1c7 Author: Ryker Fish <rfish@mines.edu> Date: Sat May 18 12:19:26 2024 -0600 combined self mobility tests into one function commit e90e9e2 Author: Ryker Fish <rfish@mines.edu> Date: Fri May 17 20:45:24 2024 -0600 bug fix in NBody with wall correction- test now passing. the fix uses the correct constants from the Swan correction paper- the previous implementation was using constants based on a self mobility constant of 1/(8*pi*eta*a) instead of 1/(6*pi*eta*a) commit 98e2ae0 Author: Ryker Fish <rfish@mines.edu> Date: Fri May 17 15:37:25 2024 -0600 changed dpstokes test to compare to dpstokes reference and added a not-working nbody test commit 7aaa0b1 Merge: 79fe6a7 d6144a5 Author: Ryker <rfish@mines.edu> Date: Mon May 13 14:25:35 2024 -0600 Merge pull request #18 from stochasticHydroTools/minor_test_edits Minor test edits commit d6144a5 Author: Ryker Fish <rfish@mines.edu> Date: Thu May 9 15:40:00 2024 -0600 removed test cases for small numbers of particles for the fluctuation dissipation test. commit 0e368d4 Author: Ryker <rfish@mines.edu> Date: Thu May 9 14:37:44 2024 -0600 Edit comment on nbody test Comment had the answer for the PSE test. Removed the periodic correction part commit 79fe6a7 Author: Ryker Fish <rfish@mines.edu> Date: Thu May 9 12:47:42 2024 -0600 now testing if mobility is 0 on the bottom wall commit 3c5002e Author: Ryker Fish <rfish@mines.edu> Date: Wed May 8 19:02:35 2024 -0600 possible correction to DPStokes code to add a buffer the top instead of the bottom when there is a wall commit c66303c Author: Ryker Fish <rfish@mines.edu> Date: Wed May 8 18:27:40 2024 -0600 selfmobility test working in comparison to refernce results for one and two wall cases commit 70aa329 Author: Ryker Fish <rfish@mines.edu> Date: Wed May 8 17:12:19 2024 -0600 mobility doesn't go to zero at the wall unless this section is commented out. commit fe2c964 Author: Ryker Fish <rfish@mines.edu> Date: Wed May 8 16:37:04 2024 -0600 first pass at computing translational self mobility components for bottom wall scenario
- Loading branch information
Showing
16 changed files
with
277 additions
and
26 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,30 @@ | ||
| #include <vector> | ||
|
|
||
| namespace dpstokes_polys{ | ||
|
|
||
| double polyEval(std::vector<double> polyCoeffs, double x){ | ||
|
|
||
| int order = polyCoeffs.size() - 1; | ||
| double accumulator = polyCoeffs[order]; | ||
| double current_x = x; | ||
| for(int i = 1; i <= order; i++){ | ||
| accumulator += polyCoeffs[order-i]*current_x; | ||
| current_x *= x; | ||
| } | ||
|
|
||
| return accumulator; | ||
| } | ||
|
|
||
| std::vector<double> cbetam_inv = { | ||
| 4131643418.193291, | ||
| -10471683395.26777, | ||
| 11833009228.6429, | ||
| -7851132955.882548, | ||
| 3388121732.651829, | ||
| -994285251.2185925, | ||
| 201183449.7086889, | ||
| -27776767.88241613, | ||
| 2515647.646492857, | ||
| -136305.2970161326, | ||
| 3445.959503226691}; | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Binary file not shown.
Binary file not shown.
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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. |
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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"), | ||
| (DPStokes, ("periodic", "periodic", "single_wall"), 1e-6, 0, "self_mobility_bw_w4.mat"), | ||
| (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, f"Self mobility does not match reference. Average error: {avgErr}" | ||
|
|
||
| @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}. Average error: {avgErr}" |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters