Skip to content

Commit

Permalink
Intrepid2: add support for weighted gradients in sum-factorized assem…
Browse files Browse the repository at this point in the history
…bly (#13391)

Intrepid2: 
This PR adds support for additional uses cases to the existing sum-factorized assembly mechanism.  On the way, it adds support for vector weights to `TransformedBasisValues`; previously, only matrix and scalar weights were supported.

This PR also makes some minor adjustments to the timings performed in `StructuredIntegrationPerformance.cpp` to ensure fair comparisons.

The new use cases are demonstrated most directly in the following:

```
packages/intrepid2/assembly-examples/VectorWeightedGRADGRADStructuredAssembly.hpp
```

which demonstrates performing an integral of the form 

$`({\mathbf a} \cdot \nabla e_i, {\mathbf b} \cdot \nabla e_j)`$

for an $H^1$ basis.

## Testing
The new use cases are well-exercised by the tests; a new set of `VectorWeightedPoisson` test cases (corresponding to the integral above) are included in the structured-versus-standard tests.

## Other notes/details:
- Added support for dot products to Intrepid2_Data; added a corresponding test.  Also hardened the test MatVec_CPDD_transpose to include a check that u' A v = v' A' u for vectors u, v.
- Added free function rank() taking BasisValues object as argument.
- Data: fixed an issue in allocateMatVec in which an incorrect variation type and/or incorrect extent could be used for the final result dimension.
- TransformedBasisValues: added support for a (C,P,D) transform, with the main use case being a dot product with a vector-valued basis evaluation.
- In allocateMatMatResult(), fixed an issue in which the wrong getUnderlyingView() method was being called; resolved by calling the one that gets a DynRankView.
- Added a test against taking the outer product of two vectors; fixed the issue that this demonstrated in Intrepid2::Data.
- Added test template StructuredVersusStandardVectorWeighted_D2_P1_P1.
- Fixed an issue with the transpose arguments to a mat-mat call.
- In setJacobianDetInv(), corrected argument name in method declaration and doxygen.
- In DataTools, broadened the use cases for multiplyByCPWeights(), and added a transposeMatrix() method.
- In TransformedBasisValues, fixed some issues with spaceDim().
- In StandardAssembly and StructuredAssembly, revised to support more vector-weighted use cases.
- Added more vector-weighted tests, covering cases when a vector field is dotted with a vector and then integrated against a scalar.
- Modified standard assembly performance tests to exclude the orientation application from "core integration" timing.
- Added VectorWeightedPoisson to allFormulationChoices.  For now, setting the "best" CUDA choices to match Poisson.  But we need to redo those calibrations regardless, and maybe switch to reading them in from file: these are for older CUDA cards.
  • Loading branch information
CamelliaDPG authored Aug 27, 2024
1 parent 3132245 commit 9c6525b
Show file tree
Hide file tree
Showing 24 changed files with 2,362 additions and 374 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -145,8 +145,11 @@ Intrepid2::ScalarView<Scalar,DeviceType> performStandardQuadratureGRADGRAD(Intre
// because structured integration performs transformations within integrate(), to get a fairer comparison here we include the transformation calls.
fstIntegrateCall->start();
FunctionSpaceTools::HGRADtransformGRAD(unorientedTransformedGradValues, jacobianInverse, basisGradValues);
// we want to exclude orientation application in the core integration timing -- this time gets reported as "Other"
fstIntegrateCall->stop();
OrientationTools<DeviceType>::modifyBasisByOrientation(transformedGradValues, unorientedTransformedGradValues,
orientationsWorkset, basis.get());
fstIntegrateCall->start();

transformIntegrateFlopCount += double(numCellsInWorkset) * double(numFields) * double(numPoints) * double(spaceDim) * (spaceDim - 1) * 2.0; // 2: one multiply, one add per (P,D) entry in the contraction.
FunctionSpaceTools::multiplyMeasure(transformedWeightedGradValues, cellMeasures, transformedGradValues);
Expand Down
6 changes: 6 additions & 0 deletions packages/intrepid2/assembly-examples/H1StandardAssembly.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -151,8 +151,11 @@ Intrepid2::ScalarView<Scalar,DeviceType> performStandardQuadratureH1(Intrepid2::
// because structured integration performs transformations within integrate(), to get a fairer comparison here we include the transformation calls.
fstIntegrateCall->start();
FunctionSpaceTools::HGRADtransformGRAD(unorientedTransformedGradValues, jacobianInverse, basisGradValues);
// we want to exclude orientation application in the core integration timing -- this time gets reported as "Other"
fstIntegrateCall->stop();
OrientationTools<DeviceType>::modifyBasisByOrientation(transformedGradValues, unorientedTransformedGradValues,
orientationsWorkset, basis.get());
fstIntegrateCall->start();
transformIntegrateFlopCount += double(numCellsInWorkset) * double(numFields) * double(numPoints) * double(spaceDim) * (spaceDim - 1) * 2.0; // 2: one multiply, one add per (P,D) entry in the contraction.
FunctionSpaceTools::multiplyMeasure(transformedWeightedGradValues, cellMeasures, transformedGradValues);
transformIntegrateFlopCount += double(numCellsInWorkset) * double(numFields) * double(numPoints) * double(spaceDim); // multiply each entry of transformedGradValues: one flop for each.
Expand All @@ -163,8 +166,11 @@ Intrepid2::ScalarView<Scalar,DeviceType> performStandardQuadratureH1(Intrepid2::
ExecutionSpace().fence();

FunctionSpaceTools::HGRADtransformVALUE(unorientedTransformedBasisValues, basisValues);
// we want to exclude orientation application in the core integration timing -- this time gets reported as "Other"
fstIntegrateCall->stop();
OrientationTools<DeviceType>::modifyBasisByOrientation(transformedBasisValues, unorientedTransformedBasisValues,
orientationsWorkset, basis.get());
fstIntegrateCall->start();
FunctionSpaceTools::multiplyMeasure(transformedWeightedBasisValues, cellMeasures, transformedBasisValues);
bool sumInto = true; // add the (value,value) integral to the (grad,grad) that we've already integrated
FunctionSpaceTools::integrate(cellStiffnessSubview, transformedBasisValues, transformedWeightedBasisValues, sumInto);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -175,8 +175,11 @@ Intrepid2::ScalarView<Scalar,DeviceType> performStandardQuadratureHCURL(Intrepid
// because structured integration performs transformations within integrate(), to get a fairer comparison here we include the transformation calls.
fstIntegrateCall->start();
FunctionSpaceTools::HCURLtransformCURL(unorientedTransformedCurlValues, jacobian, jacobianDeterminant, basisCurlValues);
// we want to exclude orientation application in the core integration timing -- this time gets reported as "Other"
fstIntegrateCall->stop();
OrientationTools<DeviceType>::modifyBasisByOrientation(transformedCurlValues, unorientedTransformedCurlValues,
orientationsWorkset, basis.get());
fstIntegrateCall->start();
transformIntegrateFlopCount += double(numCellsInWorkset) * double(numFields) * double(numPoints) * double(spaceDim) * (spaceDim - 1) * 2.0; // 2: one multiply, one add per (P,D) entry in the contraction.
FunctionSpaceTools::multiplyMeasure(transformedWeightedCurlValues, cellMeasures, transformedCurlValues);
transformIntegrateFlopCount += double(numCellsInWorkset) * double(numFields) * double(numPoints) * double(spaceDim); // multiply each entry of transformedCurlValues: one flop for each.
Expand All @@ -186,8 +189,11 @@ Intrepid2::ScalarView<Scalar,DeviceType> performStandardQuadratureHCURL(Intrepid
FunctionSpaceTools::integrate(cellStiffnessSubview, transformedCurlValues, transformedWeightedCurlValues);

FunctionSpaceTools::HCURLtransformVALUE(unorientedTransformedBasisValues, jacobianInverse, basisValues);
// we want to exclude orientation application in the core integration timing -- this time gets reported as "Other"
fstIntegrateCall->stop();
OrientationTools<DeviceType>::modifyBasisByOrientation(transformedBasisValues, unorientedTransformedBasisValues,
orientationsWorkset, basis.get());
fstIntegrateCall->start();
FunctionSpaceTools::multiplyMeasure(transformedWeightedBasisValues, cellMeasures, transformedBasisValues);
bool sumInto = true; // add the (value,value) integral to the (curl,curl) that we've already integrated
FunctionSpaceTools::integrate(cellStiffnessSubview, transformedBasisValues, transformedWeightedBasisValues, sumInto);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -151,8 +151,11 @@ Intrepid2::ScalarView<Scalar,DeviceType> performStandardQuadratureHDIV(Intrepid2
// because structured integration performs transformations within integrate(), to get a fairer comparison here we include the transformation calls.
fstIntegrateCall->start();
FunctionSpaceTools::HDIVtransformDIV(unorientedTransformedDivValues, jacobianDeterminant, basisDivValues);
// we want to exclude orientation application in the core integration timing -- this time gets reported as "Other"
fstIntegrateCall->stop();
OrientationTools<DeviceType>::modifyBasisByOrientation(transformedDivValues, unorientedTransformedDivValues,
orientationsWorkset, basis.get());
fstIntegrateCall->start();
transformIntegrateFlopCount += double(numCellsInWorkset) * double(numFields) * double(numPoints) * double(spaceDim) * (spaceDim - 1) * 2.0; // 2: one multiply, one add per (P,D) entry in the contraction.
FunctionSpaceTools::multiplyMeasure(transformedWeightedDivValues, cellMeasures, transformedDivValues);
transformIntegrateFlopCount += double(numCellsInWorkset) * double(numFields) * double(numPoints) * double(spaceDim); // multiply each entry of transformedDivValues: one flop for each.
Expand All @@ -161,10 +164,12 @@ Intrepid2::ScalarView<Scalar,DeviceType> performStandardQuadratureHDIV(Intrepid2

FunctionSpaceTools::integrate(cellStiffnessSubview, transformedDivValues, transformedWeightedDivValues);
ExecutionSpace().fence();

FunctionSpaceTools::HDIVtransformVALUE(unorientedTransformedBasisValues, jacobian, jacobianDeterminant, basisValues);
// we want to exclude orientation application in the core integration timing -- this time gets reported as "Other"
fstIntegrateCall->stop();
OrientationTools<DeviceType>::modifyBasisByOrientation(transformedBasisValues, unorientedTransformedBasisValues,
orientationsWorkset, basis.get());
fstIntegrateCall->start();
FunctionSpaceTools::multiplyMeasure(transformedWeightedBasisValues, cellMeasures, transformedBasisValues);
bool sumInto = true; // add the (value,value) integral to the (div,div) that we've already integrated
FunctionSpaceTools::integrate(cellStiffnessSubview, transformedBasisValues, transformedWeightedBasisValues, sumInto);
Expand Down
3 changes: 3 additions & 0 deletions packages/intrepid2/assembly-examples/HVOLStandardAssembly.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -139,8 +139,11 @@ Intrepid2::ScalarView<Scalar,DeviceType> performStandardQuadratureHVOL(Intrepid2
auto cellStiffnessSubview = Kokkos::subview(cellStiffness, cellRange, Kokkos::ALL(), Kokkos::ALL());

FunctionSpaceTools::HVOLtransformVALUE(unorientedTransformedBasisValues, jacobianDeterminant, basisValues);
// we want to exclude orientation application in the core integration timing -- this time gets reported as "Other"
fstIntegrateCall->stop();
OrientationTools<DeviceType>::modifyBasisByOrientation(transformedBasisValues, unorientedTransformedBasisValues,
orientationsWorkset, basis.get());
fstIntegrateCall->start();
FunctionSpaceTools::multiplyMeasure(transformedWeightedBasisValues, cellMeasures, transformedBasisValues);
bool sumInto = true; // add the (value,value) integral to the (curl,curl) that we've already integrated
FunctionSpaceTools::integrate(cellStiffnessSubview, transformedBasisValues, transformedWeightedBasisValues, sumInto);
Expand Down
Loading

0 comments on commit 9c6525b

Please sign in to comment.