diff --git a/.clang-format b/.clang-format new file mode 100644 index 00000000..25563f42 --- /dev/null +++ b/.clang-format @@ -0,0 +1,274 @@ +--- +Language: Cpp +AccessModifierOffset: -2 +AlignAfterOpenBracket: Align +AlignArrayOfStructures: None +AlignConsecutiveAssignments: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: true +AlignConsecutiveBitFields: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: false +AlignConsecutiveDeclarations: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: false +AlignConsecutiveMacros: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: false +AlignConsecutiveShortCaseStatements: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCaseArrows: false + AlignCaseColons: false +AlignConsecutiveTableGenBreakingDAGArgColons: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: false +AlignConsecutiveTableGenCondOperatorColons: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: false +AlignConsecutiveTableGenDefinitionColons: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: false +AlignEscapedNewlines: Right +AlignOperands: Align +AlignTrailingComments: + Kind: Always + OverEmptyLines: 0 +AllowAllArgumentsOnNextLine: true +AllowAllParametersOfDeclarationOnNextLine: true +AllowBreakBeforeNoexceptSpecifier: Never +AllowShortBlocksOnASingleLine: Never +AllowShortCaseExpressionOnASingleLine: true +AllowShortCaseLabelsOnASingleLine: false +AllowShortCompoundRequirementOnASingleLine: true +AllowShortEnumsOnASingleLine: true +AllowShortFunctionsOnASingleLine: All +AllowShortIfStatementsOnASingleLine: Never +AllowShortLambdasOnASingleLine: All +AllowShortLoopsOnASingleLine: false +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakBeforeMultilineStrings: false +AttributeMacros: + - __capability +BinPackArguments: true +BinPackParameters: true +BitFieldColonSpacing: Both +BraceWrapping: + AfterCaseLabel: false + AfterClass: false + AfterControlStatement: Never + AfterEnum: false + AfterExternBlock: false + AfterFunction: false + AfterNamespace: false + AfterObjCDeclaration: false + AfterStruct: false + AfterUnion: false + BeforeCatch: false + BeforeElse: false + BeforeLambdaBody: false + BeforeWhile: false + IndentBraces: false + SplitEmptyFunction: true + SplitEmptyRecord: true + SplitEmptyNamespace: true +BreakAdjacentStringLiterals: true +BreakAfterAttributes: Leave +BreakAfterJavaFieldAnnotations: false +BreakAfterReturnType: None +BreakArrays: true +BreakBeforeBinaryOperators: None +BreakBeforeConceptDeclarations: Always +BreakBeforeBraces: Attach +BreakBeforeInlineASMColon: OnlyMultiline +BreakBeforeTernaryOperators: true +BreakConstructorInitializers: BeforeColon +BreakFunctionDefinitionParameters: false +BreakInheritanceList: BeforeColon +BreakStringLiterals: true +BreakTemplateDeclarations: MultiLine +ColumnLimit: 80 +CommentPragmas: '^ IWYU pragma:' +CompactNamespaces: false +ConstructorInitializerIndentWidth: 4 +ContinuationIndentWidth: 4 +Cpp11BracedListStyle: true +DerivePointerAlignment: false +DisableFormat: false +EmptyLineAfterAccessModifier: Never +EmptyLineBeforeAccessModifier: LogicalBlock +ExperimentalAutoDetectBinPacking: false +FixNamespaceComments: true +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +IfMacros: + - KJ_IF_MAYBE +IncludeBlocks: Preserve +IncludeCategories: + - Regex: '^"(llvm|llvm-c|clang|clang-c)/' + Priority: 2 + SortPriority: 0 + CaseSensitive: false + - Regex: '^(<|"(gtest|gmock|isl|json)/)' + Priority: 3 + SortPriority: 0 + CaseSensitive: false + - Regex: '.*' + Priority: 1 + SortPriority: 0 + CaseSensitive: false +IncludeIsMainRegex: '(Test)?$' +IncludeIsMainSourceRegex: '' +IndentAccessModifiers: false +IndentCaseBlocks: false +IndentCaseLabels: false +IndentExternBlock: AfterExternBlock +IndentGotoLabels: true +IndentPPDirectives: None +IndentRequiresClause: true +IndentWidth: 2 +IndentWrappedFunctionNames: false +InsertBraces: false +InsertNewlineAtEOF: false +InsertTrailingCommas: None +IntegerLiteralSeparator: + Binary: 0 + BinaryMinDigits: 0 + Decimal: 0 + DecimalMinDigits: 0 + Hex: 0 + HexMinDigits: 0 +JavaScriptQuotes: Leave +JavaScriptWrapImports: true +KeepEmptyLines: + AtEndOfFile: false + AtStartOfBlock: true + AtStartOfFile: true +LambdaBodyIndentation: Signature +LineEnding: DeriveLF +MacroBlockBegin: '' +MacroBlockEnd: '' +MainIncludeChar: Quote +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: None +ObjCBinPackProtocolList: Auto +ObjCBlockIndentWidth: 2 +ObjCBreakBeforeNestedBlockParam: true +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: true +PackConstructorInitializers: BinPack +PenaltyBreakAssignment: 2 +PenaltyBreakBeforeFirstCallParameter: 19 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakOpenParenthesis: 0 +PenaltyBreakScopeResolution: 500 +PenaltyBreakString: 1000 +PenaltyBreakTemplateDeclaration: 10 +PenaltyExcessCharacter: 1000000 +PenaltyIndentedWhitespace: 0 +PenaltyReturnTypeOnItsOwnLine: 60 +PointerAlignment: Right +PPIndentWidth: -1 +QualifierAlignment: Leave +ReferenceAlignment: Pointer +ReflowComments: true +RemoveBracesLLVM: false +RemoveParentheses: Leave +RemoveSemicolon: false +RequiresClausePosition: OwnLine +RequiresExpressionIndentation: OuterScope +SeparateDefinitionBlocks: Leave +ShortNamespaceLines: 1 +SkipMacroDefinitionBody: false +SortIncludes: CaseSensitive +SortJavaStaticImport: Before +SortUsingDeclarations: LexicographicNumeric +SpaceAfterCStyleCast: false +SpaceAfterLogicalNot: false +SpaceAfterTemplateKeyword: true +SpaceAroundPointerQualifiers: Default +SpaceBeforeAssignmentOperators: true +SpaceBeforeCaseColon: false +SpaceBeforeCpp11BracedList: false +SpaceBeforeCtorInitializerColon: true +SpaceBeforeInheritanceColon: true +SpaceBeforeJsonColon: false +SpaceBeforeParens: ControlStatements +SpaceBeforeParensOptions: + AfterControlStatements: true + AfterForeachMacros: true + AfterFunctionDefinitionName: false + AfterFunctionDeclarationName: false + AfterIfMacros: true + AfterOverloadedOperator: false + AfterPlacementOperator: true + AfterRequiresInClause: false + AfterRequiresInExpression: false + BeforeNonEmptyParentheses: false +SpaceBeforeRangeBasedForLoopColon: true +SpaceBeforeSquareBrackets: false +SpaceInEmptyBlock: false +SpacesBeforeTrailingComments: 1 +SpacesInAngles: Never +SpacesInContainerLiterals: true +SpacesInLineCommentPrefix: + Minimum: 1 + Maximum: -1 +SpacesInParens: Never +SpacesInParensOptions: + ExceptDoubleParentheses: false + InCStyleCasts: false + InConditionalStatements: false + InEmptyParentheses: false + Other: false +SpacesInSquareBrackets: false +Standard: Latest +StatementAttributeLikeMacros: + - Q_EMIT +StatementMacros: + - Q_UNUSED + - QT_REQUIRE_VERSION +TableGenBreakInsideDAGArg: DontBreak +TabWidth: 8 +UseTab: Never +VerilogBreakBetweenInstancePorts: true +WhitespaceSensitiveMacros: + - BOOST_PP_STRINGIZE + - CF_SWIFT_NAME + - NS_SWIFT_NAME + - PP_STRINGIZE + - STRINGIZE +... + diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs new file mode 100644 index 00000000..a6204fef --- /dev/null +++ b/.git-blame-ignore-revs @@ -0,0 +1,2 @@ +e4a0546006166dab68cdceebee9b62da8865abbc +aa1c33ec08d207d68753fe7e2e16ee142e69852b \ No newline at end of file diff --git a/.github/workflows/format-python.yml b/.github/workflows/format-python.yml new file mode 100644 index 00000000..f577e7cf --- /dev/null +++ b/.github/workflows/format-python.yml @@ -0,0 +1,94 @@ +name: Python Formatting Check + +on: + pull_request: + types: [opened, synchronize, reopened] + +jobs: + python-black-check: + runs-on: ubuntu-latest + + steps: + - name: Checkout code + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: '3.11' + + - name: Install Black + run: pip install black + + - name: Run Black and extract offending files + id: black_check + continue-on-error: true + run: | + set +e + OUTPUT=$(black --check . 2>&1) + EXIT_CODE=$? + + echo "$OUTPUT" > black_output.txt + grep '^would reformat' black_output.txt | awk '{print $3}' | sed "s|$(pwd)/||" | sort -u > black_files.txt + + if [[ -s black_files.txt ]]; then + echo "py_failed=true" >> $GITHUB_ENV + echo "black_files<> $GITHUB_OUTPUT + cat black_files.txt >> $GITHUB_OUTPUT + echo "EOF" >> $GITHUB_OUTPUT + else + echo "py_failed=false" >> $GITHUB_ENV + echo "black_files=" >> $GITHUB_OUTPUT + fi + + - name: Find existing Black formatting comment + id: find_comment + uses: peter-evans/find-comment@v3 + with: + issue-number: ${{ github.event.pull_request.number }} + comment-author: 'github-actions[bot]' + body-includes: '' + + - name: Create or update PR comment - failed + if: env.py_failed == 'true' + uses: peter-evans/create-or-update-comment@v4 + with: + token: ${{ secrets.GITHUB_TOKEN }} + issue-number: ${{ github.event.pull_request.number }} + comment-id: ${{ steps.find_comment.outputs.comment-id }} + edit-mode: replace + body: | + + ⚠️ **Black formatting check failed** + + The following Python files need formatting: + + ``` + ${{ steps.black_check.outputs.black_files }} + ``` + + Please run: + ```bash + black . + ``` + And commit the changes before merging. + + - name: Create or update PR comment - success + if: env.py_failed == 'false' + uses: peter-evans/create-or-update-comment@v4 + with: + token: ${{ secrets.GITHUB_TOKEN }} + issue-number: ${{ github.event.pull_request.number }} + comment-id: ${{ steps.find_comment.outputs.comment-id }} + edit-mode: replace + body: | + + ✅ **Linter reported no issues** + + All Python files are correctly formatted with **Black**. + + - name: Fail if Black failed + if: env.py_failed == 'true' + run: | + echo "❌ Python formatting issues found." + exit 1 diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml new file mode 100644 index 00000000..5f1f785e --- /dev/null +++ b/.github/workflows/format.yml @@ -0,0 +1,90 @@ +name: Clang-Format Lint Check + +on: + pull_request: + types: [opened, synchronize, reopened] + +jobs: + clang-format-check: + runs-on: ubuntu-latest + + steps: + - name: Checkout code + uses: actions/checkout@v4 + + - name: Install clang-format + run: sudo apt-get update && sudo apt-get install -y clang-format-19 + + - name: Find unformatted files + id: format_check + run: | + FILE_LIST="unformatted_files.txt" + find solvers include ! -path 'src/third_party/*' \ + -regex '.*\.\(cpp\|hpp\|h\|c\|cc\|cu\|cuh\)' \ + -exec bash -c 'clang-format-19 -style=file -output-replacements-xml "$0" | grep -q " "$FILE_LIST" + + cat "$FILE_LIST" + + if [ -s "$FILE_LIST" ]; then + echo "failed=true" >> $GITHUB_ENV + { + echo 'files<> "$GITHUB_OUTPUT" + else + echo "failed=false" >> $GITHUB_ENV + echo 'files=' >> "$GITHUB_OUTPUT" + fi + + - name: Find existing clang-format comment + id: find_comment + uses: peter-evans/find-comment@v3 + with: + issue-number: ${{ github.event.pull_request.number }} + comment-author: 'github-actions[bot]' + body-includes: '' + + - name: Create or update PR comment - failed + if: env.failed == 'true' + uses: peter-evans/create-or-update-comment@v4 + with: + token: ${{ secrets.GITHUB_TOKEN }} + issue-number: ${{ github.event.pull_request.number }} + comment-id: ${{ steps.find_comment.outputs.comment-id }} + edit-mode: replace + body: | + + ⚠️ **Clang-format check failed** + + The following files are not correctly formatted: + + ``` + ${{ steps.format_check.outputs.files }} + ``` + + Please run: + ```bash + clang-format -i + ``` + And commit the changes before merging. + + - name: Create or update PR comment - success + if: env.failed == 'false' + uses: peter-evans/create-or-update-comment@v4 + with: + token: ${{ secrets.GITHUB_TOKEN }} + issue-number: ${{ github.event.pull_request.number }} + comment-id: ${{ steps.find_comment.outputs.comment-id }} + edit-mode: replace + body: | + + ✅ **Linter reported no issues** + + All C/C++ files are correctly formatted with **clang-format**. + + - name: Fail job if formatting issues found + if: env.failed == 'true' + run: | + echo "❌ Formatting issues detected. Failing the job." + exit 1 diff --git a/docs/source/new-solver.rst b/docs/source/new-solver.rst index eddff5a5..b48cc42f 100644 --- a/docs/source/new-solver.rst +++ b/docs/source/new-solver.rst @@ -42,15 +42,11 @@ A new solver must add a directory under the "solvers" folder, which must be the This must contain rules to create the shared library for the particular solver and its Python wrapper. The solver library should be called "lib[Solver].so", while the Python library should be called "[Solver].[Python_SOABI].so" with the correct extension suffix. -4. **example.py** +4. **README.md** - An example script using the Python bindings of the module. + The documentation for the specific module (anything relevant about the particular solver that does not fit the class docstring). -5. **README.md** - - The documentation for the specific module (stating the geometry, required arguments, etc.). - -6. **tests** +5. **tests** Add any relevant tests to the test folder using pytest. diff --git a/include/MobilityInterface/MobilityInterface.h b/include/MobilityInterface/MobilityInterface.h index 5c0107ee..bd6f6d13 100644 --- a/include/MobilityInterface/MobilityInterface.h +++ b/include/MobilityInterface/MobilityInterface.h @@ -144,8 +144,8 @@ class Mobility { "number of particles"); } // if (this->needsTorque && linear.size() != angular.size()) - // throw std::runtime_error("[libMobility] This solver requires angular " - // "velocities when configured with torques"); + // throw std::runtime_error("[libMobility] This solver requires angular " + // "velocities when configured with torques"); const auto numberElements = numberParticles + (this->includeAngular ? numberParticles : 0); if (not lanczos) { diff --git a/include/MobilityInterface/lanczos.h b/include/MobilityInterface/lanczos.h index 712984d4..aa3e99a4 100644 --- a/include/MobilityInterface/lanczos.h +++ b/include/MobilityInterface/lanczos.h @@ -3,36 +3,39 @@ */ #ifndef LIBMOBILITY_LANCZOS_ADAPTOR_H #define LIBMOBILITY_LANCZOS_ADAPTOR_H -#include"LanczosAlgorithm.h" -#include -#include -#include +#include "LanczosAlgorithm.h" #include +#include +#include +#include // This class uses the LanczosAlgorithm library to compute fluctuations. -class LanczosStochasticVelocities{ +class LanczosStochasticVelocities { using real = lanczos::real; lanczos::Solver lanczos; std::vector lanczosNoise; real lanczosTolerance; std::mt19937 engine; -public: - LanczosStochasticVelocities(real tol, std::uint64_t seed){ +public: + LanczosStochasticVelocities(real tol, std::uint64_t seed) { this->lanczosTolerance = tol; engine = std::mt19937{seed}; } - //Given a functor that applies the mobility operator, returns prefactor*(B dW). Where B is an operator that applies the square root of the provided mobility. - template - void sqrtMdotW(MobilityDot dot, real* result, int numberParticles, real prefactor = 1){ - std::normal_distribution dist {0, 1}; - auto gen = [&](){return dist(engine);}; - lanczosNoise.resize(3*numberParticles); + // Given a functor that applies the mobility operator, returns prefactor*(B + // dW). Where B is an operator that applies the square root of the provided + // mobility. + template + void sqrtMdotW(MobilityDot dot, real *result, int numberParticles, + real prefactor = 1) { + std::normal_distribution dist{0, 1}; + auto gen = [&]() { return dist(engine); }; + lanczosNoise.resize(3 * numberParticles); std::generate(lanczosNoise.begin(), lanczosNoise.end(), gen); - lanczos.run(dot, result, lanczosNoise.data(), lanczosTolerance, 3*numberParticles); + lanczos.run(dot, result, lanczosNoise.data(), lanczosTolerance, + 3 * numberParticles); } - }; #endif diff --git a/include/MobilityInterface/pythonify.h b/include/MobilityInterface/pythonify.h index e77f584e..37878e00 100644 --- a/include/MobilityInterface/pythonify.h +++ b/include/MobilityInterface/pythonify.h @@ -104,39 +104,42 @@ auto check_and_get_shape(pyarray_c &arr) { template auto setup_arrays(Solver &myself, pyarray_c &forces, pyarray_c &torques) { - size_t N = myself.getNumberParticles(); - if (forces.size() > 0 && forces.size() != 3 * N) { - throw std::runtime_error("The forces array must have size 3*N."); - } - if (torques.size() > 0 && torques.size() != 3 * N) { - throw std::runtime_error("The torques array must have size 3*N."); - } + size_t N = myself.getNumberParticles(); + if (forces.size() > 0 && forces.size() != 3 * N) { + throw std::runtime_error("The forces array must have size 3*N."); + } + if (torques.size() > 0 && torques.size() != 3 * N) { + throw std::runtime_error("The torques array must have size 3*N."); + } - auto f = cast_to_const_real(forces); - auto t = cast_to_const_real(torques); - auto framework = lp::get_framework(forces); - last_framework = framework; - int device = f.empty() ? torques.device_type() : forces.device_type(); - if (device == nb::device::none::value) - device = nb::device::cpu::value; - last_device = device; - - auto mf = nb::ndarray(); - auto mt = nb::ndarray(); - - auto f_shape = f.empty() ? check_and_get_shape(torques) : check_and_get_shape(forces); - mf = lp::create_with_framework(f_shape, device, framework); - - if (!myself.getIncludeAngular()) { - if(!t.empty()){ - throw std::runtime_error("The solver was configured without including angular velocities. " - "Set includeAngular to true when initializing if you want to use torques"); - } - } else{ - // only set mt if includeAngular is true - auto t_shape = t.empty() ? f_shape : check_and_get_shape(torques); - mt = lp::create_with_framework(t_shape, device, framework); + auto f = cast_to_const_real(forces); + auto t = cast_to_const_real(torques); + auto framework = lp::get_framework(forces); + last_framework = framework; + int device = f.empty() ? torques.device_type() : forces.device_type(); + if (device == nb::device::none::value) + device = nb::device::cpu::value; + last_device = device; + + auto mf = nb::ndarray(); + auto mt = nb::ndarray(); + + auto f_shape = + f.empty() ? check_and_get_shape(torques) : check_and_get_shape(forces); + mf = lp::create_with_framework(f_shape, device, framework); + + if (!myself.getIncludeAngular()) { + if (!t.empty()) { + throw std::runtime_error( + "The solver was configured without including angular velocities. " + "Set includeAngular to true when initializing if you want to use " + "torques"); } + } else { + // only set mt if includeAngular is true + auto t_shape = t.empty() ? f_shape : check_and_get_shape(torques); + mt = lp::create_with_framework(t_shape, device, framework); + } return std::make_tuple(f, t, mf, mt); } @@ -218,12 +221,12 @@ array_like template auto call_mdot(Solver &myself, pyarray_c &forces, pyarray_c &torques) { - auto [f, t, mf, mt] = setup_arrays(myself, forces, torques); - int N = myself.getNumberParticles(); - auto mf_ptr = cast_to_real(mf); - auto mt_ptr = cast_to_real(mt); - myself.Mdot(f, t, mf_ptr, mt_ptr); - return std::make_pair(mf, mt); + auto [f, t, mf, mt] = setup_arrays(myself, forces, torques); + int N = myself.getNumberParticles(); + auto mf_ptr = cast_to_real(mf); + auto mt_ptr = cast_to_real(mt); + myself.Mdot(f, t, mf_ptr, mt_ptr); + return std::make_pair(mf, mt); } const char *mdot_docstring = R"pbdoc( @@ -250,8 +253,8 @@ array_like )pbdoc"; template -void call_initialize(Solver &myself, real T, real eta, real a, bool includeAngular, - real tol) { +void call_initialize(Solver &myself, real T, real eta, real a, + bool includeAngular, real tol) { libmobility::Parameters par; par.temperature = T; par.viscosity = eta; @@ -267,9 +270,8 @@ template void call_setPositions(Solver &myself, pyarray_c &pos) { last_device = pos.device_type(); myself.setPositions(cast_to_const_real(pos)); last_shape.resize(pos.ndim()); - for (size_t i = 0; i < pos.ndim(); ++i) - { - last_shape[i] = pos.shape(i); + for (size_t i = 0; i < pos.ndim(); ++i) { + last_shape[i] = pos.shape(i); } } diff --git a/include/memory/allocator.h b/include/memory/allocator.h index 24d5f3fb..92a2d65a 100644 --- a/include/memory/allocator.h +++ b/include/memory/allocator.h @@ -11,7 +11,8 @@ using namespace libmobility::allocator; using float_vector = std::vector>; - using float_device_vector = thrust::device_vector>; + using float_device_vector = thrust::device_vector>; float_vector v(1000); // Allocates 1000 floats in host memory float_device_vector dv(1000); // Allocates 1000 floats in device memory diff --git a/include/memory/python_tensor.h b/include/memory/python_tensor.h index 26dd7fcb..81cd91fa 100644 --- a/include/memory/python_tensor.h +++ b/include/memory/python_tensor.h @@ -5,12 +5,12 @@ */ #pragma once +#include "allocator.h" #include #include #include #include #include -#include "allocator.h" namespace libmobility { namespace python { namespace nb = nanobind; diff --git a/setup.py b/setup.py index 33eeb6aa..024ffd58 100644 --- a/setup.py +++ b/setup.py @@ -40,9 +40,13 @@ def build_extension(self, ext): sys.prefix ) # Use sys.prefix to install in the environment's prefix - double_precision = os.environ.get("DOUBLEPRECISION", "0") # Default to "0" if not set + double_precision = os.environ.get( + "DOUBLEPRECISION", "0" + ) # Default to "0" if not set if double_precision not in {"0", "1"}: - raise ValueError("Environment variable DOUBLEPRECISION must be '0' or '1'. It is 0 by default- set to 1 to compile in double precision.") + raise ValueError( + "Environment variable DOUBLEPRECISION must be '0' or '1'. It is 0 by default- set to 1 to compile in double precision." + ) cmake_args = [ "-DCMAKE_INSTALL_PREFIX=" + install_dir, # Point to the install directory diff --git a/solvers/DPStokes/extra/poly_fits.h b/solvers/DPStokes/extra/poly_fits.h index 878fdb1f..0ea11a55 100644 --- a/solvers/DPStokes/extra/poly_fits.h +++ b/solvers/DPStokes/extra/poly_fits.h @@ -1,49 +1,45 @@ /* Ryker Fish 2024. Helpers to compute values of c^{-1} from [1]. -The ES kernel has a function c(beta) that provides a relationship between the hydrodynamic radius and the kernel -parameters alpha and beta. There is no closed form for c, but we can approximate c and c^{-1} with a polynomial fit. +The ES kernel has a function c(beta) that provides a relationship between the +hydrodynamic radius and the kernel parameters alpha and beta. There is no closed +form for c, but we can approximate c and c^{-1} with a polynomial fit. -We have rH = h*w*c(beta), where rH is the hydrodynamic radius, h is the grid spacing, and w is the width of the kernel. -Thus, we can write c^{-1}(rH/(h*w)) = beta. +We have rH = h*w*c(beta), where rH is the hydrodynamic radius, h is the grid +spacing, and w is the width of the kernel. Thus, we can write c^{-1}(rH/(h*w)) = +beta. References: -[1] Computing hydrodynamic interactions in confined doubly periodic geometries in linear time. -A. Hashemi et al. J. Chem. Phys. 158, 154101 (2023) https://doi.org/10.1063/5.0141371 +[1] Computing hydrodynamic interactions in confined doubly periodic geometries +in linear time. A. Hashemi et al. J. Chem. Phys. 158, 154101 (2023) +https://doi.org/10.1063/5.0141371 */ #include - -namespace dpstokes_polys{ - - /* Evaluates a polynomial at x with coefficients in descending order, so the highest order coefficients - are at the start of polyCoeffs. - e.g. for a polynomial of order n, this computes - polyCoeffs[n+1] + polyCoeffs[n]*x + polyCoeffs[n-2]*x^2 + ... + polyCoeffs[0]*x^n - */ - double polyEval(std::vector 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; + +namespace dpstokes_polys { + +/* Evaluates a polynomial at x with coefficients in descending order, so the +highest order coefficients are at the start of polyCoeffs. e.g. for a polynomial +of order n, this computes polyCoeffs[n+1] + polyCoeffs[n]*x + +polyCoeffs[n-2]*x^2 + ... + polyCoeffs[0]*x^n +*/ +double polyEval(std::vector 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; } - // Coefficients for the polynomial fit of c^{-1} from [1] - std::vector cbetam_inv = { - 4131643418.193291, - -10471683395.26777, - 11833009228.6429, - -7851132955.882548, - 3388121732.651829, - -994285251.2185925, - 201183449.7086889, - -27776767.88241613, - 2515647.646492857, - -136305.2970161326, - 3445.959503226691}; -} \ No newline at end of file + return accumulator; +} + +// Coefficients for the polynomial fit of c^{-1} from [1] +std::vector cbetam_inv = { + 4131643418.193291, -10471683395.26777, 11833009228.6429, + -7851132955.882548, 3388121732.651829, -994285251.2185925, + 201183449.7086889, -27776767.88241613, 2515647.646492857, + -136305.2970161326, 3445.959503226691}; +} // 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 b6a031a6..cc502fa7 100644 --- a/solvers/DPStokes/extra/uammd_interface.h +++ b/solvers/DPStokes/extra/uammd_interface.h @@ -1,67 +1,72 @@ /*Raul P. Pelaez 2021. An interface code between uammd_wrapper.cu and uammd_python.cpp. */ -#include -#include -namespace uammd_dpstokes{ - //This is in order not to use any UAMMD related includes here. - //Instead of using uammd::real I have to re define real here. +#include +#include +namespace uammd_dpstokes { +// This is in order not to use any UAMMD related includes here. +// Instead of using uammd::real I have to re define real here. #ifndef DOUBLE_PRECISION - using real = float; +using real = float; #else - using real = double; +using real = double; #endif - //This function returns either 'single' or 'double' according to the UAMMD's compiled precision. - std::string getPrecision(); +// This function returns either 'single' or 'double' according to the UAMMD's +// compiled precision. +std::string getPrecision(); - struct PyParameters{ - //The number of cells in each direction - //If -1, they will be autocomputed from the tolerance if possible (DP cannot do it, FCM can) - int nx = -1; - int ny = -1; - int nz = -1; - real dt = 0; - real viscosity; - real Lx; - real Ly; - real zmin, zmax; - //Tolerance will be ignored in DP mode, TP will use only tolerance and nxy/nz - real tolerance = 1e-5; - real w, w_d; - real hydrodynamicRadius = -1; - real beta = -1; - real beta_d = -1; - real alpha = -1; - real alpha_d = -1; - //Can be either none, bottom, slit or periodic - std::string mode; - bool allowChangingBoxSize = false; - }; +struct PyParameters { + // The number of cells in each direction + // If -1, they will be autocomputed from the tolerance if possible (DP cannot + // do it, FCM can) + int nx = -1; + int ny = -1; + int nz = -1; + real dt = 0; + real viscosity; + real Lx; + real Ly; + real zmin, zmax; + // Tolerance will be ignored in DP mode, TP will use only tolerance and nxy/nz + real tolerance = 1e-5; + real w, w_d; + real hydrodynamicRadius = -1; + real beta = -1; + real beta_d = -1; + real alpha = -1; + real alpha_d = -1; + // Can be either none, bottom, slit or periodic + std::string mode; + bool allowChangingBoxSize = false; +}; - class DPStokesUAMMD; - class DPStokesGlue{ - std::shared_ptr dpstokes; - public: - int numberParticles; +class DPStokesUAMMD; +class DPStokesGlue { + std::shared_ptr dpstokes; - //Initialize the modules with a certain set of parameters - //Reinitializes if the module was already initialized - void initialize(PyParameters pypar); +public: + int numberParticles; - //Clears all memory allocated by the module. - //This leaves the module in an unusable state until initialize is called again. - void clear(); - //Set positions to compute mobility matrix - void setPositions(const real* h_pos, int numberParticles); + // Initialize the modules with a certain set of parameters + // Reinitializes if the module was already initialized + void initialize(PyParameters pypar); - const real* getStoredPositions(); - //Compute the dot product of the mobility matrix with the forces and/or torques acting on the previously provided positions - void Mdot(const real *h_forces, const real *h_torques, real *h_MF, - real *h_MT, int numberParticles, bool includeAngular); + // Clears all memory allocated by the module. + // This leaves the module in an unusable state until initialize is called + // again. + void clear(); + // Set positions to compute mobility matrix + void setPositions(const real *h_pos, int numberParticles); - private: - void throwIfInvalid(); - }; + const real *getStoredPositions(); + // Compute the dot product of the mobility matrix with the forces and/or + // torques acting on the previously provided positions + void Mdot(const real *h_forces, const real *h_torques, real *h_MF, real *h_MT, + int numberParticles, bool includeAngular); -} +private: + void throwIfInvalid(); +}; + +} // namespace uammd_dpstokes diff --git a/solvers/DPStokes/extra/uammd_wrapper.cu b/solvers/DPStokes/extra/uammd_wrapper.cu index f63d12b7..8838e5b6 100644 --- a/solvers/DPStokes/extra/uammd_wrapper.cu +++ b/solvers/DPStokes/extra/uammd_wrapper.cu @@ -1,260 +1,263 @@ /* Raul P. Pelaez 2021. Doubly Periodic Stokes UAMMD wrapper - Allows to call the DPStokes or TP FCM modules from via a simple contained class to compute the product between the mobility tensor and a list forces and torques acting on a group of positions. + Allows to call the DPStokes or TP FCM modules from via a simple contained + class to compute the product between the mobility tensor and a list forces + and torques acting on a group of positions. - Additionally, a glue class is provided to ease separate compilation between GPU code (this source) and another code. For instance, the python wrapper in uammd_python.cpp + Additionally, a glue class is provided to ease separate compilation between + GPU code (this source) and another code. For instance, the python wrapper in + uammd_python.cpp */ #include -//Doubly Periodic FCM implementation (currently without noise) +// Doubly Periodic FCM implementation (currently without noise) #include -//Triply Periodic FCM implementation +// Triply Periodic FCM implementation +#include "uammd_interface.h" #include -#include"uammd_interface.h" // Some convenient aliases -namespace uammd_dpstokes{ - using FCM_BM = uammd::BDHI::FCM_ns::Kernels::BarnettMagland; - using FCM = uammd::BDHI::FCM_impl; - using DPStokesSlab = uammd::DPStokesSlab_ns::DPStokes; - using uammd::DPStokesSlab_ns::WallMode; - using uammd::System; +namespace uammd_dpstokes { +using FCM_BM = uammd::BDHI::FCM_ns::Kernels::BarnettMagland; +using FCM = uammd::BDHI::FCM_impl; +using DPStokesSlab = uammd::DPStokesSlab_ns::DPStokes; +using uammd::System; +using uammd::DPStokesSlab_ns::WallMode; - //Helper functions and objects - struct Real3ToReal4{ - __host__ __device__ uammd::real4 operator()(uammd::real3 i){ - auto pr4 = uammd::make_real4(i); - return pr4; - } - }; - struct Real4ToReal3{ - __host__ __device__ uammd::real3 operator()(uammd::real4 i){ - auto pr3 = uammd::make_real3(i); - return pr3; - } - }; - - struct Real3ToReal4SubstractOriginZ{ - real origin; - Real3ToReal4SubstractOriginZ(real origin):origin(origin){} - __host__ __device__ uammd::real4 operator()(uammd::real3 i){ - auto pr4 = uammd::make_real4(i); - pr4.z -= origin; - return pr4; - } - }; - - auto createFCMParameters(PyParameters pypar){ - FCM::Parameters par; - par.temperature = 0; //FCM can compute fluctuations, but they are turned off here - par.viscosity = pypar.viscosity; - 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.kernelTorque = std::make_shared(pypar.w_d, pypar.alpha_d, pypar.beta_d, pypar.Lx/pypar.nx); - return par; +// Helper functions and objects +struct Real3ToReal4 { + __host__ __device__ uammd::real4 operator()(uammd::real3 i) { + auto pr4 = uammd::make_real4(i); + return pr4; } - - WallMode stringToWallMode(std::string str){ - if(str.compare("nowall") == 0){ - return WallMode::none; - } - else if(str.compare("slit") == 0){ - return WallMode::slit; - } - else if(str.compare("bottom") == 0){ - return WallMode::bottom; - } - else return WallMode::none; +}; +struct Real4ToReal3 { + __host__ __device__ uammd::real3 operator()(uammd::real4 i) { + auto pr3 = uammd::make_real3(i); + return pr3; } +}; - auto createDPStokesParameters(PyParameters pypar){ - DPStokesSlab::Parameters par; - 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; - par.alpha_d = pypar.alpha_d; - par.mode = stringToWallMode(pypar.mode); - return par; +struct Real3ToReal4SubstractOriginZ { + real origin; + Real3ToReal4SubstractOriginZ(real origin) : origin(origin) {} + __host__ __device__ uammd::real4 operator()(uammd::real3 i) { + auto pr4 = uammd::make_real4(i); + pr4.z -= origin; + return pr4; } +}; - //Wrapper to UAMMD's TP and DP hydrodynamic modules - struct DPStokesUAMMD { - 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){ - if(fcm){ - return fcm->computeHydrodynamicDisplacements((uammd::real4*)(d_pos), - (uammd::real4*)(d_force), - (uammd::real4*)(d_torques), - numberParticles, 0.0, 0.0, st); - } - else if(dpstokes){ - return dpstokes->Mdot(reinterpret_cast(d_pos), - reinterpret_cast(d_force), - reinterpret_cast(d_torques), - numberParticles, st); - } - } - public: - std::shared_ptr dpstokes; - std::shared_ptr fcm; - cudaStream_t st; - thrust::device_vector tmp3; - thrust::device_vector force4; - thrust::device_vector torque4; - thrust::device_vector stored_positions; - thrust::device_vector stored_positions_real3; - real zOrigin; +auto createFCMParameters(PyParameters pypar) { + FCM::Parameters par; + par.temperature = + 0; // FCM can compute fluctuations, but they are turned off here + par.viscosity = pypar.viscosity; + 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.kernelTorque = std::make_shared( + pypar.w_d, pypar.alpha_d, pypar.beta_d, pypar.Lx / pypar.nx); + return par; +} - DPStokesUAMMD(PyParameters pypar){ - if(pypar.mode.compare("periodic")==0){ - auto par = createFCMParameters(pypar); - this->fcm = std::make_shared(par); - zOrigin = 0; - } - else{ - auto par = createDPStokesParameters(pypar); - this->dpstokes = std::make_shared(par); - zOrigin = pypar.zmin + par.H*0.5; - } - CudaSafeCall(cudaStreamCreate(&st)); - } +WallMode stringToWallMode(std::string str) { + if (str.compare("nowall") == 0) { + return WallMode::none; + } else if (str.compare("slit") == 0) { + return WallMode::slit; + } else if (str.compare("bottom") == 0) { + return WallMode::bottom; + } else + return WallMode::none; +} + +auto createDPStokesParameters(PyParameters pypar) { + DPStokesSlab::Parameters par; + 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; + par.alpha_d = pypar.alpha_d; + par.mode = stringToWallMode(pypar.mode); + return par; +} - //Copy positions to UAMMD's ParticleData - void setPositions(const real* d_pos, int numberParticles){ - stored_positions_real3.resize(numberParticles); - // Copy onto stored_positions_real3 - thrust::copy(thrust::cuda::par.on(st), - reinterpret_cast(d_pos), - reinterpret_cast(d_pos) + numberParticles, - stored_positions_real3.begin()); +// Wrapper to UAMMD's TP and DP hydrodynamic modules +struct DPStokesUAMMD { +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) { + if (fcm) { + return fcm->computeHydrodynamicDisplacements( + (uammd::real4 *)(d_pos), (uammd::real4 *)(d_force), + (uammd::real4 *)(d_torques), numberParticles, 0.0, 0.0, st); + } else if (dpstokes) { + return dpstokes->Mdot(reinterpret_cast(d_pos), + reinterpret_cast(d_force), + reinterpret_cast(d_torques), + numberParticles, st); } + } + +public: + std::shared_ptr dpstokes; + std::shared_ptr fcm; + cudaStream_t st; + thrust::device_vector tmp3; + thrust::device_vector force4; + thrust::device_vector torque4; + thrust::device_vector stored_positions; + thrust::device_vector stored_positions_real3; + real zOrigin; - const real* getStoredPositions(){ - return reinterpret_cast(stored_positions_real3.data().get()); + DPStokesUAMMD(PyParameters pypar) { + if (pypar.mode.compare("periodic") == 0) { + auto par = createFCMParameters(pypar); + this->fcm = std::make_shared(par); + zOrigin = 0; + } else { + auto par = createDPStokesParameters(pypar); + this->dpstokes = std::make_shared(par); + zOrigin = pypar.zmin + par.H * 0.5; } + CudaSafeCall(cudaStreamCreate(&st)); + } - //Compute the hydrodynamic displacements due to a series of forces and/or torques acting on the particles - void Mdot(const real *h_forces, const real *h_torques, real *h_MF, - real *h_MT, int numberParticles, bool includeAngular) { - bool givenForces = h_forces; - bool givenTorques = h_torques; - force4.resize(numberParticles); - if(givenForces){ - thrust::transform( - thrust::cuda::par.on(st), - reinterpret_cast(h_forces), - reinterpret_cast(h_forces) + numberParticles, - force4.begin(), Real3ToReal4()); - } else { - thrust::fill(thrust::cuda::par.on(st), force4.begin(), force4.end(), uammd::make_real4(0, 0, 0, 0)); - } + // Copy positions to UAMMD's ParticleData + void setPositions(const real *d_pos, int numberParticles) { + stored_positions_real3.resize(numberParticles); + // Copy onto stored_positions_real3 + thrust::copy( + thrust::cuda::par.on(st), reinterpret_cast(d_pos), + reinterpret_cast(d_pos) + numberParticles, + stored_positions_real3.begin()); + } - if (givenTorques) { - torque4.resize(numberParticles); - thrust::transform(thrust::cuda::par.on(st), - reinterpret_cast(h_torques), - reinterpret_cast(h_torques) + - numberParticles, - torque4.begin(), Real3ToReal4()); - } else if (includeAngular) { - // If we are including angular forces but no torques were given, we fill - // with zeros to generate angular velocities - torque4.resize(numberParticles); - thrust::fill(thrust::cuda::par.on(st), torque4.begin(), torque4.end(), - uammd::make_real4(0, 0, 0, 0)); - } + const real *getStoredPositions() { + return reinterpret_cast(stored_positions_real3.data().get()); + } - stored_positions.resize(numberParticles); + // Compute the hydrodynamic displacements due to a series of forces and/or + // torques acting on the particles + void Mdot(const real *h_forces, const real *h_torques, real *h_MF, real *h_MT, + int numberParticles, bool includeAngular) { + bool givenForces = h_forces; + bool givenTorques = h_torques; + force4.resize(numberParticles); + if (givenForces) { thrust::transform(thrust::cuda::par.on(st), - stored_positions_real3.begin(), stored_positions_real3.end(), - stored_positions.begin(), Real3ToReal4SubstractOriginZ(zOrigin)); - auto mob = this->computeHydrodynamicDisplacements( - stored_positions.data().get(), force4.data().get(), - (includeAngular) ? torque4.data().get() : nullptr, numberParticles, - 0.0, 0.0, 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); - if (includeAngular) { - thrust::copy(thrust::cuda::par.on(st), mob.second.begin(), - mob.second.end(), (uammd::real3 *)h_MT); - } + reinterpret_cast(h_forces), + reinterpret_cast(h_forces) + + numberParticles, + force4.begin(), Real3ToReal4()); + } else { + thrust::fill(thrust::cuda::par.on(st), force4.begin(), force4.end(), + uammd::make_real4(0, 0, 0, 0)); } - ~DPStokesUAMMD(){ - cudaDeviceSynchronize(); - cudaStreamDestroy(st); + if (givenTorques) { + torque4.resize(numberParticles); + thrust::transform(thrust::cuda::par.on(st), + reinterpret_cast(h_torques), + reinterpret_cast(h_torques) + + numberParticles, + torque4.begin(), Real3ToReal4()); + } else if (includeAngular) { + // If we are including angular forces but no torques were given, we fill + // with zeros to generate angular velocities + torque4.resize(numberParticles); + thrust::fill(thrust::cuda::par.on(st), torque4.begin(), torque4.end(), + uammd::make_real4(0, 0, 0, 0)); } - }; - - //Initialize the modules with a certain set of parameters - //Reinitializes if the module was already initialized - void DPStokesGlue::initialize(PyParameters pypar){ - dpstokes = std::make_shared(pypar); + stored_positions.resize(numberParticles); + thrust::transform(thrust::cuda::par.on(st), stored_positions_real3.begin(), + stored_positions_real3.end(), stored_positions.begin(), + Real3ToReal4SubstractOriginZ(zOrigin)); + auto mob = this->computeHydrodynamicDisplacements( + stored_positions.data().get(), force4.data().get(), + (includeAngular) ? torque4.data().get() : nullptr, numberParticles, 0.0, + 0.0, 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); + if (includeAngular) { + thrust::copy(thrust::cuda::par.on(st), mob.second.begin(), + mob.second.end(), (uammd::real3 *)h_MT); + } } - //Clears all memory allocated by the module. - //This leaves the module in an unusable state until initialize is called again. - void DPStokesGlue::clear(){ - dpstokes.reset(); + ~DPStokesUAMMD() { + cudaDeviceSynchronize(); + cudaStreamDestroy(st); } +}; - //Set positions to compute mobility matrix - void DPStokesGlue::setPositions(const real* h_pos, int numberParticles){ - throwIfInvalid(); - dpstokes->setPositions(h_pos, numberParticles); - } +// Initialize the modules with a certain set of parameters +// Reinitializes if the module was already initialized +void DPStokesGlue::initialize(PyParameters pypar) { + dpstokes = std::make_shared(pypar); +} - //Get the positions of the particles in the system - const real* DPStokesGlue::getStoredPositions(){ - throwIfInvalid(); - return dpstokes->getStoredPositions(); - } +// Clears all memory allocated by the module. +// This leaves the module in an unusable state until initialize is called again. +void DPStokesGlue::clear() { dpstokes.reset(); } - //Compute the dot product of the mobility matrix with the forces and/or torques acting on the previously provided positions - void DPStokesGlue::Mdot(const real *h_forces, const real *h_torques, - real *h_MF, real *h_MT, int numberParticles, - bool includeAngular) { - throwIfInvalid(); - dpstokes->Mdot(h_forces, h_torques, h_MF, h_MT, numberParticles, - includeAngular); - } +// Set positions to compute mobility matrix +void DPStokesGlue::setPositions(const real *h_pos, int numberParticles) { + throwIfInvalid(); + dpstokes->setPositions(h_pos, numberParticles); +} - void DPStokesGlue::throwIfInvalid(){ - if(not dpstokes){ - throw std::runtime_error("DPStokes is not initialized. Call Initialize first"); - } +// Get the positions of the particles in the system +const real *DPStokesGlue::getStoredPositions() { + throwIfInvalid(); + return dpstokes->getStoredPositions(); +} + +// Compute the dot product of the mobility matrix with the forces and/or torques +// acting on the previously provided positions +void DPStokesGlue::Mdot(const real *h_forces, const real *h_torques, real *h_MF, + real *h_MT, int numberParticles, bool includeAngular) { + throwIfInvalid(); + dpstokes->Mdot(h_forces, h_torques, h_MF, h_MT, numberParticles, + includeAngular); +} + +void DPStokesGlue::throwIfInvalid() { + if (not dpstokes) { + throw std::runtime_error( + "DPStokes is not initialized. Call Initialize first"); } +} - std::string getPrecision() { +std::string getPrecision() { #ifndef DOUBLE_PRECISION - return "single"; + return "single"; #else - return "double"; + return "double"; #endif - } +} - struct precision_type{ +struct precision_type { #ifndef DOUBLE_PRECISION - using type = float; + using type = float; #else - using type = double; + using type = double; #endif - }; +}; -} +} // namespace uammd_dpstokes diff --git a/solvers/DPStokes/python_wrapper.cu b/solvers/DPStokes/python_wrapper.cu index 379e67f2..ce0013aa 100644 --- a/solvers/DPStokes/python_wrapper.cu +++ b/solvers/DPStokes/python_wrapper.cu @@ -22,20 +22,19 @@ allowChangingBoxSize : bool Whether the periodic extents Lx & Ly can be modified during parameter selection. Default: false. )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); - , docstring); +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); + , docstring); diff --git a/solvers/NBody/example.py b/solvers/NBody/example.py deleted file mode 100644 index 62e83a3c..00000000 --- a/solvers/NBody/example.py +++ /dev/null @@ -1,3 +0,0 @@ -// Donev: This seems too trivial to be useful -import NBody -help(NBody) diff --git a/solvers/NBody/extra/NbodyRPY.cu b/solvers/NBody/extra/NbodyRPY.cu index 24dfb31e..df12b81f 100644 --- a/solvers/NBody/extra/NbodyRPY.cu +++ b/solvers/NBody/extra/NbodyRPY.cu @@ -23,12 +23,11 @@ namespace nbody_rpy { // computation. Threads will tipically read one value from global memory but // blockDim.x from shared memory. template -__global__ void computeRPYBatchedFastGPU(const vecType *pos, - const vecType *forces, - const vecType *torques, real3 *Mv, - real3 *Mw, Constants constants, - blocksToCompute blocks, - int Nbatches, int NperBatch) { +__global__ void +computeRPYBatchedFastGPU(const vecType *pos, const vecType *forces, + const vecType *torques, real3 *Mv, real3 *Mw, + Constants constants, blocksToCompute blocks, + int Nbatches, int NperBatch) { const int tid = blockIdx.x * blockDim.x + threadIdx.x; const int N = Nbatches * NperBatch; const bool active = tid < N; @@ -79,7 +78,8 @@ __global__ void computeRPYBatchedFastGPU(const vecType *pos, const real3 pj = shPos[counter]; const real3 fj = forces ? shForce[counter] : real3(); const real3 tj = torques ? shTorque[counter] : real3(); - const mdot_result dot = HydrodynamicKernel::dotProduct(pi, pj, fj, tj, constants, blocks); + const mdot_result dot = + HydrodynamicKernel::dotProduct(pi, pj, fj, tj, constants, blocks); MF += dot.MF; MT += dot.MT; } @@ -88,7 +88,7 @@ __global__ void computeRPYBatchedFastGPU(const vecType *pos, __syncthreads(); } if (active) { - if(Mv) + if (Mv) Mv[id] = MF; if (Mw) Mw[id] = MT; @@ -98,14 +98,14 @@ __global__ void computeRPYBatchedFastGPU(const vecType *pos, template void computeRPYBatchedFast(const vecType *pos, const vecType *force, const vecType *torque, real3 *Mv, real3 *Mw, - int Nbatches, int NperBatch, - Constants &constants, blocksToCompute &blocks) { + int Nbatches, int NperBatch, Constants &constants, + blocksToCompute &blocks) { int N = Nbatches * NperBatch; int nearestWarpMultiple = ((NperBatch + 16) / 32) * 32; int minBlockSize = std::max(nearestWarpMultiple, 32); int Nthreads = std::min(std::min(minBlockSize, N), 256); int Nblocks = (N + Nthreads - 1) / Nthreads; - int sharedMemoryFactor = 1; // always allocate shared memory for positions + int sharedMemoryFactor = 1; // always allocate shared memory for positions // allocate shared memory for forces & torques separately if provided sharedMemoryFactor += force != nullptr; sharedMemoryFactor += torque != nullptr; @@ -116,12 +116,11 @@ void computeRPYBatchedFast(const vecType *pos, const vecType *force, template // Naive N^2 algorithm (looks like x20 times slower than the fast kernel -__global__ void computeRPYBatchedNaiveGPU(const vecType *pos, - const vecType *forces, - const vecType *torques, real3 *Mv, - real3 *Mw, Constants constants, - blocksToCompute blocks, - int Nbatches, int NperBatch) { +__global__ void +computeRPYBatchedNaiveGPU(const vecType *pos, const vecType *forces, + const vecType *torques, real3 *Mv, real3 *Mw, + Constants constants, blocksToCompute blocks, + int Nbatches, int NperBatch) { const int tid = blockIdx.x * blockDim.x + threadIdx.x; if (tid >= Nbatches * NperBatch) return; @@ -140,8 +139,10 @@ __global__ void computeRPYBatchedNaiveGPU(const vecType *pos, MF += dot.MF; MT += dot.MT; } - if (Mv) Mv[tid] = MF; - if (Mw) Mw[tid] = MT; + if (Mv) + Mv[tid] = MF; + if (Mw) + Mw[tid] = MT; } template @@ -159,10 +160,11 @@ void computeRPYBatchedNaive(const vecType *pos, const vecType *force, template // NaiveBlock N^2 algorithm (looks like x20 times slower than the fast kernel -__global__ void computeRPYBatchedNaiveBlockGPU( - const vecType *pos, const vecType *forces, const vecType *torque, real3 *Mv, - real3 *Mw, Constants constants, blocksToCompute blocks, int Nbatches, - int NperBatch) { +__global__ void +computeRPYBatchedNaiveBlockGPU(const vecType *pos, const vecType *forces, + const vecType *torque, real3 *Mv, real3 *Mw, + Constants constants, blocksToCompute blocks, + int Nbatches, int NperBatch) { const int tid = blockIdx.x; if (tid >= Nbatches * NperBatch) return; @@ -178,22 +180,29 @@ __global__ void computeRPYBatchedNaiveBlockGPU( // real3 fj = make_real3(forces[i]); const real3 fj = forces ? make_real3(forces[i]) : real3(); const real3 tj = torque ? make_real3(torque[i]) : real3(); - mdot_result dot = HydrodynamicKernel::dotProduct(pi, pj, fj, tj, constants, blocks); + mdot_result dot = + HydrodynamicKernel::dotProduct(pi, pj, fj, tj, constants, blocks); MF += dot.MF; MT += dot.MT; } - if (Mv) sharedMemory[threadIdx.x] = MF; - if (Mw) sharedMemory[threadIdx.x + blockDim.x] = MT; + if (Mv) + sharedMemory[threadIdx.x] = MF; + if (Mw) + sharedMemory[threadIdx.x + blockDim.x] = MT; __syncthreads(); if (threadIdx.x == 0) { auto MFTot = real3(); auto MTTot = real3(); for (int i = 0; i < blockDim.x; i++) { - if (Mv) MFTot += sharedMemory[i]; - if (Mw) MTTot += sharedMemory[i + blockDim.x]; + if (Mv) + MFTot += sharedMemory[i]; + if (Mw) + MTTot += sharedMemory[i + blockDim.x]; } - if (Mv) Mv[tid] = MFTot; - if (Mw) Mw[tid] = MTTot; + if (Mv) + Mv[tid] = MFTot; + if (Mw) + Mw[tid] = MTTot; } } @@ -219,7 +228,8 @@ template void batchedNBody(device_span ipos, device_span iforces, device_span itorques, device_span iMF, device_span iMT, int Nbatches, int NperBatch, - Constants& constants, blocksToCompute& blocks, algorithm alg) { + Constants &constants, blocksToCompute &blocks, + algorithm alg) { if (ipos.size() < Nbatches * NperBatch * 3) throw std::runtime_error("Not enough space in pos"); device_adapter pos(ipos, device::cuda); @@ -264,12 +274,12 @@ void callBatchedNBody(device_span pos, } if (kernel == kernel_type::bottom_wall) { - batchedNBody( - pos, forces, torques, MF, MT, Nbatches, NperBatch, constants, blocks, alg); + batchedNBody(pos, forces, torques, MF, MT, Nbatches, NperBatch, + constants, blocks, alg); } else if (kernel == kernel_type::open_rpy) { - batchedNBody( - pos, forces, torques, MF, MT, Nbatches, NperBatch, constants, blocks, alg); + batchedNBody(pos, forces, torques, MF, MT, Nbatches, + NperBatch, constants, blocks, alg); } else { throw std::runtime_error("Unknown kernel type"); } diff --git a/solvers/NBody/extra/hydrodynamicKernels.cuh b/solvers/NBody/extra/hydrodynamicKernels.cuh index f83dc962..ffcfe8dc 100644 --- a/solvers/NBody/extra/hydrodynamicKernels.cuh +++ b/solvers/NBody/extra/hydrodynamicKernels.cuh @@ -16,27 +16,25 @@ and angular forces (torques). For example, dotProduct_UF computes U = MF. namespace nbody_rpy { - // constants for each block of the mobility tensor and the hydrodynamic radius - struct Constants - { - const real t0; // trans-trans - const real r0; // rot-rot - const real rt0; //rot-trans and trans-rot - const real rh; // hydro radius - }; - - struct blocksToCompute{ - bool useUF = false; - bool useWT = false; - bool useUT = false; - bool useWF = false; - }; - - struct mdot_result - { - real3 MF = real3(0, 0, 0); - real3 MT = real3(0, 0, 0); - }; +// constants for each block of the mobility tensor and the hydrodynamic radius +struct Constants { + const real t0; // trans-trans + const real r0; // rot-rot + const real rt0; // rot-trans and trans-rot + const real rh; // hydro radius +}; + +struct blocksToCompute { + bool useUF = false; + bool useWT = false; + bool useUT = false; + bool useWF = false; +}; + +struct mdot_result { + real3 MF = real3(0, 0, 0); + real3 MT = real3(0, 0, 0); +}; // RPY = (1/(6*pi*viscosity*rh))*(f*I + g* r\diadic r/r^2). rh is hydrodynamic // radius. This function returns {f, g/r^2} @@ -123,58 +121,55 @@ __device__ real3 RPY_WF(real3 rij, real r, real rh) { class OpenBoundary { public: - static __device__ mdot_result dotProduct(real3 pi, real3 pj, real3 fj, real3 tj, Constants &c, blocksToCompute &blocks) - { - mdot_result result; - real3 rij = make_real3(pi) - make_real3(pj); - const real r = sqrt(dot(rij, rij)); - - if (blocks.useUF) - result.MF += c.t0 * dotProduct_UF(rij, r, fj, c.rh); - if (blocks.useUT) - result.MF += c.rt0 * dotProduct_UT(rij, r, tj, c.rh); - if (blocks.useWT) - result.MT += c.r0 * dotProduct_WT(rij, r, tj, c.rh); - if (blocks.useWF) - result.MT += c.rt0 * dotProduct_WF(rij, r, fj, c.rh); - return result; - } + static __device__ mdot_result dotProduct(real3 pi, real3 pj, real3 fj, + real3 tj, Constants &c, + blocksToCompute &blocks) { + mdot_result result; + real3 rij = make_real3(pi) - make_real3(pj); + const real r = sqrt(dot(rij, rij)); + + if (blocks.useUF) + result.MF += c.t0 * dotProduct_UF(rij, r, fj, c.rh); + if (blocks.useUT) + result.MF += c.rt0 * dotProduct_UT(rij, r, tj, c.rh); + if (blocks.useWT) + result.MT += c.r0 * dotProduct_WT(rij, r, tj, c.rh); + if (blocks.useWF) + result.MT += c.rt0 * dotProduct_WF(rij, r, fj, c.rh); + return result; + } // Computes M(ri, rj)*vj - static __device__ real3 dotProduct_UF(real3 rij, real r, real3 fj, real rh) - { - const real2 c12 = RPY_UF(r, rh); - const real f = c12.x; - const real gdivr2 = c12.y; - const real gv = gdivr2 * dot(rij, fj); - const real3 Mv_t = f * fj + (r > real(0) ? gv * rij : real3()); - return Mv_t; + static __device__ real3 dotProduct_UF(real3 rij, real r, real3 fj, real rh) { + const real2 c12 = RPY_UF(r, rh); + const real f = c12.x; + const real gdivr2 = c12.y; + const real gv = gdivr2 * dot(rij, fj); + const real3 Mv_t = f * fj + (r > real(0) ? gv * rij : real3()); + return Mv_t; } - static __device__ real3 dotProduct_WT(real3 rij, real r, real3 tj, real rh) - { - const real2 c12 = RPY_WT(r, rh); - const real f = c12.x; - const real gdivr2 = c12.y; - const real gv = gdivr2 * dot(rij, tj); - const real3 Mv_t = f * tj + (r > real(0) ? gv * rij : real3()); - return Mv_t; + static __device__ real3 dotProduct_WT(real3 rij, real r, real3 tj, real rh) { + const real2 c12 = RPY_WT(r, rh); + const real f = c12.x; + const real gdivr2 = c12.y; + const real gv = gdivr2 * dot(rij, tj); + const real3 Mv_t = f * tj + (r > real(0) ? gv * rij : real3()); + return Mv_t; } - static __device__ real3 dotProduct_UT(real3 rij, real r, real3 tj, real rh) - { - const real3 m = RPY_UT(rij, r, rh); // (M_xy, M_xz, M_yz) - const real3 Mv_t = {m.x * tj.y + m.y * tj.z, -m.x * tj.x + m.z * tj.z, - -m.y * tj.x - m.z * tj.y}; - return Mv_t; + static __device__ real3 dotProduct_UT(real3 rij, real r, real3 tj, real rh) { + const real3 m = RPY_UT(rij, r, rh); // (M_xy, M_xz, M_yz) + const real3 Mv_t = {m.x * tj.y + m.y * tj.z, -m.x * tj.x + m.z * tj.z, + -m.y * tj.x - m.z * tj.y}; + return Mv_t; } - static __device__ real3 dotProduct_WF(real3 rij, real r, real3 fj, real rh) - { - const real3 m = RPY_WF(rij, r, rh); // (M_xy, M_xz, M_yz) - const real3 Mv_t = {m.x * fj.y + m.y * fj.z, -m.x * fj.x + m.z * fj.z, - -m.y * fj.x - m.z * fj.y}; - return Mv_t; + static __device__ real3 dotProduct_WF(real3 rij, real r, real3 fj, real rh) { + const real3 m = RPY_WF(rij, r, rh); // (M_xy, M_xz, M_yz) + const real3 Mv_t = {m.x * fj.y + m.y * fj.z, -m.x * fj.x + m.z * fj.z, + -m.y * fj.x - m.z * fj.y}; + return Mv_t; } }; @@ -189,7 +184,8 @@ class BottomWall { // located at z=0 rij: distance between particles rij.z: This component // contains ((pi.z-pj.z) + 2*pj.z)/rh self: self interaction hj: height of the // particle j vj: quantity (i.e force) of particle j - static __device__ real3 wallCorrection_UF(real3 rij, bool self, real hj, real3 vj) { + static __device__ real3 wallCorrection_UF(real3 rij, bool self, real hj, + real3 vj) { real3 correction = real3(0, 0, 0); if (self) { // B1*vj in [1] real invZi = real(1.0) / hj; @@ -254,7 +250,8 @@ class BottomWall { // NOTE: normalized by 8 pi eta a**3. [1] normalizes by 6 pi et a**3 // so, all coeffs are multiplied by 4/3 if compared to [1] - static __device__ real3 wallCorrection_WT(real3 rij, bool self, real hj, real3 vj) { + static __device__ real3 wallCorrection_WT(real3 rij, bool self, real hj, + real3 vj) { real3 correction = real3(0, 0, 0); if (self) { // B3*vj in [1] real invZi = real(1.0) / hj; @@ -287,7 +284,8 @@ class BottomWall { // note: [1] seemingly uses the left-hand rule for the cross product, so the // signs are flipped on expressions in that paper that use the Levi-Civita // symbol - static __device__ real3 wallCorrection_UT(real3 rij, bool self, real h, real3 vj) { + static __device__ real3 wallCorrection_UT(real3 rij, bool self, real h, + real3 vj) { real3 correction = real3(0, 0, 0); if (self) { // B2^T*vj in [1]. ^T denotes transpose. real invZi = real(1.0) / h; @@ -323,7 +321,8 @@ class BottomWall { // note: [1] seemingly uses the left-hand rule for the cross product, so the // signs are flipped on expressions in that paper that use the Levi-Civita // symbol - static __device__ real3 wallCorrection_WF(real3 rij, bool self, real h, real3 vj) { + static __device__ real3 wallCorrection_WF(real3 rij, bool self, real h, + real3 vj) { real3 correction = real3(0, 0, 0); if (self) { // B2*fj in [1]. real invZi = real(1.0) / h; @@ -360,37 +359,38 @@ class BottomWall { } public: - - static __device__ mdot_result dotProduct(real3 pi, real3 pj, real3 fj, real3 tj, Constants &c, blocksToCompute &blocks) - { - // implements damping from Appendix 1 in [2] so the matrix is positive - // definite when a particle overlaps the wall - mdot_result result; - real bi = min(pi.z / c.rh, real(1.0)); - bi = max(bi, real(0.0)); - real bj = min(pj.z / c.rh, real(1.0)); - bj = max(bj, real(0.0)); - real bij = bi * bj; - - pi.z = max(pi.z, c.rh); - pj.z = max(pj.z, c.rh); - - real3 rij = make_real3(pi) - make_real3(pj); - const real r = sqrt(dot(rij, rij)); - - if (blocks.useUF) - result.MF += bij * c.t0 * dotProduct_UF(rij, r, fj, pj.z, c.rh); - if (blocks.useUT) - result.MF += bij * c.rt0 * dotProduct_UT(rij, r, tj, pi.z, c.rh); - if (blocks.useWT) - result.MT += bij * c.r0 * dotProduct_WT(rij, r, tj, pj.z, c.rh); - if (blocks.useWF) - result.MT += bij * c.rt0 * dotProduct_WF(rij, r, fj, pj.z, c.rh); + static __device__ mdot_result dotProduct(real3 pi, real3 pj, real3 fj, + real3 tj, Constants &c, + blocksToCompute &blocks) { + // implements damping from Appendix 1 in [2] so the matrix is positive + // definite when a particle overlaps the wall + mdot_result result; + real bi = min(pi.z / c.rh, real(1.0)); + bi = max(bi, real(0.0)); + real bj = min(pj.z / c.rh, real(1.0)); + bj = max(bj, real(0.0)); + real bij = bi * bj; + + pi.z = max(pi.z, c.rh); + pj.z = max(pj.z, c.rh); + + real3 rij = make_real3(pi) - make_real3(pj); + const real r = sqrt(dot(rij, rij)); + + if (blocks.useUF) + result.MF += bij * c.t0 * dotProduct_UF(rij, r, fj, pj.z, c.rh); + if (blocks.useUT) + result.MF += bij * c.rt0 * dotProduct_UT(rij, r, tj, pi.z, c.rh); + if (blocks.useWT) + result.MT += bij * c.r0 * dotProduct_WT(rij, r, tj, pj.z, c.rh); + if (blocks.useWF) + result.MT += bij * c.rt0 * dotProduct_WF(rij, r, fj, pj.z, c.rh); return result; } - static __device__ real3 dotProduct_UF(real3 rij, real r, real3 fj, real hj, real rh) { + static __device__ real3 dotProduct_UF(real3 rij, real r, real3 fj, real hj, + real rh) { const real2 c12 = RPY_UF(r, rh); const real f = c12.x; const real gdivr2 = c12.y; @@ -401,7 +401,8 @@ public: return Mv_t; } - static __device__ real3 dotProduct_WT(real3 rij, real r, real3 tj, real hj, real rh) { + static __device__ real3 dotProduct_WT(real3 rij, real r, real3 tj, real hj, + real rh) { const real2 c12 = RPY_WT(r, rh); const real f = c12.x; const real gdivr2 = c12.y; @@ -418,7 +419,8 @@ public: // calling loop. so, we call wallCorrection_UT with R = -rij = pj - pi and h = // pi.z, i.e. we flip the order of the (ij) arguments so that we get M_{UT, // ij} - static __device__ real3 dotProduct_UT(real3 rij, real r, real3 tj, real hi, real rh) { + static __device__ real3 dotProduct_UT(real3 rij, real r, real3 tj, real hi, + real rh) { const real3 m = RPY_UT(rij, r, rh); // (M_xy, M_xz, M_yz) real3 Mv_t = {m.x * tj.y + m.y * tj.z, -m.x * tj.x + m.z * tj.z, -m.y * tj.x - m.z * tj.y}; @@ -428,7 +430,8 @@ public: return Mv_t; } - static __device__ real3 dotProduct_WF(real3 rij, real r, real3 fj, real hj, real rh) { + static __device__ real3 dotProduct_WF(real3 rij, real r, real3 fj, real hj, + real rh) { const real3 m = RPY_WF(rij, r, rh); // (M_xy, M_xz, M_yz) real3 Mv_t = {m.x * fj.y + m.y * fj.z, -m.x * fj.x + m.z * fj.z, -m.y * fj.x - m.z * fj.y}; diff --git a/solvers/NBody/extra/interface.h b/solvers/NBody/extra/interface.h index bf486370..8b48972c 100644 --- a/solvers/NBody/extra/interface.h +++ b/solvers/NBody/extra/interface.h @@ -1,8 +1,8 @@ #ifndef NBODYRPY_INTERFACE_H #define NBODYRPY_INTERFACE_H -#include "memory/container.h" #include "MobilityInterface/defines.h" +#include "memory/container.h" namespace nbody_rpy { enum class kernel_type { open_rpy, bottom_wall }; using namespace libmobility; @@ -13,8 +13,8 @@ void callBatchedNBody(device_span pos, device_span torques, device_span MF, device_span MT, int Nbatches, int NperBatch, real transMobility, real rotMobility, - real transRotMobility, real hydrodynamicRadius, bool needsTorque, - algorithm alg, kernel_type kernel); + real transRotMobility, real hydrodynamicRadius, + bool needsTorque, algorithm alg, kernel_type kernel); } // namespace nbody_rpy #endif diff --git a/solvers/NBody/extra/vector.cuh b/solvers/NBody/extra/vector.cuh index 5d4b3cbe..7d0e7d68 100644 --- a/solvers/NBody/extra/vector.cuh +++ b/solvers/NBody/extra/vector.cuh @@ -3,23 +3,23 @@ */ #ifndef VECTOR_OVERLOADS_H #define VECTOR_OVERLOADS_H -//Include built in ones +// Include built in ones #include "cuda_runtime.h" -#include #include +#include -namespace nbody_rpy{ - using real = libmobility::real; +namespace nbody_rpy { +using real = libmobility::real; #ifdef DOUBLE_PRECISION - using real2 = double2; - using real3 = double3; - using real4 = double4; +using real2 = double2; +using real3 = double3; +using real4 = double4; #else - using real2 = float2; - using real3 = float3; - using real4 = float4; +using real2 = float2; +using real3 = float3; +using real4 = float4; #endif -} +} // namespace nbody_rpy typedef unsigned short ushort; typedef unsigned int uint; @@ -29,920 +29,763 @@ typedef unsigned long long int ullint; /////////////////////FLOAT2/////////////////////////////// -VECATTR int2 make_int2(float2 a){return make_int2((int)a.x, (int)a.y);} -VECATTR float2 make_float2(float a){return make_float2(a, a);} +VECATTR int2 make_int2(float2 a) { return make_int2((int)a.x, (int)a.y); } +VECATTR float2 make_float2(float a) { return make_float2(a, a); } -VECATTR float2 make_float2(int2 a){return make_float2(a.x, a.y);} -VECATTR float2 make_float2(float2 a){return make_float2(a.x, a.y);} +VECATTR float2 make_float2(int2 a) { return make_float2(a.x, a.y); } +VECATTR float2 make_float2(float2 a) { return make_float2(a.x, a.y); } -VECATTR float2 operator +(const float2 &a, const float2 &b){return make_float2( - a.x + b.x, - a.y + b.y); +VECATTR float2 operator+(const float2 &a, const float2 &b) { + return make_float2(a.x + b.x, a.y + b.y); } -VECATTR void operator +=(float2 &a, const float2 &b){ +VECATTR void operator+=(float2 &a, const float2 &b) { a.x += b.x; a.y += b.y; } -VECATTR float2 operator +(const float2 &a, const float &b){return make_float2( - a.x + b, - a.y + b); +VECATTR float2 operator+(const float2 &a, const float &b) { + return make_float2(a.x + b, a.y + b); } -VECATTR float2 operator +(const float &b, const float2 &a){return a+b;} -VECATTR void operator +=(float2 &a, const float &b){ +VECATTR float2 operator+(const float &b, const float2 &a) { return a + b; } +VECATTR void operator+=(float2 &a, const float &b) { a.x += b; a.y += b; } -VECATTR float2 operator -(const float2 &a, const float2 &b){return make_float2( - a.x - b.x, - a.y - b.y); +VECATTR float2 operator-(const float2 &a, const float2 &b) { + return make_float2(a.x - b.x, a.y - b.y); } -VECATTR void operator -=(float2 &a, const float2 &b){ +VECATTR void operator-=(float2 &a, const float2 &b) { a.x -= b.x; a.y -= b.y; } -VECATTR float2 operator -(const float2 &a, const float &b){return make_float2( - a.x - b, - a.y - b); +VECATTR float2 operator-(const float2 &a, const float &b) { + return make_float2(a.x - b, a.y - b); } -VECATTR float2 operator -(const float &b, const float2 &a){return make_float2( - b-a.x, - b-a.y); +VECATTR float2 operator-(const float &b, const float2 &a) { + return make_float2(b - a.x, b - a.y); } -VECATTR void operator -=(float2 &a, const float &b){ +VECATTR void operator-=(float2 &a, const float &b) { a.x -= b; a.y -= b; } -VECATTR float2 operator *(const float2 &a, const float2 &b){ - return make_float2(a.x * b.x, - a.y * b.y); +VECATTR float2 operator*(const float2 &a, const float2 &b) { + return make_float2(a.x * b.x, a.y * b.y); } -VECATTR void operator *=(float2 &a, const float2 &b){ +VECATTR void operator*=(float2 &a, const float2 &b) { a.x *= b.x; a.y *= b.y; } -VECATTR float2 operator *(const float2 &a, const float &b){ - return make_float2(a.x * b, - a.y * b); +VECATTR float2 operator*(const float2 &a, const float &b) { + return make_float2(a.x * b, a.y * b); } -VECATTR float2 operator *(const float &b, const float2 &a){ - return make_float2(a.x * b, - a.y * b); +VECATTR float2 operator*(const float &b, const float2 &a) { + return make_float2(a.x * b, a.y * b); } -VECATTR void operator *=(float2 &a, const float &b){ +VECATTR void operator*=(float2 &a, const float &b) { a.x *= b; a.y *= b; } -VECATTR float2 operator /(const float2 &a, const float2 &b){ - return make_float2(a.x / b.x, - a.y / b.y); +VECATTR float2 operator/(const float2 &a, const float2 &b) { + return make_float2(a.x / b.x, a.y / b.y); } -VECATTR void operator /=(float2 &a, const float2 &b){ +VECATTR void operator/=(float2 &a, const float2 &b) { a.x /= b.x; a.y /= b.y; } -VECATTR float2 operator /(const float2 &a, const float &b){ - return (1.0f/b)*a; -} -VECATTR float2 operator /(const float &b, const float2 &a){ - return make_float2(b / a.x, - b / a.y); +VECATTR float2 operator/(const float2 &a, const float &b) { + return (1.0f / b) * a; } -VECATTR void operator /=(float2 &a, const float &b){ - a *= 1.0f/b; +VECATTR float2 operator/(const float &b, const float2 &a) { + return make_float2(b / a.x, b / a.y); } - - +VECATTR void operator/=(float2 &a, const float &b) { a *= 1.0f / b; } /////////////////////FLOAT3/////////////////////////////// -VECATTR int3 make_int3(float3 a){return make_int3((int)a.x, (int)a.y, (int)a.z);} -VECATTR float3 make_float3(float a){return make_float3(a, a, a);} +VECATTR int3 make_int3(float3 a) { + return make_int3((int)a.x, (int)a.y, (int)a.z); +} +VECATTR float3 make_float3(float a) { return make_float3(a, a, a); } -VECATTR float3 make_float3(int3 a){return make_float3(a.x, a.y, a.z);} -VECATTR float3 make_float3(float3 a){return make_float3(a.x, a.y, a.z);} -VECATTR float3 make_float3(float4 a){return make_float3(a.x, a.y, a.z);} +VECATTR float3 make_float3(int3 a) { return make_float3(a.x, a.y, a.z); } +VECATTR float3 make_float3(float3 a) { return make_float3(a.x, a.y, a.z); } +VECATTR float3 make_float3(float4 a) { return make_float3(a.x, a.y, a.z); } -VECATTR float3 operator +(const float3 &a, const float3 &b){return make_float3( - a.x + b.x, - a.y + b.y, - a.z + b.z); +VECATTR float3 operator+(const float3 &a, const float3 &b) { + return make_float3(a.x + b.x, a.y + b.y, a.z + b.z); } -VECATTR void operator +=(float3 &a, const float3 &b){ +VECATTR void operator+=(float3 &a, const float3 &b) { a.x += b.x; a.y += b.y; a.z += b.z; } -VECATTR float3 operator +(const float3 &a, const float &b){return make_float3( - a.x + b, - a.y + b, - a.z + b); +VECATTR float3 operator+(const float3 &a, const float &b) { + return make_float3(a.x + b, a.y + b, a.z + b); } -VECATTR float3 operator +(const float &b, const float3 &a){return a+b;} -VECATTR void operator +=(float3 &a, const float &b){ +VECATTR float3 operator+(const float &b, const float3 &a) { return a + b; } +VECATTR void operator+=(float3 &a, const float &b) { a.x += b; a.y += b; a.z += b; } -VECATTR float3 operator -(const float3 &a, const float3 &b){return make_float3( - a.x - b.x, - a.y - b.y, - a.z - b.z); +VECATTR float3 operator-(const float3 &a, const float3 &b) { + return make_float3(a.x - b.x, a.y - b.y, a.z - b.z); } -VECATTR void operator -=(float3 &a, const float3 &b){ +VECATTR void operator-=(float3 &a, const float3 &b) { a.x -= b.x; a.y -= b.y; a.z -= b.z; } -VECATTR float3 operator -(const float3 &a, const float &b){return make_float3( - a.x - b, - a.y - b, - a.z - b); +VECATTR float3 operator-(const float3 &a, const float &b) { + return make_float3(a.x - b, a.y - b, a.z - b); } -VECATTR float3 operator -(const float &b, const float3 &a){return make_float3( - b-a.x, - b-a.y, - b-a.z); +VECATTR float3 operator-(const float &b, const float3 &a) { + return make_float3(b - a.x, b - a.y, b - a.z); } -VECATTR void operator -=(float3 &a, const float &b){ +VECATTR void operator-=(float3 &a, const float &b) { a.x -= b; a.y -= b; a.z -= b; } -VECATTR float3 operator *(const float3 &a, const float3 &b){ - return make_float3( - a.x * b.x, - a.y * b.y, - a.z * b.z - ); +VECATTR float3 operator*(const float3 &a, const float3 &b) { + return make_float3(a.x * b.x, a.y * b.y, a.z * b.z); } -VECATTR void operator *=(float3 &a, const float3 &b){ +VECATTR void operator*=(float3 &a, const float3 &b) { a.x *= b.x; a.y *= b.y; a.z *= b.z; } -VECATTR float3 operator *(const float3 &a, const float &b){ - return make_float3( - a.x * b, - a.y * b, - a.z * b - ); -} -VECATTR float3 operator *(const float &b, const float3 &a){ - return make_float3( - a.x * b, - a.y * b, - a.z * b - ); -} -VECATTR void operator *=(float3 &a, const float &b){ +VECATTR float3 operator*(const float3 &a, const float &b) { + return make_float3(a.x * b, a.y * b, a.z * b); +} +VECATTR float3 operator*(const float &b, const float3 &a) { + return make_float3(a.x * b, a.y * b, a.z * b); +} +VECATTR void operator*=(float3 &a, const float &b) { a.x *= b; a.y *= b; a.z *= b; } -VECATTR float3 operator /(const float3 &a, const float3 &b){ - return make_float3( - a.x / b.x, - a.y / b.y, - a.z / b.z - ); +VECATTR float3 operator/(const float3 &a, const float3 &b) { + return make_float3(a.x / b.x, a.y / b.y, a.z / b.z); } -VECATTR void operator /=(float3 &a, const float3 &b){ +VECATTR void operator/=(float3 &a, const float3 &b) { a.x /= b.x; a.y /= b.y; a.z /= b.z; } -VECATTR float3 operator /(const float3 &a, const float &b){ - return (1.0f/b)*a; -} -VECATTR float3 operator /(const float &b, const float3 &a){ - return make_float3( - b / a.x, - b / a.y, - b / a.z - ); +VECATTR float3 operator/(const float3 &a, const float &b) { + return (1.0f / b) * a; } -VECATTR void operator /=(float3 &a, const float &b){ - a *= 1.0f/b; +VECATTR float3 operator/(const float &b, const float3 &a) { + return make_float3(b / a.x, b / a.y, b / a.z); } +VECATTR void operator/=(float3 &a, const float &b) { a *= 1.0f / b; } - -VECATTR float3 floorf(const float3 &a){return make_float3(floorf(a.x), floorf(a.y), floorf(a.z));} +VECATTR float3 floorf(const float3 &a) { + return make_float3(floorf(a.x), floorf(a.y), floorf(a.z)); +} /////////////////////FLOAT4/////////////////////////////// +VECATTR float4 make_float4(float a) { return make_float4(a, a, a, a); } -VECATTR float4 make_float4(float a){return make_float4(a,a,a,a);} - -VECATTR float4 make_float4(float3 a){return make_float4(a.x, a.y, a.z, 0);} -VECATTR float4 make_float4(float4 a){return make_float4(a.x, a.y, a.z, a.w);} +VECATTR float4 make_float4(float3 a) { return make_float4(a.x, a.y, a.z, 0); } +VECATTR float4 make_float4(float4 a) { return make_float4(a.x, a.y, a.z, a.w); } -VECATTR float4 operator +(const float4 &a, const float4 &b){return make_float4( - a.x + b.x, - a.y + b.y, - a.z + b.z, - a.w + b.w); +VECATTR float4 operator+(const float4 &a, const float4 &b) { + return make_float4(a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w); } -VECATTR void operator +=(float4 &a, const float4 &b){ +VECATTR void operator+=(float4 &a, const float4 &b) { a.x += b.x; a.y += b.y; a.z += b.z; a.w += b.w; } -VECATTR float4 operator +(const float4 &a, const float &b){return make_float4( - a.x + b, - a.y + b, - a.z + b, - a.w + b); +VECATTR float4 operator+(const float4 &a, const float &b) { + return make_float4(a.x + b, a.y + b, a.z + b, a.w + b); } -VECATTR float4 operator +(const float &b, const float4 &a){return a+b;} -VECATTR void operator +=(float4 &a, const float &b){ +VECATTR float4 operator+(const float &b, const float4 &a) { return a + b; } +VECATTR void operator+=(float4 &a, const float &b) { a.x += b; a.y += b; a.z += b; a.w += b; } -VECATTR float4 operator -(const float4 &a, const float4 &b){return make_float4( - a.x - b.x, - a.y - b.y, - a.z - b.z, - a.w - b.w); +VECATTR float4 operator-(const float4 &a, const float4 &b) { + return make_float4(a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w); } -VECATTR void operator -=(float4 &a, const float4 &b){ +VECATTR void operator-=(float4 &a, const float4 &b) { a.x -= b.x; a.y -= b.y; a.z -= b.z; a.w -= b.w; } -VECATTR float4 operator -(const float4 &a, const float &b){return make_float4( - a.x - b, - a.y - b, - a.z - b, - a.w - b); +VECATTR float4 operator-(const float4 &a, const float &b) { + return make_float4(a.x - b, a.y - b, a.z - b, a.w - b); } -VECATTR float4 operator -(const float &b, const float4 &a){return make_float4( - b-a.x, - b-a.y, - b-a.z, - b-a.w); +VECATTR float4 operator-(const float &b, const float4 &a) { + return make_float4(b - a.x, b - a.y, b - a.z, b - a.w); } -VECATTR void operator -=(float4 &a, const float &b){ +VECATTR void operator-=(float4 &a, const float &b) { a.x -= b; a.y -= b; a.z -= b; a.w -= b; } -VECATTR float4 operator *(const float4 &a, const float4 &b){ - return make_float4(a.x * b.x, - a.y * b.y, - a.z * b.z, - a.w * b.w); +VECATTR float4 operator*(const float4 &a, const float4 &b) { + return make_float4(a.x * b.x, a.y * b.y, a.z * b.z, a.w * b.w); } -VECATTR void operator *=(float4 &a, const float4 &b){ +VECATTR void operator*=(float4 &a, const float4 &b) { a.x *= b.x; a.y *= b.y; a.z *= b.z; a.w *= b.w; } -VECATTR float4 operator *(const float4 &a, const float &b){ - return make_float4(a.x * b, - a.y * b, - a.z * b, - a.w * b); +VECATTR float4 operator*(const float4 &a, const float &b) { + return make_float4(a.x * b, a.y * b, a.z * b, a.w * b); } -VECATTR float4 operator *(const float &b, const float4 &a){ - return make_float4(a.x * b, - a.y * b, - a.z * b, - a.w * b); +VECATTR float4 operator*(const float &b, const float4 &a) { + return make_float4(a.x * b, a.y * b, a.z * b, a.w * b); } -VECATTR void operator *=(float4 &a, const float &b){ +VECATTR void operator*=(float4 &a, const float &b) { a.x *= b; a.y *= b; a.z *= b; a.w *= b; } -VECATTR float4 operator /(const float4 &a, const float4 &b){ - return make_float4(a.x / b.x, - a.y / b.y, - a.z / b.z, - a.w / b.w); +VECATTR float4 operator/(const float4 &a, const float4 &b) { + return make_float4(a.x / b.x, a.y / b.y, a.z / b.z, a.w / b.w); } -VECATTR void operator /=(float4 &a, const float4 &b){ +VECATTR void operator/=(float4 &a, const float4 &b) { a.x /= b.x; a.y /= b.y; a.z /= b.z; a.w /= b.w; } -VECATTR float4 operator /(const float4 &a, const float &b){ - return (1.0f/b)*a; +VECATTR float4 operator/(const float4 &a, const float &b) { + return (1.0f / b) * a; } -VECATTR float4 operator /(const float &b, const float4 &a){ - return make_float4(b / a.x, - b / a.y, - b / a.z, - b / a.w); +VECATTR float4 operator/(const float &b, const float4 &a) { + return make_float4(b / a.x, b / a.y, b / a.z, b / a.w); } -VECATTR void operator /=(float4 &a, const float &b){ - a *= 1.0f/b; -} - +VECATTR void operator/=(float4 &a, const float &b) { a *= 1.0f / b; } -VECATTR float4 floorf(const float4 &a){ +VECATTR float4 floorf(const float4 &a) { return make_float4(floorf(a.x), floorf(a.y), floorf(a.z), floorf(a.w)); } -VECATTR float dot(float4 a, float4 b){return a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w;} - - +VECATTR float dot(float4 a, float4 b) { + return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w; +} -namespace nbody_rpy{ +namespace nbody_rpy { /////////////////REAL4//////////////////////////////// -VECATTR real4 make_real4(real x, real y, real z, real w){ - #ifdef SINGLE_PRECISION - return make_float4(x,y,z,w); - #else - return make_double4(x,y,z,w); - #endif +VECATTR real4 make_real4(real x, real y, real z, real w) { +#ifdef SINGLE_PRECISION + return make_float4(x, y, z, w); +#else + return make_double4(x, y, z, w); +#endif } -VECATTR real4 make_real4(real s){return make_real4(s, s, s, s);} - VECATTR real4 make_real4(real3 a){ return make_real4(a.x, a.y, a.z, real(0.0));} -VECATTR real4 make_real4(real3 a, real w){ return make_real4(a.x, a.y, a.z, w);} +VECATTR real4 make_real4(real s) { return make_real4(s, s, s, s); } +VECATTR real4 make_real4(real3 a) { + return make_real4(a.x, a.y, a.z, real(0.0)); +} +VECATTR real4 make_real4(real3 a, real w) { + return make_real4(a.x, a.y, a.z, w); +} - VECATTR real4 make_real4(real2 a){ return make_real4(a.x, a.y, real(0.0), real(0.0));} +VECATTR real4 make_real4(real2 a) { + return make_real4(a.x, a.y, real(0.0), real(0.0)); +} #ifdef SINGLE_PRECISION -VECATTR real4 make_real4(double3 a, real w){return make_real4(a.x, a.y, a.z, w);} +VECATTR real4 make_real4(double3 a, real w) { + return make_real4(a.x, a.y, a.z, w); +} #else -VECATTR real4 make_real4(float3 a, real w){ return make_real4(a.x, a.y, a.z, w);} +VECATTR real4 make_real4(float3 a, real w) { + return make_real4(a.x, a.y, a.z, w); +} #endif -VECATTR real4 make_real4(int4 a){ return make_real4(real(a.x), real(a.y), real(a.z), real(a.w));} -VECATTR real4 make_real4(uint4 a){return make_real4(real(a.x), real(a.y), real(a.z), real(a.w));} - +VECATTR real4 make_real4(int4 a) { + return make_real4(real(a.x), real(a.y), real(a.z), real(a.w)); +} +VECATTR real4 make_real4(uint4 a) { + return make_real4(real(a.x), real(a.y), real(a.z), real(a.w)); +} //////////////////REAL3/////////////////////////// - -VECATTR real3 make_real3(real x, real y, real z){ +VECATTR real3 make_real3(real x, real y, real z) { #ifdef SINGLE_PRECISION - return make_float3(x,y,z); + return make_float3(x, y, z); #else - return make_double3(x,y,z); + return make_double3(x, y, z); #endif } -VECATTR real3 make_real3(real s){ return make_real3(s, s, s);} -VECATTR real3 make_real3(real3 a){return make_real3(a.x, a.y, a.z);} +VECATTR real3 make_real3(real s) { return make_real3(s, s, s); } +VECATTR real3 make_real3(real3 a) { return make_real3(a.x, a.y, a.z); } #ifdef SINGLE_PRECISION -VECATTR real3 make_real3(double3 a){return make_real3(a.x, a.y, a.z);} -VECATTR real3 make_real3(double4 a){return make_real3(a.x, a.y, a.z);} +VECATTR real3 make_real3(double3 a) { return make_real3(a.x, a.y, a.z); } +VECATTR real3 make_real3(double4 a) { return make_real3(a.x, a.y, a.z); } #else - template - VECATTR real3 make_real3(float2 a, T b){return make_real3(a.x, a.y, b);} -VECATTR real3 make_real3(float3 a){return make_real3(a.x, a.y, a.z);} -VECATTR real3 make_real3(float4 a){return make_real3(a.x, a.y, a.z);} +template VECATTR real3 make_real3(float2 a, T b) { + return make_real3(a.x, a.y, b); +} +VECATTR real3 make_real3(float3 a) { return make_real3(a.x, a.y, a.z); } +VECATTR real3 make_real3(float4 a) { return make_real3(a.x, a.y, a.z); } #endif -VECATTR real3 make_real3(real4 a){ return make_real3(a.x, a.y, a.z);} - -VECATTR real3 make_real3(real2 a, real z){return make_real3(a.x, a.y, z);} -VECATTR real3 make_real3(int3 a){ return make_real3(real(a.x), real(a.y), real(a.z));} -VECATTR real3 make_real3(uint3 a){return make_real3(real(a.x), real(a.y), real(a.z));} - +VECATTR real3 make_real3(real4 a) { return make_real3(a.x, a.y, a.z); } +VECATTR real3 make_real3(real2 a, real z) { return make_real3(a.x, a.y, z); } +VECATTR real3 make_real3(int3 a) { + return make_real3(real(a.x), real(a.y), real(a.z)); +} +VECATTR real3 make_real3(uint3 a) { + return make_real3(real(a.x), real(a.y), real(a.z)); +} //////////////////REAL2/////////////////////////// - -VECATTR real2 make_real2(real x, real y){ +VECATTR real2 make_real2(real x, real y) { #ifdef SINGLE_PRECISION - return make_float2(x,y); + return make_float2(x, y); #else - return make_double2(x,y); + return make_double2(x, y); #endif } #ifdef SINGLE_PRECISION - VECATTR real2 make_real2(double2 a){return make_real2(a.x, a.y);} +VECATTR real2 make_real2(double2 a) { return make_real2(a.x, a.y); } #else - VECATTR real2 make_real2(float2 a){return make_real2(a.x, a.y);} +VECATTR real2 make_real2(float2 a) { return make_real2(a.x, a.y); } #endif -VECATTR real2 make_real2(real s){ return make_real2(s, s);} -VECATTR real2 make_real2(real2 a){return make_real2(a.x, a.y);} -VECATTR real2 make_real2(real3 a){return make_real2(a.x, a.y);} -VECATTR real2 make_real2(real4 a){return make_real2(a.x, a.y);} -VECATTR real2 make_real2(int2 a){ return make_real2(real(a.x), real(a.y));} -VECATTR real2 make_real2(int3 a){ return make_real2(real(a.x), real(a.y));} -VECATTR real2 make_real2(uint3 a){return make_real2(real(a.x), real(a.y));} - - VECATTR real dot(real2 a, real2 b){ return a.x*b.x + a.y*b.y;} +VECATTR real2 make_real2(real s) { return make_real2(s, s); } +VECATTR real2 make_real2(real2 a) { return make_real2(a.x, a.y); } +VECATTR real2 make_real2(real3 a) { return make_real2(a.x, a.y); } +VECATTR real2 make_real2(real4 a) { return make_real2(a.x, a.y); } +VECATTR real2 make_real2(int2 a) { return make_real2(real(a.x), real(a.y)); } +VECATTR real2 make_real2(int3 a) { return make_real2(real(a.x), real(a.y)); } +VECATTR real2 make_real2(uint3 a) { return make_real2(real(a.x), real(a.y)); } +VECATTR real dot(real2 a, real2 b) { return a.x * b.x + a.y * b.y; } -} +} // namespace nbody_rpy ////////////////DOUBLE PRECISION////////////////////// #ifdef SINGLE_PRECISION -VECATTR double3 make_double3(nbody_rpy::real4 a){return make_double3(a.x, a.y, a.z);} +VECATTR double3 make_double3(nbody_rpy::real4 a) { + return make_double3(a.x, a.y, a.z); +} #else -VECATTR double3 make_double3(nbody_rpy::real3 a){return make_double3(a.x, a.y, a.z);} -VECATTR double3 make_double3(nbody_rpy::real4 a){return make_double3(a.x, a.y, a.z);} +VECATTR double3 make_double3(nbody_rpy::real3 a) { + return make_double3(a.x, a.y, a.z); +} +VECATTR double3 make_double3(nbody_rpy::real4 a) { + return make_double3(a.x, a.y, a.z); +} #endif -VECATTR float4 make_float4(double4 a){return make_float4(float(a.x), float(a.y), float(a.z), float(a.w));} +VECATTR float4 make_float4(double4 a) { + return make_float4(float(a.x), float(a.y), float(a.z), float(a.w)); +} -VECATTR double4 make_double4(double s){ return make_double4(s, s, s, s);} -VECATTR double4 make_double4(double3 a){return make_double4(a.x, a.y, a.z, 0.0f);} -VECATTR double4 make_double4(double3 a, double w){return make_double4(a.x, a.y, a.z, w);} -VECATTR double4 make_double4(int4 a){return make_double4(double(a.x), double(a.y), double(a.z), double(a.w));} -VECATTR double4 make_double4(uint4 a){return make_double4(double(a.x), double(a.y), double(a.z), double(a.w));} -VECATTR double4 make_double4(float4 a){return make_double4(double(a.x), double(a.y), double(a.z), double(a.w));} +VECATTR double4 make_double4(double s) { return make_double4(s, s, s, s); } +VECATTR double4 make_double4(double3 a) { + return make_double4(a.x, a.y, a.z, 0.0f); +} +VECATTR double4 make_double4(double3 a, double w) { + return make_double4(a.x, a.y, a.z, w); +} +VECATTR double4 make_double4(int4 a) { + return make_double4(double(a.x), double(a.y), double(a.z), double(a.w)); +} +VECATTR double4 make_double4(uint4 a) { + return make_double4(double(a.x), double(a.y), double(a.z), double(a.w)); +} +VECATTR double4 make_double4(float4 a) { + return make_double4(double(a.x), double(a.y), double(a.z), double(a.w)); +} //////DOUBLE4/////////////// -VECATTR double4 operator +(const double4 &a, const double4 &b){ - return make_double4(a.x + b.x, - a.y + b.y, - a.z + b.z, - a.w + b.w - ); -} -VECATTR void operator +=(double4 &a, const double4 &b){ +VECATTR double4 operator+(const double4 &a, const double4 &b) { + return make_double4(a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w); +} +VECATTR void operator+=(double4 &a, const double4 &b) { a.x += b.x; a.y += b.y; a.z += b.z; a.w += b.w; } -VECATTR double4 operator +(const double4 &a, const double &b){ - return make_double4( - a.x + b, - a.y + b, - a.z + b, - a.w + b - ); +VECATTR double4 operator+(const double4 &a, const double &b) { + return make_double4(a.x + b, a.y + b, a.z + b, a.w + b); } -VECATTR double4 operator +(const double &b, const double4 &a){ - return a+b; -} -VECATTR void operator +=(double4 &a, const double &b){ +VECATTR double4 operator+(const double &b, const double4 &a) { return a + b; } +VECATTR void operator+=(double4 &a, const double &b) { a.x += b; a.y += b; a.z += b; a.w += b; } -VECATTR double4 operator -(const double4 &a, const double4 &b){ - return make_double4( - a.x - b.x, - a.y - b.y, - a.z - b.z, - a.w - b.w - ); +VECATTR double4 operator-(const double4 &a, const double4 &b) { + return make_double4(a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w); } -VECATTR void operator -=(double4 &a, const double4 &b){ +VECATTR void operator-=(double4 &a, const double4 &b) { a.x -= b.x; a.y -= b.y; a.z -= b.z; a.w -= b.w; } -VECATTR double4 operator -(const double4 &a, const double &b){ - return make_double4( - a.x - b, - a.y - b, - a.z - b, - a.w - b - ); -} -VECATTR double4 operator -(const double &b, const double4 &a){ - return make_double4( - b - a.x, - b - a.y, - b - a.z, - b - a.w - ); -} -VECATTR void operator -=(double4 &a, const double &b){ +VECATTR double4 operator-(const double4 &a, const double &b) { + return make_double4(a.x - b, a.y - b, a.z - b, a.w - b); +} +VECATTR double4 operator-(const double &b, const double4 &a) { + return make_double4(b - a.x, b - a.y, b - a.z, b - a.w); +} +VECATTR void operator-=(double4 &a, const double &b) { a.x -= b; a.y -= b; a.z -= b; a.w -= b; } -VECATTR double4 operator *(const double4 &a, const double4 &b){ - return make_double4( - a.x * b.x, - a.y * b.y, - a.z * b.z, - a.w * b.w - ); +VECATTR double4 operator*(const double4 &a, const double4 &b) { + return make_double4(a.x * b.x, a.y * b.y, a.z * b.z, a.w * b.w); } -VECATTR void operator *=(double4 &a, const double4 &b){ +VECATTR void operator*=(double4 &a, const double4 &b) { a.x *= b.x; a.y *= b.y; a.z *= b.z; a.w *= b.w; } -VECATTR double4 operator *(const double4 &a, const double &b){ - return make_double4( - a.x * b, - a.y * b, - a.z * b, - a.w * b - ); +VECATTR double4 operator*(const double4 &a, const double &b) { + return make_double4(a.x * b, a.y * b, a.z * b, a.w * b); } -VECATTR double4 operator *(const double &b, const double4 &a){ - return a*b; -} -VECATTR void operator *=(double4 &a, const double &b){ +VECATTR double4 operator*(const double &b, const double4 &a) { return a * b; } +VECATTR void operator*=(double4 &a, const double &b) { a.x *= b; a.y *= b; a.z *= b; a.w *= b; } -VECATTR double4 operator /(const double4 &a, const double4 &b){ - return make_double4( - a.x / b.x, - a.y / b.y, - a.z / b.z, - a.w / b.w - ); +VECATTR double4 operator/(const double4 &a, const double4 &b) { + return make_double4(a.x / b.x, a.y / b.y, a.z / b.z, a.w / b.w); } -VECATTR void operator /=(double4 &a, const double4 &b){ +VECATTR void operator/=(double4 &a, const double4 &b) { a.x /= b.x; a.y /= b.y; a.z /= b.z; a.w /= b.w; } -VECATTR double4 operator /(const double4 &a, const double &b){return (1.0/b)*a;} -VECATTR double4 operator /(const double &b, const double4 &a){ - return make_double4( - b / a.x, - b / a.y, - b / a.z, - b / a.w - ); +VECATTR double4 operator/(const double4 &a, const double &b) { + return (1.0 / b) * a; } -VECATTR void operator /=(double4 &a, const double &b){ - a *= 1.0/b; +VECATTR double4 operator/(const double &b, const double4 &a) { + return make_double4(b / a.x, b / a.y, b / a.z, b / a.w); } +VECATTR void operator/=(double4 &a, const double &b) { a *= 1.0 / b; } -VECATTR double dot(double4 a, double4 b) -{ +VECATTR double dot(double4 a, double4 b) { return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w; } -VECATTR double length(double4 v) -{ - return sqrt(dot(v, v)); -} -VECATTR double4 normalize(double4 v) -{ - double invLen = 1.0/sqrt(dot(v, v)); +VECATTR double length(double4 v) { return sqrt(dot(v, v)); } +VECATTR double4 normalize(double4 v) { + double invLen = 1.0 / sqrt(dot(v, v)); return v * invLen; } -VECATTR double4 floorf(double4 v) -{ +VECATTR double4 floorf(double4 v) { return make_double4(floor(v.x), floor(v.y), floor(v.z), floor(v.w)); } /////////////////////DOUBLE3/////////////////////////////// -VECATTR int3 make_int3(double3 a){ +VECATTR int3 make_int3(double3 a) { return make_int3((int)a.x, (int)a.y, (int)a.z); } -VECATTR double3 make_double3(double a){ - return make_double3(a, a, a); -} +VECATTR double3 make_double3(double a) { return make_double3(a, a, a); } -VECATTR double3 make_double3(double2 xy, double z){ +VECATTR double3 make_double3(double2 xy, double z) { return make_double3(xy.x, xy.y, z); } -VECATTR double3 make_double3(float2 xy, double z){ +VECATTR double3 make_double3(float2 xy, double z) { return make_double3(xy.x, xy.y, z); } -VECATTR double3 make_double3(double x, double2 yz){ +VECATTR double3 make_double3(double x, double2 yz) { return make_double3(x, yz.x, yz.y); } #ifdef SINGLE_PRECISION -VECATTR double3 make_double3(double3 a){ - return a; -} +VECATTR double3 make_double3(double3 a) { return a; } #endif +VECATTR double3 make_double3(int3 a) { return make_double3(a.x, a.y, a.z); } +VECATTR double3 make_double3(float3 a) { return make_double3(a.x, a.y, a.z); } -VECATTR double3 make_double3(int3 a){ - return make_double3(a.x, a.y, a.z); -} -VECATTR double3 make_double3(float3 a){ - return make_double3(a.x, a.y, a.z); -} - -VECATTR double3 operator +(const double3 &a, const double3 &b){ - return make_double3( - a.x + b.x, - a.y + b.y, - a.z + b.z - ); +VECATTR double3 operator+(const double3 &a, const double3 &b) { + return make_double3(a.x + b.x, a.y + b.y, a.z + b.z); } -VECATTR void operator +=(double3 &a, const double3 &b){ +VECATTR void operator+=(double3 &a, const double3 &b) { a.x += b.x; a.y += b.y; a.z += b.z; } -VECATTR double3 operator +(const double3 &a, const double &b){ - return make_double3( - a.x + b, - a.y + b, - a.z + b - ); +VECATTR double3 operator+(const double3 &a, const double &b) { + return make_double3(a.x + b, a.y + b, a.z + b); } -VECATTR double3 operator +(const double &b, const double3 &a){ - return a+b; -} -VECATTR void operator +=(double3 &a, const double &b){ +VECATTR double3 operator+(const double &b, const double3 &a) { return a + b; } +VECATTR void operator+=(double3 &a, const double &b) { a.x += b; a.y += b; a.z += b; } -VECATTR double3 operator -(const double3 &a, const double3 &b){ - return make_double3( - a.x - b.x, - a.y - b.y, - a.z - b.z - ); +VECATTR double3 operator-(const double3 &a, const double3 &b) { + return make_double3(a.x - b.x, a.y - b.y, a.z - b.z); } -VECATTR void operator -=(double3 &a, const double3 &b){ +VECATTR void operator-=(double3 &a, const double3 &b) { a.x -= b.x; a.y -= b.y; a.z -= b.z; } -VECATTR double3 operator -(const double3 &a, const double &b){ - return make_double3( - a.x - b, - a.y - b, - a.z - b - ); -} -VECATTR double3 operator -(const double &b, const double3 &a){ - return make_double3( - b-a.x, - b-a.y, - b-a.z - ); -} -VECATTR void operator -=(double3 &a, const double &b){ +VECATTR double3 operator-(const double3 &a, const double &b) { + return make_double3(a.x - b, a.y - b, a.z - b); +} +VECATTR double3 operator-(const double &b, const double3 &a) { + return make_double3(b - a.x, b - a.y, b - a.z); +} +VECATTR void operator-=(double3 &a, const double &b) { a.x -= b; a.y -= b; a.z -= b; } -VECATTR double3 operator *(const double3 &a, const double3 &b){ - return make_double3( - a.x * b.x, - a.y * b.y, - a.z * b.z - ); +VECATTR double3 operator*(const double3 &a, const double3 &b) { + return make_double3(a.x * b.x, a.y * b.y, a.z * b.z); } -VECATTR void operator *=(double3 &a, const double3 &b){ +VECATTR void operator*=(double3 &a, const double3 &b) { a.x *= b.x; a.y *= b.y; a.z *= b.z; } -VECATTR double3 operator *(const double3 &a, const double &b){ - return make_double3( - a.x * b, - a.y * b, - a.z * b - ); -} -VECATTR double3 operator *(const double &b, const double3 &a){ - return a*b; +VECATTR double3 operator*(const double3 &a, const double &b) { + return make_double3(a.x * b, a.y * b, a.z * b); } -VECATTR void operator *=(double3 &a, const double &b){ +VECATTR double3 operator*(const double &b, const double3 &a) { return a * b; } +VECATTR void operator*=(double3 &a, const double &b) { a.x *= b; a.y *= b; a.z *= b; } -VECATTR double3 operator /(const double3 &a, const double3 &b){ - return make_double3( - a.x / b.x, - a.y / b.y, - a.z / b.z - ); +VECATTR double3 operator/(const double3 &a, const double3 &b) { + return make_double3(a.x / b.x, a.y / b.y, a.z / b.z); } -VECATTR void operator /=(double3 &a, const double3 &b){ +VECATTR void operator/=(double3 &a, const double3 &b) { a.x /= b.x; a.y /= b.y; a.z /= b.z; } -VECATTR double3 operator /(const double3 &a, const double &b){return (1.0/b)*a;} - -VECATTR double3 operator /(const double &b, const double3 &a){ - return make_double3( - b / a.x, - b / a.y, - b / a.z - ); +VECATTR double3 operator/(const double3 &a, const double &b) { + return (1.0 / b) * a; } -VECATTR void operator /=(double3 &a, const double &b){ - - a *= 1.0/b; +VECATTR double3 operator/(const double &b, const double3 &a) { + return make_double3(b / a.x, b / a.y, b / a.z); } +VECATTR void operator/=(double3 &a, const double &b) { a *= 1.0 / b; } -//DOUBLE2 - +// DOUBLE2 -VECATTR double2 operator -(const double2 &a, const double2 &b){ - return make_double2( - a.x - b.x, - a.y - b.y - ); +VECATTR double2 operator-(const double2 &a, const double2 &b) { + return make_double2(a.x - b.x, a.y - b.y); } -VECATTR void operator -=(double2 &a, const double2 &b){ +VECATTR void operator-=(double2 &a, const double2 &b) { a.x -= b.x; a.y -= b.y; } -VECATTR double2 operator -(const double2 &a, const double &b){ - return make_double2( - a.x - b, - a.y - b - ); +VECATTR double2 operator-(const double2 &a, const double &b) { + return make_double2(a.x - b, a.y - b); } -VECATTR double2 operator -(const double &b, const double2 &a){ - return make_double2( - b - a.x, - b - a.y - ); +VECATTR double2 operator-(const double &b, const double2 &a) { + return make_double2(b - a.x, b - a.y); +} +VECATTR void operator-=(double2 &a, const double &b) { + a.x -= b; + a.y -= b; } -VECATTR void operator -=(double2 &a, const double &b){a.x -= b; a.y -= b;} - - - -VECATTR double2 operator +(const double2 &a, const double2 &b){ - return make_double2( - a.x + b.x, - a.y + b.y - ); +VECATTR double2 operator+(const double2 &a, const double2 &b) { + return make_double2(a.x + b.x, a.y + b.y); } -VECATTR void operator +=(double2 &a, const double2 &b){ +VECATTR void operator+=(double2 &a, const double2 &b) { a.x += b.x; a.y += b.y; } -VECATTR double2 operator +(const double2 &a, const double &b){ - return make_double2( - a.x + b, - a.y + b - ); +VECATTR double2 operator+(const double2 &a, const double &b) { + return make_double2(a.x + b, a.y + b); +} +VECATTR double2 operator+(const double &b, const double2 &a) { return a + b; } +VECATTR void operator+=(double2 &a, const double &b) { + a.x += b; + a.y += b; } -VECATTR double2 operator +(const double &b, const double2 &a){ return a+b;} -VECATTR void operator +=(double2 &a, const double &b){a.x += b; a.y += b;} - -VECATTR double2 operator *(const double2 &a, const double2 &b){ +VECATTR double2 operator*(const double2 &a, const double2 &b) { return make_double2(a.x * b.x, a.y * b.y); } -VECATTR void operator *=(double2 &a, const double2 &b){ +VECATTR void operator*=(double2 &a, const double2 &b) { a.x *= b.x; a.y *= b.y; } -VECATTR double2 operator *(const double2 &a, const double &b){ +VECATTR double2 operator*(const double2 &a, const double &b) { return make_double2(a.x * b, a.y * b); } -VECATTR double2 operator *(const double &b, const double2 &a){ - return a*b; -} -VECATTR void operator *=(double2 &a, const double &b){ +VECATTR double2 operator*(const double &b, const double2 &a) { return a * b; } +VECATTR void operator*=(double2 &a, const double &b) { a.x *= b; a.y *= b; } - - -VECATTR double2 operator /(const double2 &a, const double2 &b){ +VECATTR double2 operator/(const double2 &a, const double2 &b) { return make_double2(a.x / b.x, a.y / b.y); } -VECATTR void operator /=(double2 &a, const double2 &b){ +VECATTR void operator/=(double2 &a, const double2 &b) { a.x /= b.x; a.y /= b.y; } -VECATTR double2 operator /(const double2 &a, const double &b){ +VECATTR double2 operator/(const double2 &a, const double &b) { return make_double2(a.x / b, a.y / b); } -VECATTR double2 operator /(const double &b, const double2 &a){ - return make_double2(b/a.x, b/a.y); +VECATTR double2 operator/(const double &b, const double2 &a) { + return make_double2(b / a.x, b / a.y); } -VECATTR void operator /=(double2 &a, const double &b){ +VECATTR void operator/=(double2 &a, const double &b) { a.x /= b; a.y /= b; } - - //////////////////////////// -VECATTR double3 floorf(double3 v){return make_double3(floor(v.x), floor(v.y), floor(v.z));} - - - -VECATTR double dot(const double3 &a, const double3 &b){return a.x * b.x + a.y * b.y + a.z * b.z;} -VECATTR float dot(const float3 &a, const float3 &b){return a.x * b.x + a.y * b.y + a.z * b.z;} -VECATTR int dot(const int3 &a, const int3 &b){return a.x * b.x + a.y * b.y + a.z * b.z;} -VECATTR int dot(const int2 &a, const int2 &b){return a.x * b.x + a.y * b.y;} +VECATTR double3 floorf(double3 v) { + return make_double3(floor(v.x), floor(v.y), floor(v.z)); +} +VECATTR double dot(const double3 &a, const double3 &b) { + return a.x * b.x + a.y * b.y + a.z * b.z; +} +VECATTR float dot(const float3 &a, const float3 &b) { + return a.x * b.x + a.y * b.y + a.z * b.z; +} +VECATTR int dot(const int3 &a, const int3 &b) { + return a.x * b.x + a.y * b.y + a.z * b.z; +} +VECATTR int dot(const int2 &a, const int2 &b) { return a.x * b.x + a.y * b.y; } -VECATTR double length(double3 v){return sqrt(dot(v, v));} -VECATTR double3 normalize(double3 v) -{ - double invLen = 1.0/sqrt(dot(v, v)); +VECATTR double length(double3 v) { return sqrt(dot(v, v)); } +VECATTR double3 normalize(double3 v) { + double invLen = 1.0 / sqrt(dot(v, v)); return v * invLen; } -VECATTR double3 cross(double3 a, double3 b){ - return make_double3(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x); +VECATTR double3 cross(double3 a, double3 b) { + return make_double3(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, + a.x * b.y - a.y * b.x); } -VECATTR float3 cross(float3 a, float3 b){ - return make_float3(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x); +VECATTR float3 cross(float3 a, float3 b) { + return make_float3(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, + a.x * b.y - a.y * b.x); } - -VECATTR float3 sqrt(const float3 &a){ return {sqrt(a.x), sqrt(a.y), sqrt(a.z)};} -VECATTR double3 sqrt(const double3 &a){ return {sqrt(a.x), sqrt(a.y), sqrt(a.z)};} - +VECATTR float3 sqrt(const float3 &a) { + return {sqrt(a.x), sqrt(a.y), sqrt(a.z)}; +} +VECATTR double3 sqrt(const double3 &a) { + return {sqrt(a.x), sqrt(a.y), sqrt(a.z)}; +} ////////////////////////////////////////////////////////// - /****************************************************************************************/ - - ///////////INT2///////////////// -VECATTR int2 make_int2(int3 a){return make_int2(a.x, a.y);} +VECATTR int2 make_int2(int3 a) { return make_int2(a.x, a.y); } ///////////INT3///////////////// -VECATTR int3 make_int3(int a){return make_int3(a,a,a);} -VECATTR int3 make_int3(int2 a, int b){return make_int3(a.x,a.y,b);} +VECATTR int3 make_int3(int a) { return make_int3(a, a, a); } +VECATTR int3 make_int3(int2 a, int b) { return make_int3(a.x, a.y, b); } -VECATTR int3 operator /(int3 a, int3 b){ - return make_int3( a.x/b.x, a.y/b.y, a.z/b.z); +VECATTR int3 operator/(int3 a, int3 b) { + return make_int3(a.x / b.x, a.y / b.y, a.z / b.z); } -VECATTR int3 operator /(int3 a, int b){ - return make_int3( a.x/b, a.y/b, a.z/b); +VECATTR int3 operator/(int3 a, int b) { + return make_int3(a.x / b, a.y / b, a.z / b); } -VECATTR int3 operator /(int a, int3 b){ - return make_int3( a/b.x, a/b.y, a/b.z); +VECATTR int3 operator/(int a, int3 b) { + return make_int3(a / b.x, a / b.y, a / b.z); } -VECATTR int3 operator +(const int3 &a, const int3 &b){ +VECATTR int3 operator+(const int3 &a, const int3 &b) { return make_int3(a.x + b.x, a.y + b.y, a.z + b.z); } -VECATTR int3 operator +(const int3 &a, const int &b){ +VECATTR int3 operator+(const int3 &a, const int &b) { return make_int3(a.x + b, a.y + b, a.z + b); } -VECATTR int3 operator -(const int3 &a, const int3 &b){ +VECATTR int3 operator-(const int3 &a, const int3 &b) { return make_int3(a.x - b.x, a.y - b.y, a.z - b.z); } -VECATTR void operator -=(int3 &a, const int3 &b){ - a.x -= b.x; a.y -= b.y; a.z -= b.z; +VECATTR void operator-=(int3 &a, const int3 &b) { + a.x -= b.x; + a.y -= b.y; + a.z -= b.z; } -VECATTR void operator +=(int3 &a, const int3 &b){ - a.x += b.x; a.y += b.y; a.z += b.z; +VECATTR void operator+=(int3 &a, const int3 &b) { + a.x += b.x; + a.y += b.y; + a.z += b.z; } -VECATTR void operator /=(int3 &a, const int &b){ - a.x /= b; a.y /= b; a.z /= b; +VECATTR void operator/=(int3 &a, const int &b) { + a.x /= b; + a.y /= b; + a.z /= b; } -VECATTR int3 operator -(const int3 &a, const int &b){ +VECATTR int3 operator-(const int3 &a, const int &b) { return make_int3(a.x - b, a.y - b, a.z - b); } -VECATTR int3 operator -(const int &b, const int3 &a){ +VECATTR int3 operator-(const int &b, const int3 &a) { return make_int3(b - a.x, b - a.y, b - a.z); } - -VECATTR int3 operator *(const int3 &a, const int3 &b){return make_int3(a.x*b.x, a.y*b.y, a.z*b.z);} -VECATTR void operator *=(int3 &a, const int3 &b) {a = a*b;} -VECATTR int3 operator *(const int3 &a, const int &b){return make_int3(a.x*b, a.y*b, a.z*b);} -VECATTR int3 operator *(const int &b, const int3 &a){return a*b;} +VECATTR int3 operator*(const int3 &a, const int3 &b) { + return make_int3(a.x * b.x, a.y * b.y, a.z * b.z); +} +VECATTR void operator*=(int3 &a, const int3 &b) { a = a * b; } +VECATTR int3 operator*(const int3 &a, const int &b) { + return make_int3(a.x * b, a.y * b, a.z * b); +} +VECATTR int3 operator*(const int &b, const int3 &a) { return a * b; } #endif diff --git a/solvers/NBody/python_wrapper.cu b/solvers/NBody/python_wrapper.cu index d955553b..ae4b9af3 100644 --- a/solvers/NBody/python_wrapper.cu +++ b/solvers/NBody/python_wrapper.cu @@ -39,12 +39,12 @@ MOBILITY_PYTHONIFY_WITH_EXTRA_CODE( NBody, solver.def( "setParameters", - [](NBody &myself, std::string algo, int NBatch, int NperBatch, std::optional wallHeight) - { - myself.setParametersNBody( - {nbody_rpy::string2NBodyAlgorithm(algo), NBatch, NperBatch, wallHeight}); + [](NBody &myself, std::string algo, int NBatch, int NperBatch, + std::optional wallHeight) { + myself.setParametersNBody({nbody_rpy::string2NBodyAlgorithm(algo), + NBatch, NperBatch, wallHeight}); }, - docstringSetParameters, - "algorithm"_a = "advise", "Nbatch"_a = -1, "NperBatch"_a = -1, "wallHeight"_a = std::nullopt); + docstringSetParameters, "algorithm"_a = "advise", "Nbatch"_a = -1, + "NperBatch"_a = -1, "wallHeight"_a = std::nullopt); , "This module computes the RPY mobility using an N^2 algorithm in the " "GPU. Different hydrodynamic kernels can be chosen."); diff --git a/solvers/NBody/test.py b/solvers/NBody/test.py deleted file mode 100644 index 90126253..00000000 --- a/solvers/NBody/test.py +++ /dev/null @@ -1,20 +0,0 @@ -from NBody import * -import numpy as np - -#Constructor requires periodicity in each dimension -nb = NBody("open", "open", "open") - -numberParticles = 10 -precision = np.float32 if nb.precision=="float" else np.float64 -pos = np.linspace(-10, 10, 3*numberParticles).astype(precision) -force = np.linspace(-1, 1, 3*numberParticles).astype(precision) -result = np.zeros(3*numberParticles).astype(precision) - -par = Parameters(temperature = 0.0, viscosity = 1.0, hydrodynamicRadius = 1.0) -nb.initialize(par); -#Batched parameters to -1 turns the functionality off (default) -nb.setParametersNBody(algorithm="advise", Nbatch=-1, NperBatch=-1); -nb.setPositions(pos, numberParticles) -nb.Mdot(forces = force, result = result) -nb.clean() -print(result) diff --git a/solvers/PSE/example.py b/solvers/PSE/example.py deleted file mode 100644 index 7127af73..00000000 --- a/solvers/PSE/example.py +++ /dev/null @@ -1,19 +0,0 @@ -#Raul P. Pelaez 2021. PSE Mobility python example. -from PSE import * -#help(PSE) -import numpy as np - -pse = PSE() - -numberParticles = 10 -precision = np.float32 if pse.precision=="float" else np.float64 -pos = np.linspace(-10, 10, 3*numberParticles).astype(precision) -force = np.linspace(-1, 1, 3*numberParticles).astype(precision) -result = np.zeros(3*numberParticles).astype(precision) - -par = Parameters(temperature = 1.0, viscosity = 1.0, hydrodynamicRadius = 1.0, boxSize=BoxSize(128,128,128), shearStrain = 0) -pse.initialize(par); -pse.setPositions(pos, numberParticles) -pse.Mdot(forces = force, result = result) -pse.clean() -print(result) diff --git a/solvers/PSE/extra/uammd_interface.h b/solvers/PSE/extra/uammd_interface.h index ea3d73b2..4e0031a5 100644 --- a/solvers/PSE/extra/uammd_interface.h +++ b/solvers/PSE/extra/uammd_interface.h @@ -3,48 +3,50 @@ */ #ifndef UAMMD_INTERFACE_H #define UAMMD_INTERFACE_H -#include -#include -namespace uammd_pse{ - //This is in order not to use any UAMMD related includes here. - //Instead of using uammd::real I have to re define real here. +#include +#include +namespace uammd_pse { +// This is in order not to use any UAMMD related includes here. +// Instead of using uammd::real I have to re define real here. #ifndef DOUBLE_PRECISION - using real = float; +using real = float; #else - using real = double; +using real = double; #endif - //This function returns either 'single' or 'double' according to the UAMMD's compiled precision. - namespace uammd_wrapper{ - std::string getPrecision(); - } +// This function returns either 'single' or 'double' according to the UAMMD's +// compiled precision. +namespace uammd_wrapper { +std::string getPrecision(); +} - struct PyParameters{ - real viscosity; - real hydrodynamicRadius; - real Lx, Ly, Lz; - real tolerance; - real psi; - real shearStrain; - }; +struct PyParameters { + real viscosity; + real hydrodynamicRadius; + real Lx, Ly, Lz; + real tolerance; + real psi; + real shearStrain; +}; - class UAMMD_PSE; - class UAMMD_PSE_Glue{ - std::shared_ptr pse; - public: +class UAMMD_PSE; +class UAMMD_PSE_Glue { + std::shared_ptr pse; - UAMMD_PSE_Glue(PyParameters pypar, int numberParticles); +public: + UAMMD_PSE_Glue(PyParameters pypar, int numberParticles); - void MdotNearField(const real* h_pos, const real* h_F, real* h_MF); + void MdotNearField(const real *h_pos, const real *h_F, real *h_MF); - void MdotFarField(const real* h_pos, const real* h_F, real* h_MF); + void MdotFarField(const real *h_pos, const real *h_F, real *h_MF); - void computeHydrodynamicDisplacements(const real* h_pos, const real* h_F, real* h_MF, - real temperature, real prefactor); + void computeHydrodynamicDisplacements(const real *h_pos, const real *h_F, + real *h_MF, real temperature, + real prefactor); - void setShearStrain(real newStrain); + void setShearStrain(real newStrain); - void clean(); - }; -} + void clean(); +}; +} // namespace uammd_pse #endif diff --git a/solvers/PSE/extra/uammd_wrapper.cu b/solvers/PSE/extra/uammd_wrapper.cu index 4f34d478..bc00c05c 100644 --- a/solvers/PSE/extra/uammd_wrapper.cu +++ b/solvers/PSE/extra/uammd_wrapper.cu @@ -9,7 +9,7 @@ #include namespace uammd_pse { - using namespace uammd; +using namespace uammd; class Gaussian { IBM_kernels::Gaussian kern; diff --git a/solvers/PSE/python_wrapper.cu b/solvers/PSE/python_wrapper.cu index 51bf190a..8b02a740 100644 --- a/solvers/PSE/python_wrapper.cu +++ b/solvers/PSE/python_wrapper.cu @@ -28,8 +28,6 @@ MOBILITY_PYTHONIFY_WITH_EXTRA_CODE( [](PSE &self, real psi, real Lx, real Ly, real Lz, real shearStrain) { self.setParametersPSE({psi, Lx, Ly, Lz, shearStrain}); }, - docstringSetParameters, - "psi"_a, "Lx"_a, - "Ly"_a, "Lz"_a, + docstringSetParameters, "psi"_a, "Lx"_a, "Ly"_a, "Lz"_a, "shearStrain"_a); , docstring); diff --git a/solvers/PSE/test.py b/solvers/PSE/test.py deleted file mode 100644 index e2b6aaae..00000000 --- a/solvers/PSE/test.py +++ /dev/null @@ -1,22 +0,0 @@ -#Raul P. Pelaez 2021-2022. PSE Mobility test. -from PSE import * -import numpy as np - -pse = PSE("periodic", "periodic", "periodic") - -numberParticles = 10 -precision = np.float32 if pse.precision=="float" else np.float64 -pos = np.linspace(-10, 10, 3*numberParticles).astype(precision) -force = np.linspace(-1, 1, 3*numberParticles).astype(precision) -result = np.zeros(3*numberParticles).astype(precision) - - -pse.setParametersPSE(psi=1.0, Lx=32, Ly=32, Lz=32, shearStrain=0) -pse.initialize(temperature=1.0, viscosity=1.0, hydrodynamicRadius=1, numberParticles=numberParticles); -#If the second call to setParametersPSE+initialize only differs in the shear strain the module is not reinitialized. Only the shear strain is updated -#pse.setParametersPSE(psi=1.0, Lx=32, Ly=32, Lz=32, shearStrain=1) -#pse.initialize(temperature=1.0, viscosity=1.0, hydrodynamicRadius=1, numberParticles=numberParticles); -pse.setPositions(pos) -pse.Mdot(forces = force, result = result) -pse.clean() -print(result) diff --git a/solvers/SelfMobility/python_wrapper.cu b/solvers/SelfMobility/python_wrapper.cu index f6b26817..6b992dce 100644 --- a/solvers/SelfMobility/python_wrapper.cu +++ b/solvers/SelfMobility/python_wrapper.cu @@ -1,14 +1,18 @@ /*Raul P. Pelaez 2021-2022. Python wrapper for the SelfMobility module -*/ -#include"mobility.h" + */ +#include "mobility.h" #include -MOBILITY_PYTHONIFY_WITH_EXTRA_CODE(SelfMobility, - solver.def("setParameters", - [](SelfMobility &self, real parameter){ - self.setParametersSelfMobility(parameter); - }, "parameter"_a, "Some example parameter");, - "This module ignores hydrodynamic interactions, AKA the mobility matrix is simply (1/(6*pi*eta*a))*I"); +MOBILITY_PYTHONIFY_WITH_EXTRA_CODE( + SelfMobility, solver.def( + "setParameters", + [](SelfMobility &self, real parameter) { + self.setParametersSelfMobility(parameter); + }, + "parameter"_a, "Some example parameter"); + , "This module ignores hydrodynamic interactions, AKA the mobility matrix " + "is simply (1/(6*pi*eta*a))*I"); -//A module that does not use setParametersSolver can pythonify like this instead -//MOBILITY_PYTHONIFY(SelfMobility, "This module ignores hydrodynamic interactions, AKA the mobility matrix is simply (1/(6*pi*eta*a))*I"); +// A module that does not use setParametersSolver can pythonify like this +// instead MOBILITY_PYTHONIFY(SelfMobility, "This module ignores hydrodynamic +// interactions, AKA the mobility matrix is simply (1/(6*pi*eta*a))*I"); diff --git a/tests/test_array_origins.py b/tests/test_array_origins.py index 6227c440..240029ab 100644 --- a/tests/test_array_origins.py +++ b/tests/test_array_origins.py @@ -64,7 +64,9 @@ def setup_inputs(framework, use_torques, numberParticles, precision): @pytest.mark.parametrize(("Solver", "periodicity"), solver_configs_all) @pytest.mark.parametrize("framework", ["numpy", "torch", "cupy", "jax", "tensorflow"]) @pytest.mark.parametrize("use_torques", [False, True]) -@pytest.mark.parametrize("method", ["Mdot", "thermalDrift", "sqrtMdotW", "hydrodynamicVelocities"]) +@pytest.mark.parametrize( + "method", ["Mdot", "thermalDrift", "sqrtMdotW", "hydrodynamicVelocities"] +) def test_framework(Solver, periodicity, framework, use_torques, method): if use_torques and Solver.__name__ in "PSE": pytest.skip("PSE does not support torques") diff --git a/tests/test_interface.py b/tests/test_interface.py index 7aeb2d1c..ab75a654 100644 --- a/tests/test_interface.py +++ b/tests/test_interface.py @@ -1,4 +1,5 @@ import pytest + # from libMobility import * import numpy as np from utils import ( @@ -87,6 +88,7 @@ def check_results(mf, mt, ivec, Solver): mf, mt = solver.Mdot(torques=torques) check_results(mf, mt, torques, Solver) + @pytest.mark.parametrize(("Solver", "periodicity"), solver_configs_torques) def test_forces_with_includeAngular_returns_both(Solver, periodicity): @@ -152,6 +154,7 @@ def check_results(mf, mt, fvec, tvec): u, w = solver.Mdot(forces, torques) check_results(u, w, forces, torques) + @pytest.mark.parametrize(("Solver", "periodicity"), solver_configs_all) def test_returns_sqrtM(Solver, periodicity): numberParticles = 1 diff --git a/tests/test_thermal_drift.py b/tests/test_thermal_drift.py index 23c614e0..e0d266db 100644 --- a/tests/test_thermal_drift.py +++ b/tests/test_thermal_drift.py @@ -31,8 +31,10 @@ def thermal_drift_func(): _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 + _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 solver.setPositions(positions) @@ -74,7 +76,7 @@ def test_thermal_drift_does_not_change_positions(Solver, periodicity): @pytest.mark.parametrize("numberParticles", [1, 2, 3, 10]) @pytest.mark.parametrize("includeAngular", [False, True]) def test_thermal_drift_is_zero( - Solver, periodicity, hydrodynamicRadius, numberParticles, includeAngular + Solver, periodicity, hydrodynamicRadius, numberParticles, includeAngular ): if not np.all(np.array(periodicity) == "open") and not np.all( np.array(periodicity) == "periodic" @@ -82,7 +84,7 @@ def test_thermal_drift_is_zero( pytest.skip( "Only periodic and open boundary conditions have zero thermal drift" ) - if(Solver.__name__ == "PSE" and 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) @@ -112,15 +114,14 @@ def test_thermal_drift_is_zero( ), 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]) @pytest.mark.parametrize("includeAngular", [False, True]) def test_thermal_drift_returns_different_numbers( - Solver, periodicity, hydrodynamicRadius, numberParticles, includeAngular + Solver, periodicity, hydrodynamicRadius, numberParticles, includeAngular ): - if(Solver.__name__ == "PSE" and includeAngular): + if Solver.__name__ == "PSE" and includeAngular: pytest.skip("PSE does not support torques") temperature = 1.2 @@ -138,10 +139,10 @@ def test_thermal_drift_returns_different_numbers( generate_positions_in_box(parameters, numberParticles).astype(precision) * 0.8 ) solver.setPositions(positions) - rfd1,_ = solver.thermalDrift() + rfd1, _ = solver.thermalDrift() if np.all(rfd1 == 0): pytest.skip("RFD is zero, skipping test") - rfd2,_ = solver.thermalDrift() + rfd2, _ = solver.thermalDrift() assert np.any( np.abs(rfd1 - rfd2) > 1e-5 ), f"RFD is not different: {np.max(np.abs(rfd1 - rfd2))}" @@ -152,9 +153,9 @@ def test_thermal_drift_returns_different_numbers( @pytest.mark.parametrize("numberParticles", [1, 10]) @pytest.mark.parametrize("includeAngular", [False, True]) def test_thermal_drift_matches_rfd( - Solver, periodicity, hydrodynamicRadius, numberParticles, includeAngular + Solver, periodicity, hydrodynamicRadius, numberParticles, includeAngular ): - if(Solver.__name__ == "PSE" and includeAngular): + if Solver.__name__ == "PSE" and includeAngular: pytest.skip("PSE does not support torques") temperature = 1.2 precision = np.float32 if Solver.precision == "float" else np.float64