From c757eac283d297e615ff39174cb51f9ba0ceb062 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 16:38:33 +0200 Subject: [PATCH 01/24] style: Add a clang-format --- .clang-format | 274 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 274 insertions(+) create mode 100644 .clang-format 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 +... + From d137322c311cee0767eb9e253666010cf127ea32 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 16:53:55 +0200 Subject: [PATCH 02/24] Remove old file --- solvers/NBody/example.py | 3 --- 1 file changed, 3 deletions(-) delete mode 100644 solvers/NBody/example.py 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) From 6b279b5b317b50c24419c78b640c012ffe6d8b8f Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 16:55:18 +0200 Subject: [PATCH 03/24] Remove more old files --- solvers/NBody/test.py | 20 -------------------- solvers/PSE/example.py | 19 ------------------- solvers/PSE/test.py | 22 ---------------------- 3 files changed, 61 deletions(-) delete mode 100644 solvers/NBody/test.py delete mode 100644 solvers/PSE/example.py delete mode 100644 solvers/PSE/test.py 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/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) From 4321a9a9683bcc87e4a14fe2cc51c4af1fceb0ea Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 17:00:11 +0200 Subject: [PATCH 04/24] style: add a linter check CI --- .github/workflows/format.yml | 78 ++++++++++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100644 .github/workflows/format.yml diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml new file mode 100644 index 00000000..9f854652 --- /dev/null +++ b/.github/workflows/format.yml @@ -0,0 +1,78 @@ +name: Clang and Black Format Lint Check + +on: + pull_request: + types: [opened, synchronize, reopened] + +jobs: + format-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 dependencies + run: | + sudo apt-get update && sudo apt-get install -y clang-format-19 + pip install black + + - name: Check C++ formatting with clang-format + id: clang_check + run: | + FILE_LIST="unformatted_cpp.txt" + find src ! -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 'cpp_files<> "$GITHUB_OUTPUT" + else + echo 'cpp_files=' >> "$GITHUB_OUTPUT" + fi + + - name: Check Python formatting with Black + id: black_check + run: | + OUTPUT=$(black --check --diff python || true) + echo "$OUTPUT" + if echo "$OUTPUT" | grep -q "would reformat"; then + { + echo 'py_output<> "$GITHUB_OUTPUT" + else + echo 'py_output=' >> "$GITHUB_OUTPUT" + fi + + - name: Comment on PR if any formatting issues found + if: steps.clang_check.outputs.cpp_files != '' || steps.black_check.outputs.py_output != '' + uses: peter-evans/create-or-update-comment@v4 + with: + issue-number: ${{ github.event.pull_request.number }} + body: | + ⚠️ **Formatting check failed** + + ${{ steps.clang_check.outputs.cpp_files != '' && format('### C++ files not properly formatted:\n```\n{0}\n```', steps.clang_check.outputs.cpp_files) || '' }} + + ${{ steps.black_check.outputs.py_output != '' && format('### Python files not properly formatted:\n```\n{0}\n```', steps.black_check.outputs.py_output) || '' }} + + Please format these files using: + - `clang-format -i` for C++ files + - `black python/` for Python files + + - name: Fail job if any formatting issues found + if: steps.clang_check.outputs.cpp_files != '' || steps.black_check.outputs.py_output != '' + run: | + echo "Formatting issues detected. Failing the job." + exit 1 From ebf402a589740f9355836538aa0483782aef039f Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 17:03:44 +0200 Subject: [PATCH 05/24] doc: fix incorrect info in new solver instructions --- docs/source/new-solver.rst | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) 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. From ded9d90bfb96f2879522441bf39d40bf00457cf7 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 17:07:33 +0200 Subject: [PATCH 06/24] style: fix path --- .github/workflows/format.yml | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index 9f854652..712566f8 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -26,8 +26,7 @@ jobs: id: clang_check run: | FILE_LIST="unformatted_cpp.txt" - find src ! -path 'src/third_party/*' \ - -regex '.*\.\(cpp\|hpp\|h\|c\|cc\|cu\|cuh\)' \ + find include solvers -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 @@ -43,7 +42,7 @@ jobs: - name: Check Python formatting with Black id: black_check run: | - OUTPUT=$(black --check --diff python || true) + OUTPUT=$(black --check --diff . || true) echo "$OUTPUT" if echo "$OUTPUT" | grep -q "would reformat"; then { From caf53f3696e447dfee038cfda54f14b9235d1936 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 17:38:54 +0200 Subject: [PATCH 07/24] style: try another workflow --- .github/workflows/format.yml | 104 +++++++++++++++++------------------ 1 file changed, 50 insertions(+), 54 deletions(-) diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index 712566f8..ef40ba96 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -1,77 +1,73 @@ -name: Clang and Black Format Lint Check +name: Code Formatting Check -on: - pull_request: - types: [opened, synchronize, reopened] +on: [push, pull_request] jobs: - format-check: + lint-format: runs-on: ubuntu-latest steps: - name: Checkout code uses: actions/checkout@v4 + # -------------------- + # C++ Formatting Check + # -------------------- + - name: Run C++ format and lint check + id: cpp_lint + uses: cpp-linter/cpp-linter-action@v2 + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + with: + style: 'file' + tidy-checks: '' + thread-comments: ${{ github.event_name == 'pull_request' && 'update' }} + + - name: Record C++ status + run: echo "cpp_failed=${{ steps.cpp_lint.outputs.checks-failed }}" >> $GITHUB_ENV + + # ---------------------- + # Python Formatting Check + # ---------------------- - name: Set up Python uses: actions/setup-python@v5 with: python-version: '3.11' - - name: Install dependencies - run: | - sudo apt-get update && sudo apt-get install -y clang-format-19 - pip install black - - - name: Check C++ formatting with clang-format - id: clang_check - run: | - FILE_LIST="unformatted_cpp.txt" - find include solvers -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 'cpp_files<> "$GITHUB_OUTPUT" - else - echo 'cpp_files=' >> "$GITHUB_OUTPUT" - fi + - name: Install Black + run: pip install black - - name: Check Python formatting with Black - id: black_check + - name: Run Black and capture output + id: py_black + continue-on-error: true run: | - OUTPUT=$(black --check --diff . || true) + echo "Running black --check --diff..." + OUTPUT=$(black --check --diff src tests 2>&1 || true) echo "$OUTPUT" + echo "$OUTPUT" > black_diff.txt + if echo "$OUTPUT" | grep -q "would reformat"; then - { - echo 'py_output<> "$GITHUB_OUTPUT" + echo "py_failed=true" >> $GITHUB_ENV else - echo 'py_output=' >> "$GITHUB_OUTPUT" + echo "py_failed=false" >> $GITHUB_ENV fi - - name: Comment on PR if any formatting issues found - if: steps.clang_check.outputs.cpp_files != '' || steps.black_check.outputs.py_output != '' - uses: peter-evans/create-or-update-comment@v4 - with: - issue-number: ${{ github.event.pull_request.number }} - body: | - ⚠️ **Formatting check failed** - - ${{ steps.clang_check.outputs.cpp_files != '' && format('### C++ files not properly formatted:\n```\n{0}\n```', steps.clang_check.outputs.cpp_files) || '' }} - - ${{ steps.black_check.outputs.py_output != '' && format('### Python files not properly formatted:\n```\n{0}\n```', steps.black_check.outputs.py_output) || '' }} - - Please format these files using: - - `clang-format -i` for C++ files - - `black python/` for Python files + - name: Print unformatted Python files + if: env.py_failed == 'true' + run: | + echo "❌ Python formatting issues detected:" + grep '^would reformat' black_diff.txt | cut -d' ' -f3 | sort -u - - name: Fail job if any formatting issues found - if: steps.clang_check.outputs.cpp_files != '' || steps.black_check.outputs.py_output != '' + # ---------------------- + # Final Step: Global Fail + # ---------------------- + - name: Fail if either check failed run: | - echo "Formatting issues detected. Failing the job." - exit 1 + echo "C++ checks failed: $cpp_failed" + echo "Python checks failed: $py_failed" + if [[ "$cpp_failed" != "0" || "$py_failed" == "true" ]]; then + echo "❌ One or more format checks failed." + exit 1 + else + echo "✅ All format checks passed." + fi From 24878a8827fc37f134448c9a1b1c548c72b788fa Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 17:40:30 +0200 Subject: [PATCH 08/24] style: try another workflow --- .github/workflows/format.yml | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index ef40ba96..bd5d1a00 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -10,9 +10,6 @@ jobs: - name: Checkout code uses: actions/checkout@v4 - # -------------------- - # C++ Formatting Check - # -------------------- - name: Run C++ format and lint check id: cpp_lint uses: cpp-linter/cpp-linter-action@v2 @@ -20,7 +17,7 @@ jobs: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} with: style: 'file' - tidy-checks: '' + tidy-checks: '-*' thread-comments: ${{ github.event_name == 'pull_request' && 'update' }} - name: Record C++ status @@ -42,7 +39,7 @@ jobs: continue-on-error: true run: | echo "Running black --check --diff..." - OUTPUT=$(black --check --diff src tests 2>&1 || true) + OUTPUT=$(black --check --diff . 2>&1 || true) echo "$OUTPUT" echo "$OUTPUT" > black_diff.txt From 8a4cce5784e916a79e80fcab5b64b1057b4ce02b Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 17:45:18 +0200 Subject: [PATCH 09/24] style: try another workflow --- .github/workflows/format.yml | 84 ++++++++++-------------------------- 1 file changed, 22 insertions(+), 62 deletions(-) diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index bd5d1a00..40d0500c 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -1,70 +1,30 @@ -name: Code Formatting Check - -on: [push, pull_request] +name: cpp-linter +on: + pull_request: + branches: [main, master, develop] + paths: ['**.c', '**.cu', '**.cuh', '**.cpp', '**.h', '**.hpp', '**.cxx', '**.hxx', '**.cc', '**.hh', '**CMakeLists.txt', 'meson.build', '**.cmake'] + push: + branches: [main, master, develop] + paths: ['**.c', '**.cu', '**.cuh', '**.cpp', '**.h', '**.hpp', '**.cxx', '**.hxx', '**.cc', '**.hh', '**CMakeLists.txt', 'meson.build', '**.cmake'] jobs: - lint-format: + cpp-linter: runs-on: ubuntu-latest - + permissions: + pull-requests: write + steps: - - name: Checkout code - uses: actions/checkout@v4 - - - name: Run C++ format and lint check - id: cpp_lint - uses: cpp-linter/cpp-linter-action@v2 + - uses: actions/checkout@v4 + - uses: cpp-linter/cpp-linter-action@v2 + id: linter env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} with: - style: 'file' - tidy-checks: '-*' + style: 'file' # Use .clang-format config file. + tidy-checks: '-*' # disable clang-tidy checks. + # only 'update' a single comment in a pull request's thread. thread-comments: ${{ github.event_name == 'pull_request' && 'update' }} - - - name: Record C++ status - run: echo "cpp_failed=${{ steps.cpp_lint.outputs.checks-failed }}" >> $GITHUB_ENV - - # ---------------------- - # Python Formatting Check - # ---------------------- - - 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 capture output - id: py_black - continue-on-error: true - run: | - echo "Running black --check --diff..." - OUTPUT=$(black --check --diff . 2>&1 || true) - echo "$OUTPUT" - echo "$OUTPUT" > black_diff.txt - - if echo "$OUTPUT" | grep -q "would reformat"; then - echo "py_failed=true" >> $GITHUB_ENV - else - echo "py_failed=false" >> $GITHUB_ENV - fi - - - name: Print unformatted Python files - if: env.py_failed == 'true' - run: | - echo "❌ Python formatting issues detected:" - grep '^would reformat' black_diff.txt | cut -d' ' -f3 | sort -u - - # ---------------------- - # Final Step: Global Fail - # ---------------------- - - name: Fail if either check failed - run: | - echo "C++ checks failed: $cpp_failed" - echo "Python checks failed: $py_failed" - if [[ "$cpp_failed" != "0" || "$py_failed" == "true" ]]; then - echo "❌ One or more format checks failed." - exit 1 - else - echo "✅ All format checks passed." - fi + + - name: Fail fast?! + if: steps.linter.outputs.clang-format-checks-failed > 0 + run: exit 1 From 58c1cdc3761ae3722048283c13804301e9c87e8f Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 17:47:34 +0200 Subject: [PATCH 10/24] style: try another workflow --- .github/workflows/format.yml | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index 40d0500c..3fa3252c 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -1,11 +1,7 @@ name: cpp-linter on: pull_request: - branches: [main, master, develop] - paths: ['**.c', '**.cu', '**.cuh', '**.cpp', '**.h', '**.hpp', '**.cxx', '**.hxx', '**.cc', '**.hh', '**CMakeLists.txt', 'meson.build', '**.cmake'] - push: - branches: [main, master, develop] - paths: ['**.c', '**.cu', '**.cuh', '**.cpp', '**.h', '**.hpp', '**.cxx', '**.hxx', '**.cc', '**.hh', '**CMakeLists.txt', 'meson.build', '**.cmake'] + branches: [ main ] jobs: cpp-linter: From 5be0d509bff657388a243522a17784fcf8ea94a3 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 17:49:41 +0200 Subject: [PATCH 11/24] style: try another workflow --- .github/workflows/format.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index 3fa3252c..1648f3d3 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -17,7 +17,8 @@ jobs: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} with: style: 'file' # Use .clang-format config file. - tidy-checks: '-*' # disable clang-tidy checks. + tidy-checks: '-*' # disable clang-tidy checks. + extensions: 'cu,cuh,c,h,C,H,cpp,hpp,cc,hh,c++,h++,cxx,hxx' # only 'update' a single comment in a pull request's thread. thread-comments: ${{ github.event_name == 'pull_request' && 'update' }} From f4e7046ae7c3d9e77f1c9896d814a01409729fc7 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 17:51:10 +0200 Subject: [PATCH 12/24] style: try another workflow --- .github/workflows/format.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index 1648f3d3..7cac0382 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -19,6 +19,7 @@ jobs: style: 'file' # Use .clang-format config file. tidy-checks: '-*' # disable clang-tidy checks. extensions: 'cu,cuh,c,h,C,H,cpp,hpp,cc,hh,c++,h++,cxx,hxx' + files-changed-only: false # only 'update' a single comment in a pull request's thread. thread-comments: ${{ github.event_name == 'pull_request' && 'update' }} From 900f046b8dddccd73eb48557114d3d20c7dbfba0 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 17:57:08 +0200 Subject: [PATCH 13/24] style: try another workflow --- .github/workflows/format.yml | 70 +++++++++++++++++++++++++----------- 1 file changed, 49 insertions(+), 21 deletions(-) diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index 7cac0382..0ce51992 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -1,28 +1,56 @@ -name: cpp-linter +name: Clang-Format Lint Check + on: pull_request: - branches: [ main ] + types: [opened, synchronize, reopened] jobs: - cpp-linter: + clang-format-check: runs-on: ubuntu-latest - permissions: - pull-requests: write - + steps: - - uses: actions/checkout@v4 - - uses: cpp-linter/cpp-linter-action@v2 - id: linter - env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + - 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 src ! -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 'files<> "$GITHUB_OUTPUT" + else + echo 'files=' >> "$GITHUB_OUTPUT" + fi + + - name: Comment on PR if formatting issues found + if: steps.format_check.outputs.files != '' + uses: peter-evans/create-or-update-comment@v4 with: - style: 'file' # Use .clang-format config file. - tidy-checks: '-*' # disable clang-tidy checks. - extensions: 'cu,cuh,c,h,C,H,cpp,hpp,cc,hh,c++,h++,cxx,hxx' - files-changed-only: false - # only 'update' a single comment in a pull request's thread. - thread-comments: ${{ github.event_name == 'pull_request' && 'update' }} - - - name: Fail fast?! - if: steps.linter.outputs.clang-format-checks-failed > 0 - run: exit 1 + issue-number: ${{ github.event.pull_request.number }} + body: | + ⚠️ **Clang-format check failed** + + The following files are not correctly formatted: + ``` + "${{ steps.format_check.outputs.files }}" + ``` + Please run `clang-format -i` on these files and commit the changes before merging. + + + - name: Fail job if formatting issues found + if: steps.format_check.outputs.files != '' + run: | + echo "Formatting issues detected. Failing the job." + exit 1 + From b83d0283ee198de89d3273f9b7f9370db0d6a3b9 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 18:01:12 +0200 Subject: [PATCH 14/24] style: try another workflow --- .github/workflows/format.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index 0ce51992..fe78be4f 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -19,7 +19,7 @@ jobs: id: format_check run: | FILE_LIST="unformatted_files.txt" - find src ! -path 'src/third_party/*' \ + 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" From b4f45c95c6e6edd7ceed08a39ec298880a7f1969 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 18:05:07 +0200 Subject: [PATCH 15/24] style: add CI for python formatting --- .github/workflows/format-python.yml | 66 +++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) create mode 100644 .github/workflows/format-python.yml diff --git a/.github/workflows/format-python.yml b/.github/workflows/format-python.yml new file mode 100644 index 00000000..e905b97d --- /dev/null +++ b/.github/workflows/format-python.yml @@ -0,0 +1,66 @@ +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 capture output + id: black_check + continue-on-error: true + run: | + set +e + OUTPUT=$(black --check --diff . 2>&1) + EXIT_CODE=$? + + echo "$OUTPUT" + echo "$OUTPUT" > black_output.txt + echo "exit_code=$EXIT_CODE" >> $GITHUB_OUTPUT + + if [[ "$EXIT_CODE" -ne 0 ]]; then + echo "py_failed=true" >> $GITHUB_ENV + else + echo "py_failed=false" >> $GITHUB_ENV + fi + + - name: Comment on PR if Black 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 }} + body: | + ⚠️ **Black formatting check failed** + + The following Python files are not correctly formatted: + + ``` + ${{ steps.black_check.outputs.black_output }} + ``` + + Please run: + ```bash + black . + ``` + And commit the changes before merging. + + - name: Fail if Black failed + if: env.py_failed == 'true' + run: | + echo "❌ Python formatting issues found." + exit 1 From 38cb8d1e5589c70708280f0b1d869b09307af345 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 18:10:56 +0200 Subject: [PATCH 16/24] style: add CI for python formatting --- .github/workflows/format-python.yml | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/.github/workflows/format-python.yml b/.github/workflows/format-python.yml index e905b97d..0b74f029 100644 --- a/.github/workflows/format-python.yml +++ b/.github/workflows/format-python.yml @@ -20,22 +20,26 @@ jobs: - name: Install Black run: pip install black - - name: Run Black and capture output + - name: Run Black and extract offending files id: black_check continue-on-error: true run: | set +e - OUTPUT=$(black --check --diff . 2>&1) + OUTPUT=$(black --check . 2>&1) EXIT_CODE=$? - echo "$OUTPUT" echo "$OUTPUT" > black_output.txt - echo "exit_code=$EXIT_CODE" >> $GITHUB_OUTPUT - if [[ "$EXIT_CODE" -ne 0 ]]; then + grep '^would reformat' black_output.txt | awk '{print $3}' | 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: Comment on PR if Black failed @@ -47,10 +51,10 @@ jobs: body: | ⚠️ **Black formatting check failed** - The following Python files are not correctly formatted: + The following Python files need formatting: ``` - ${{ steps.black_check.outputs.black_output }} + ${{ steps.black_check.outputs.black_files }} ``` Please run: From e4a0546006166dab68cdceebee9b62da8865abbc Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 18:17:13 +0200 Subject: [PATCH 17/24] style: black formatting --- setup.py | 8 ++++++-- tests/test_array_origins.py | 4 +++- tests/test_interface.py | 3 +++ tests/test_thermal_drift.py | 23 ++++++++++++----------- 4 files changed, 24 insertions(+), 14 deletions(-) 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/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 From aa1c33ec08d207d68753fe7e2e16ee142e69852b Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Wed, 25 Jun 2025 18:18:26 +0200 Subject: [PATCH 18/24] style: clang-format formatting --- include/MobilityInterface/MobilityInterface.h | 4 +- include/MobilityInterface/lanczos.h | 33 +- include/MobilityInterface/pythonify.h | 86 +- include/memory/allocator.h | 3 +- include/memory/python_tensor.h | 2 +- solvers/DPStokes/extra/poly_fits.h | 74 +- solvers/DPStokes/extra/uammd_interface.h | 111 ++- solvers/DPStokes/extra/uammd_wrapper.cu | 437 ++++---- solvers/DPStokes/python_wrapper.cu | 33 +- solvers/NBody/extra/NbodyRPY.cu | 80 +- solvers/NBody/extra/hydrodynamicKernels.cuh | 201 ++-- solvers/NBody/extra/interface.h | 6 +- solvers/NBody/extra/vector.cuh | 939 ++++++++---------- solvers/NBody/python_wrapper.cu | 12 +- solvers/PSE/extra/uammd_interface.h | 66 +- solvers/PSE/extra/uammd_wrapper.cu | 2 +- solvers/PSE/python_wrapper.cu | 4 +- solvers/SelfMobility/python_wrapper.cu | 24 +- 18 files changed, 993 insertions(+), 1124 deletions(-) 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/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/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/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/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"); From ede5b97c220f133010846f70e8358c2ce120cff2 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Thu, 26 Jun 2025 10:48:09 +0200 Subject: [PATCH 19/24] Add ignore blame revs --- .git-blame-ignore-revs | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 .git-blame-ignore-revs 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 From df63b3b2f6d189e76154a7f18749547083106d12 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Thu, 26 Jun 2025 11:03:43 +0200 Subject: [PATCH 20/24] style: update python linter CI --- .github/workflows/format-python.yml | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/.github/workflows/format-python.yml b/.github/workflows/format-python.yml index 0b74f029..5d22aa42 100644 --- a/.github/workflows/format-python.yml +++ b/.github/workflows/format-python.yml @@ -42,13 +42,24 @@ jobs: echo "black_files=" >> $GITHUB_OUTPUT fi - - name: Comment on PR if Black failed + - 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: @@ -63,6 +74,20 @@ jobs: ``` 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: | From 0f08fb1fdd540b6db05f382257719a7c5b25a7be Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Thu, 26 Jun 2025 11:05:01 +0200 Subject: [PATCH 21/24] Add test file to trip the python linter --- test.py | 1 + 1 file changed, 1 insertion(+) create mode 100644 test.py diff --git a/test.py b/test.py new file mode 100644 index 00000000..024a9fec --- /dev/null +++ b/test.py @@ -0,0 +1 @@ +import numpy From fd00512b91022a74f355ea939f43146a22b5aaef Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Thu, 26 Jun 2025 11:08:05 +0200 Subject: [PATCH 22/24] style: Remove absolute path --- .github/workflows/format-python.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.github/workflows/format-python.yml b/.github/workflows/format-python.yml index 5d22aa42..f577e7cf 100644 --- a/.github/workflows/format-python.yml +++ b/.github/workflows/format-python.yml @@ -29,8 +29,7 @@ jobs: EXIT_CODE=$? echo "$OUTPUT" > black_output.txt - - grep '^would reformat' black_output.txt | awk '{print $3}' | sort -u > black_files.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 From 85c01fa9369748a1cd81793d409a0a4bdab4a550 Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Thu, 26 Jun 2025 11:09:10 +0200 Subject: [PATCH 23/24] Delete test file to trip the formatting ci --- test.py | 1 - 1 file changed, 1 deletion(-) delete mode 100644 test.py diff --git a/test.py b/test.py deleted file mode 100644 index 024a9fec..00000000 --- a/test.py +++ /dev/null @@ -1 +0,0 @@ -import numpy From c0792ae9f6994bf28a12e041e8b53826a01be53d Mon Sep 17 00:00:00 2001 From: "Raul P. Pelaez" Date: Thu, 26 Jun 2025 11:11:13 +0200 Subject: [PATCH 24/24] style: update clang linter CI --- .github/workflows/format.yml | 54 +++++++++++++++++++++++++++++------- 1 file changed, 44 insertions(+), 10 deletions(-) diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index fe78be4f..5f1f785e 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -22,35 +22,69 @@ jobs: 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 + 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: Comment on PR if formatting issues found - if: steps.format_check.outputs.files != '' + - 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 }}" + ${{ steps.format_check.outputs.files }} ``` - Please run `clang-format -i` on these files and commit the changes before merging. - + + 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: steps.format_check.outputs.files != '' + if: env.failed == 'true' run: | - echo "Formatting issues detected. Failing the job." + echo "❌ Formatting issues detected. Failing the job." exit 1 -