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.

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

Implement relative permeabilities with tables on device #1277

Merged
merged 14 commits into from
Feb 19, 2021

Conversation

francoishamon
Copy link
Contributor

@francoishamon francoishamon commented Jan 19, 2021

This PR adds the option to compute relative permeabilities on device with tables provided as input in the xml file.

The design follows the plan outlined by @klevzoff in issue #1232:

  • The PR adds a TableFunction::KernelWrapper nested in TableFunction.hpp to encapsulates the views necessary to interpolate in the table. This wrapper has a device-callable compute function that interpolates in the table and returns the value and the derivatives.
  • The PR modifies the TableFunction::evaluate function, which is now just calling TableFunction::KernelWrapper::compute, then returns the interpolated value and discards the derivatives.
  • The PR adds a TableRelativePermeability class that has an instance of TableFunction::KernelWrapper for each phase, and uses it to interpolate in the table to get the relperm and its derivative w.r.t. phase volume fraction.

To make the implementation suitable for the specific use case that we have (three-phase Black-Oil), I will have to add a three-phase relperm interpolation (Baker, Stone, etc). This may require some restructuring in the relperm folder to more easily reuse what has already been implemented for the other models.

  • Initial implementation working (checked on Quartz only).
  • Add three-phase interpolation method (Baker).
  • Add unit test for derivatives in TableFunction and TableRelativePermeability.
  • Add integrated test.

Resolves #1232
Resolves #1226

Related to https://github.com/GEOSX/integratedTests/pull/117

Copy link
Contributor

@klevzoff klevzoff left a comment

Choose a reason for hiding this comment

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

Looks awesome so far!

Comment on lines 183 to 187
// Create the kernel wrapper
// Warning! This function must be called AFTER all the input parameters have been set (interpolation method, etc)
// TODO: see if there is a better way, this seems quite error prone
m_kernelWrapper =
std::make_unique< TableFunctionKernelWrapper >( m_interpolationMethod,
Copy link
Contributor

Choose a reason for hiding this comment

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

I don't think there's a way to avoid having to recreate the kernel wrapper every time the table object is modified. Luckily, this doesn't happen much outside of unit tests (I think?). If only the coordinates and values change (but not the number of points in each dir), then array views will see the updated data automatically, but on any array resize the views need to be recreated.

You could get rid of the unique_ptr if you wanted to (if you gave a default constructor to TableFunctionKernelWrapper) and keep the wrapper object as a direct member of TableFunction. This would avoid an extra indirection in Evaluate. I'm ok with either solution 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.

Yes you are right I only had this problem in the testFunctions unit test. Now I recreate the kernel wrapper after each call to the setters.

After adding the default constructor to TableFunctionKernelWrapper, I did not find a good way to have the wrapper object as a direct member of TableFunction to remove the pointer:

  • Currently, in the hpp file, I do a forward declaration of class TableFunctionKernelWrapper, and then there is the definition of the class TableFunction, and then that of the class TableFunctionKernelWrapper. So when I have the wrapper object as a direct member of TableFunction, the compiler says that the field 'm_kernelWrapper' has incomplete type.
  • If I try the other approach: forward declare class TableFunction, then define the class TableFunctionKernelWrapper, and finally the class TableFunction, I get similar compilation problems with some types and variables needed in TableFunctionKernelWrapper, like TableFunction::InterpolationType and TableFunction::maxDimensions.

Is there a simple workaround for these errors? (I left the pointer in the current version for now).

Copy link
Contributor

Choose a reason for hiding this comment

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

Hmm, I see. Indeed, we have a circular dependency. The dependency of TableFunction on TableFunctionKernelWrapper is strong (does require the full definition if contains a member variable), but the other direction is weak (it only needs some nested constants and typedefs). Normally it would be resolved one of two ways:

  • taking those constants and typedefs out of TableFunction, which I already advised against (I definitely like the proper scoping of those);
  • or making TableFunctionKernelWrapper a nested (inner) class of TableFunction, which seems like a good idea to me, except we currently don't follow that pattern in any of the constitutive models that have a dedicated "update"/"kernel wrapper" class. In that case, it can lose the TableFunction part of its name and become simply KernelWrapper.

If you feel like trying the second option, go ahead. I can't think of a better solution yet. Leaving everything "as is" is fine too, I don't think fetching a unique_ptr is likely to have an observable effect, given the complexity of the function.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I implemented the second option to deal with the circular dependency, with KernelWrapper now a nested class of TableFunction. I removed the unique_ptr and the KernelWrapper object is now a direct member of TableFunction. Thanks!

@francoishamon
Copy link
Contributor Author

Thanks @klevzoff for looking at this PR so quickly! I am now going to add the TableRelativePermeability relperms to one of the integrated tests to make sure everything works on Lassen (so far I have only checked that it compiles there), and I will let you know when the code is ready for the final review.

@francoishamon francoishamon marked this pull request as ready for review January 26, 2021 05:28
Comment on lines 175 to 181
for( localIndex ip = 0; ip < NP; ++ip )
{
for( localIndex jp = 0; jp < NP; ++jp )
{
dPhaseRelPerm_dPhaseVolFrac[ip][jp] = 0.0;
}
}
Copy link
Contributor Author

Choose a reason for hiding this comment

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

@corbett5 To follow-up on this morning's question, what would be a good syntax to set the values in this arraySlice2d to zero? Here NP (soon changed to nPhases) is dynamic so I cannot use LvArray::tensorOps::fill. Is there a looping function in sliceHelpers.hpp that I could use here?

Copy link
Contributor

Choose a reason for hiding this comment

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

I would just recommend doing

for( double & val : dPhaseRelPerm_dPhaseVolFrac )
{ val = 0.0; }

This only works when dPhaseRelPerm_dPhaseVolFrac has continues values (otherwise it's a runtime error), but that appears to be the case here.

If you're operating on an Array or ArrayView there's a setValues method that takes a policy template parameter.

Copy link
Contributor

@TotoGaz TotoGaz left a comment

Choose a reason for hiding this comment

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

A few things here and there, nothing smart

Comment on lines +390 to +411
array1d< real64_array > coordinates;
coordinates.resize( Ndim );
coordinates[0].resize( Nx );
coordinates[0][0] = -1.0;
coordinates[0][1] = 0.0;
coordinates[0][2] = 1.0;
coordinates[1].resize( Ny );
coordinates[1][0] = -1.0;
coordinates[1][1] = 0.0;
coordinates[1][2] = 0.5;
coordinates[1][3] = 1.0;
coordinates[2].resize( Nz );
coordinates[2][0] = -1.0;
coordinates[2][1] = -0.4;
coordinates[2][2] = 0.3;
coordinates[2][3] = 0.5;
coordinates[2][4] = 1.0;
coordinates[3].resize( Nt );
coordinates[3][0] = -1.0;
coordinates[3][1] = 0.34;
coordinates[3][2] = 0.5;
coordinates[3][3] = 1.0;
Copy link
Contributor

Choose a reason for hiding this comment

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

Hmm, LvArray should allow std::initializer_list assignments!

@francoishamon
Copy link
Contributor Author

A few things here and there, nothing smart

Thanks @TotoGaz for looking at this PR

GEOSX_ERROR_IF( increment != m_values.size(), "Table dimensions do not match!" );
if( m_coordinates.size() > 0 && m_values.size() > 0 ) // coordinates and values have been set
{
GEOSX_ERROR_IF( increment != m_values.size(), "Table dimensions do not match!" );
Copy link
Contributor

Choose a reason for hiding this comment

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

Is there a check in the original/new version to make sure that the axes are monotonic (maybe I'm missing it)? That could certainly cause some unclear error messages...

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@cssherman Good idea, there was not check in TableFunction originally, but I added one in the relperm constitutive model that uses the tables. In the latest commit I moved it to TableFunction, so the check will be applied for all the tables, not all relperm tables. Thanks!

Copy link
Contributor

@klevzoff klevzoff left a comment

Choose a reason for hiding this comment

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

Looks good to me!

@francoishamon francoishamon added ci: run CUDA builds Allows to triggers (costly) CUDA jobs and removed flag: ready for review labels Feb 19, 2021
@francoishamon francoishamon merged commit b9fb6bc into develop Feb 19, 2021
@francoishamon francoishamon deleted the feature/hamon/relpermTable branch February 19, 2021 18:16
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
ci: run CUDA builds Allows to triggers (costly) CUDA jobs
Projects
None yet
Development

Successfully merging this pull request may close these issues.

[Feature] TableFunctions applied on device [Feature] Customized realtive permeability table
5 participants