Skip to content
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

refactor: Use matrix-matrix interface in LAPACK solve #3234

Merged
merged 9 commits into from
Aug 26, 2024
Merged

refactor: Use matrix-matrix interface in LAPACK solve #3234

Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
1e3f671
Use matrix interface
dkachuma Jul 23, 2024
6b8e3ed
uncrustify
dkachuma Jul 23, 2024
1bdbda2
Remove unrelated changes
dkachuma Jul 23, 2024
562d8b8
Remove unncessesary assignment
dkachuma Jul 23, 2024
c2c8fab
Merge branch 'develop' into refactor/dkachuma/use-lapack-matrix-matri…
dkachuma Aug 10, 2024
e5867d4
Merge branch 'develop' into refactor/dkachuma/use-lapack-matrix-matri…
dkachuma Aug 12, 2024
006e81e
Merge branch 'develop' into refactor/dkachuma/use-lapack-matrix-matri…
dkachuma Aug 13, 2024
5fee5af
Merge branch 'develop' into refactor/dkachuma/use-lapack-matrix-matri…
dkachuma Aug 16, 2024
bdf0c0d
Merge branch 'develop' into refactor/dkachuma/use-lapack-matrix-matri…
CusiniM Aug 23, 2024
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
72 changes: 33 additions & 39 deletions ...omponents/constitutive/fluid/multifluid/compositional/functions/NegativeTwoPhaseFlash.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -195,23 +195,21 @@ struct NegativeTwoPhaseFlash

/**
* @brief Solve the lineat system for the derivatives of the flash
* @param[in] A the coefficient matrix
* @param[in] b the rhs
* @param[out] x the solution
* @param[in/out] A the coefficient matrix. Destroyed after call
* @param[in/out] X the rhs and solution
* @return @c true if the problem is well solved @c false otherwise
*/
template< int USD >
GEOS_HOST_DEVICE
static bool solveLinearSystem( arraySlice2d< real64 const > const & A,
arraySlice1d< real64 const > const & b,
arraySlice1d< real64 > const & x )
static bool solveLinearSystem( arraySlice2d< real64, USD > const & A,
Copy link
Collaborator

Choose a reason for hiding this comment

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

I personally preferred having the rhs and the solution into two separate arrays. I am guessing that this is the lapack interface though... We could keep the separation on the GEOS side though.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The separate rhs/solution interface on BlasLapackLA still exists. This change destroys the rhs but avoids having to allocate extra space for solution and rhs knowing that we are not going to need the rhs or the matrix after the solve.

Copy link
Collaborator

Choose a reason for hiding this comment

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

is the allocation the same? A single array2d instead of 2 array1d?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

My mistake. The allocation was on the matrix. With the previous interface, everytime we call the solve, additional space had to be allocated for the LU factors of the matrix. This way there is no extra space allocated since it's a single solve and the space allocated for the matrix is used for the LU factors (destroying the matrix but once we have solved we don't care anymore).

With the single column by column solve it would have been nice to factorize the matrix and then use the LU factors on each column but we don't have that interface. Instead I added an interface where you can do an in-place matrix-matrix solve.

arraySlice2d< real64, USD > const & X )
{
#if defined(GEOS_DEVICE_COMPILE)
GEOS_UNUSED_VAR( A );
GEOS_UNUSED_VAR( b );
GEOS_UNUSED_VAR( x );
GEOS_UNUSED_VAR( X );
return false;
#else
BlasLapackLA::solveLinearSystem( A, b, x );
BlasLapackLA::solveLinearSystem( A, X );
return true;
#endif
}
Expand Down Expand Up @@ -432,12 +430,11 @@ void NegativeTwoPhaseFlash::computeDerivatives(

constexpr integer maxNumVals = 2*MultiFluidConstants::MAX_NUM_COMPONENTS+1;
integer const numVals = 2*numComps;
stackArray1d< real64, maxNumVals > b( numVals + 1 );
stackArray1d< real64, maxNumVals > x( numVals + 1 );
stackArray2d< real64, maxNumVals * maxNumVals > A( numVals + 1, numVals + 1 );
StackArray< real64, 2, maxNumVals * maxNumVals, MatrixLayout::COL_MAJOR_PERM > A( numVals + 1, numVals + 1 );
StackArray< real64, 2, maxNumVals * maxNumVals, MatrixLayout::COL_MAJOR_PERM > X( numVals + 1, numVals + 1 );

LvArray::forValuesInSlice( A.toSlice(), setZero );
LvArray::forValuesInSlice( b.toSlice(), setZero );
LvArray::forValuesInSlice( X.toSlice(), setZero );

for( integer ic = 0; ic < numComps; ++ic )
{
Expand Down Expand Up @@ -469,45 +466,42 @@ void NegativeTwoPhaseFlash::computeDerivatives(
A( e, xi ) = -1.0;
A( e, yi ) = 1.0;
}

// Pressure and temperature derivatives
for( integer const pc : {Deriv::dP, Deriv::dT} )
for( integer ic = 0; ic < numComps; ++ic )
{
for( integer ic = 0; ic < numComps; ++ic )
{
real64 const phiL = exp( logLiquidFugacity( ic ) );
real64 const phiV = exp( logVapourFugacity( ic ) );
b( ic ) = 0.0;
b( ic + numComps ) = -liquidComposition[ic] * phiL * logLiquidFugacityDerivs( ic, pc )
+ vapourComposition[ic] * phiV * logVapourFugacityDerivs( ic, pc );
}
b( numVals ) = 0.0;
solveLinearSystem( A, b, x );
for( integer ic = 0; ic < numComps; ++ic )
{
liquidCompositionDerivs( ic, pc ) = x( ic );
vapourCompositionDerivs( ic, pc ) = x( ic + numComps );
}
vapourFractionDerivs( pc ) = x( numVals );
real64 const phiL = -liquidComposition[ic] * exp( logLiquidFugacity( ic ) );
real64 const phiV = vapourComposition[ic] * exp( logVapourFugacity( ic ) );
X( ic + numComps, Deriv::dP ) = phiL * logLiquidFugacityDerivs( ic, Deriv::dP ) +
phiV * logVapourFugacityDerivs( ic, Deriv::dP );
X( ic + numComps, Deriv::dT ) = phiL * logLiquidFugacityDerivs( ic, Deriv::dT ) +
phiV * logVapourFugacityDerivs( ic, Deriv::dT );
}

// Composition derivatives
for( integer kc = 0; kc < numComps; ++kc )
{
integer const pc = Deriv::dC + kc;
integer const idof = Deriv::dC + kc;

for( integer ic = 0; ic < numComps; ++ic )
{
b( ic ) = -composition[ic];
b( ic + numComps ) = 0.0;
X( ic, idof ) = -composition[ic];
}
b( kc ) += 1.0;
b( numVals ) = 0.0;
solveLinearSystem( A, b, x );
X( kc, idof ) += 1.0;
}

// Solve linear system
solveLinearSystem( A.toSlice(), X.toSlice() );

// Fill in the derivatives
for( integer idof = 0; idof < numDofs; ++idof )
{
for( integer ic = 0; ic < numComps; ++ic )
{
liquidCompositionDerivs( ic, pc ) = x( ic );
vapourCompositionDerivs( ic, pc ) = x( ic + numComps );
liquidCompositionDerivs( ic, idof ) = X( ic, idof );
vapourCompositionDerivs( ic, idof ) = X( ic + numComps, idof );
}
vapourFractionDerivs( pc ) = x( numVals );
vapourFractionDerivs( idof ) = X( numVals, idof );
}
}
}
Expand Down