Skip to content

Combined multifluid derivatives to reduce fluid kernel wrapper size #1779

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
francoishamon merged 10 commits into from
Feb 23, 2022
Merged

Combined multifluid derivatives to reduce fluid kernel wrapper size #1779

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
c928262
combined multifluid derivatives
francoishamon Feb 15, 2022
73e01aa
fixed the testMultiFluid unit test
francoishamon Feb 17, 2022
33f5c99
Merge branch 'develop' into refactor/hamon/multifluidDerivatives
francoishamon Feb 17, 2022
6b2d7ca
Merge branch 'develop' into refactor/hamon/multifluidDerivatives
francoishamon Feb 17, 2022
1609e7c
addressed review comments; no numerical diffs in integratedTests
francoishamon Feb 18, 2022
a5b5bd9
fixed unit test
francoishamon Feb 18, 2022
3e72aad
fixed debug output
francoishamon Feb 19, 2022
fbf4f5d
Merge branch 'develop' into refactor/hamon/multifluidDerivatives
francoishamon Feb 22, 2022
3910b37
Merge branch 'develop' into refactor/hamon/multifluidDerivatives
francoishamon Feb 23, 2022
b9b88e3
updated integratedTests submodule
francoishamon Feb 23, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion integratedTests
Open in desktop
Submodule integratedTests updated from 2e4590 to 25a069
50 changes: 43 additions & 7 deletions src/coreComponents/codingUtilities/Utilities.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -186,43 +186,79 @@ void forEachArgInTuple( std::tuple< Ts ... > const & tuple, F && func )

// The code below should work with any subscriptable vector/matrix types

/**
* @brief Utility function to copy a vector into another vector
* @tparam VEC1 the type of the source vector
* @tparam VEC2 the type of the destination vector
* @param[in] N the number of values to copy
* @param[in] v1 the source vector
* @param[out] v2 the destination vector
* @param[in] offset the first index of v2 at which we start copying values
*/
template< typename VEC1, typename VEC2 >
GEOSX_HOST_DEVICE
void copy( integer const N, VEC1 const & v1, VEC2 const & v2 )
void copy( integer const N,
Copy link
Contributor

@TotoGaz TotoGaz Feb 17, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry to bother you with this, but documenting why we need this firstDeriv could be helpful (in which context we need this, e.g.)

More important, instead of providing firstDeriv, couldn't v2 be a range/slice instead? That would keep this function more mathematical and less implementation aware.

Copy link
Contributor Author

@francoishamon francoishamon Feb 19, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, as we discussed, it would be nice to passe a slice or range, but I don't know how to do it.
Here for example:
https://github.com/GEOSX/GEOSX/blob/3e72aadf62b7083dd7f434052f3d90c8f3470230/src/coreComponents/constitutive/fluid/MultiFluidBase.hpp#L722
my phaseFrac.derivs[ip] is already a slice of size numComp+2, but I don't need first two entries (derivatives with respect to pressure and temperature). I would need a way to skip them, maybe by doing:

applyChainRuleInPlace( numComp, dCompMoleFrac_dCompMassFrac, phaseFrac.derivs[ip][Deriv::dC], work );

Now that I am writing this, I should probably try that to see if it works.

For now, I added the Doxygen documentation in the applyChainRule function to explain why we have to pass Deriv::dC.

Copy link
Contributor

@klevzoff klevzoff Feb 19, 2022 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In this case v2 represents a general single-subscriptable entity, which may be a local array or an ArraySlice. The latter may be not contiguous in memory, and in that case precomputing an offset pointer is not safe or correct. So we have to rely on the type's operator[] to do the job at all times.

Theoretically it's possible to write a SubSlice type which represents a sub-selection of an ArraySlice along one or more dimensions. The type would store the offset(s) internally and apply as needed. But we don't have it yet.

VEC1 const & v1,
VEC2 const & v2,
integer const offset = 0 )
{
for( integer i = 0; i < N; ++i )
{
v2[i] = v1[i];
v2[offset+i] = v1[i];
}
}

/**
* @brief Utility function to apply the chain rule to compute df_dx as a function of df_dy and dy_dx
* @tparam MATRIX the type of the matrix of derivatives
* @tparam VEC1 the type of the source vector
* @tparam VEC2 the type of the destination vector
* @param[in] N the number of derivative values
* @param[in] dy_dx the matrix of derivatives of y(x) wrt x
* @param[in] df_dy the derivatives of f wrt y
* @param[out] df_dx the computed derivatives of f wrt x
* @param[in] firstDerivativeOffset the first derivative offset in df_dy
*/
template< typename MATRIX, typename VEC1, typename VEC2 >
GEOSX_HOST_DEVICE
void applyChainRule( integer const N,
MATRIX const & dy_dx,
VEC1 const & df_dy,
VEC2 && df_dx )
VEC2 && df_dx,
integer const firstDerivativeOffset = 0 )
{
// this could use some dense linear algebra
for( integer i = 0; i < N; ++i )
{
df_dx[i] = 0.0;
for( integer j = 0; j < N; ++j )
{
df_dx[i] += df_dy[j] * dy_dx[j][i];
df_dx[i] += df_dy[firstDerivativeOffset+j] * dy_dx[j][i];
}
}
}

/**
* @brief Utility function to apply the chain rule to compute df_dxy in place
* @tparam MATRIX the type of the matrix of derivatives
* @tparam VEC1 the type of the source vector
* @tparam VEC2 the type of the utility vector
* @param[in] N the number of derivative values
* @param[in] dy_dx the matrix of derivatives of y(x) wrt x
* @param[out] df_dxy the derivatives of f wrt y
* @param[out] work utility array
* @param[in] offset the first derivative offset in df_dy
*/
template< typename MATRIX, typename VEC1, typename VEC2 >
GEOSX_HOST_DEVICE
void applyChainRuleInPlace( integer const N,
MATRIX const & dy_dx,
VEC1 && df_dxy,
VEC2 && work )
VEC2 && work,
integer const firstDeriv = 0 )
{
applyChainRule( N, dy_dx, df_dxy, work );
copy( N, work, df_dxy );
applyChainRule( N, dy_dx, df_dxy, work, firstDeriv );
copy( N, work, df_dxy, firstDeriv );
}

} // namespace geosx
Expand Down
3 changes: 1 addition & 2 deletions src/coreComponents/constitutive/fluid/BlackOilFluid.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ void BlackOilFluid::readInputDataFromPVTFiles()
fillWaterData( data );

// gas data
fillHydrocarbonData( PT::GAS, boTables.getGasTable() );
fillHydrocarbonData( m_phaseOrder[PT::GAS], boTables.getGasTable() );

// for the Black-Oil model, the oil PVT is treated differently from gas
fillPVTOData( boTables.getOilTable(),
Expand Down Expand Up @@ -110,7 +110,6 @@ void BlackOilFluid::fillPVTOData( array1d< array1d< real64 > > const & oilTable,
// if oilTable.size() == 4 (saturated case): Rs, bubble point pressure, Bo, viscosity
// if ollTable.size() == 3 (unsaturated case): unsaturated pressure, Bo, viscosity


// Step 1: count the number of saturated points by looping through oilTable, and resize tables accordingly

integer numSaturatedPoints = 0;
Expand Down
Loading