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

feat: Introduce full tensor permeability #2822

Open
wants to merge 11 commits into
base: develop
Choose a base branch
from
Open

feat: Introduce full tensor permeability #2822

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
b8f4c46
boundary stencil
CusiniM Nov 15, 2023
ea3fd3b
Merge branch 'develop' into feature/cus-cas/fullTensorPerm
CusiniM Nov 15, 2023
638a410
missing ;
CusiniM Nov 15, 2023
b740ec4
it compiles
CusiniM Nov 17, 2023
eb1fa9a
Full tensor permebility test
castelletto1 Nov 27, 2023
ce54941
Merge remote-tracking branch 'origin/develop' into feature/cus-cas/fu…
CusiniM Mar 13, 2024
5274772
Merge branch 'feature/cus-cas/fullTensorPerm' of github.com:GEOS-DEV/…
CusiniM Mar 13, 2024
5a10e64
Merge branch 'refs/heads/develop' into feature/cus-cas/fullTensorPerm
jafranc May 23, 2024
a9314e2
Merge branch 'develop' into feature/cus-cas/fullTensorPerm
rrsettgast Dec 17, 2024
e4c929a
Post merge fixes
castelletto1 Dec 20, 2024
cc53d0d
Post merge fixes
castelletto1 Dec 20, 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
4 changes: 3 additions & 1 deletion src/coreComponents/constitutive/permeability/CarmanKozenyPermeability.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,8 @@ CarmanKozenyPermeability::CarmanKozenyPermeability( string const & name, Group *
setDefaultValue( m_anisotropy ).
setDescription( "Anisotropy factors for three permeability components." );

registerWrapper( viewKeyStruct::dPerm_dPorosityString(), &m_dPerm_dPorosity );
registerWrapper( viewKeyStruct::dPerm_dPorosityString(), &m_dPerm_dPorosity ).
setApplyDefaultValue(0.0);
}

std::unique_ptr< ConstitutiveBase >
Expand All @@ -60,6 +61,7 @@ void CarmanKozenyPermeability::allocateConstitutiveData( dataRepository::Group &
// NOTE: enforcing 1 quadrature point
m_dPerm_dPorosity.resize( 0, 1, 3 );
PermeabilityBase::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex );
getWrapper( PermeabilityBase::viewKeyStruct::permeabilityString() ).setApplyDefaultValue(0.0);
}

REGISTER_CATALOG_ENTRY( ConstitutiveBase, CarmanKozenyPermeability, string const &, Group * const )
Expand Down
2 changes: 1 addition & 1 deletion src/coreComponents/constitutive/permeability/CarmanKozenyPermeability.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -145,7 +145,7 @@ void CarmanKozenyPermeabilityUpdate::compute( real64 const & porosity,

real64 const dPerm_dPorValue = -constant * ( (porosity - 3) * pow( porosity, 2 ) / pow( (1-porosity), 3 ) );

for( localIndex i = 0; i < permeability.size(); ++i )
for( localIndex i = 0; i < m_anisotropy[i].size(); ++i )
{
permeability[i] = permValue * m_anisotropy[i];
dPerm_dPorosity[i] = dPerm_dPorValue * m_anisotropy[i];
Expand Down
53 changes: 40 additions & 13 deletions src/coreComponents/constitutive/permeability/ConstantPermeability.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,10 +30,38 @@ namespace constitutive
ConstantPermeability::ConstantPermeability( string const & name, Group * const parent ):
PermeabilityBase( name, parent )
{
registerWrapper( viewKeyStruct::permeabilityComponentsString(), &m_permeabilityComponents ).
setInputFlag( InputFlags::REQUIRED ).
registerWrapper( viewKeyStruct::diagonalPermeabilityTensorString(), &m_diagonalPermeabilityTensor ).
setInputFlag( InputFlags::OPTIONAL ).
setRestartFlags( RestartFlags::NO_WRITE ).
setApplyDefaultValue(-1.0)
setDescription( "xx, yy and zz components of a diagonal permeability tensor." );

registerWrapper( viewKeyStruct::symmetricFullPermabilityTensorString(), &m_symmetricFullPermeabilityTensor ).
setInputFlag( InputFlags::OPTIONAL ).
setRestartFlags( RestartFlags::NO_WRITE ).
setApplyDefaultValue(-1.0).
setDescription( "xx, yy and zz, xz, yz, xy components of a symmetric permeability tensor." );
}

ConstantPermeability::postProcessInput()
{
GEOS_ERROR_IF( m_diagonalPermeabilityTensor[0] < 0.0 && m_symmetricFullPermeabilityTensor[0] < 0.0,
"Either a diagonal permeability tensor or a full tensor must be provided.");
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Suggested change
"Either a diagonal permeability tensor or a full tensor must be provided.");
"Either a diagonal permeability tensor or a symmetric full tensor must be provided.");


GEOS_ERROR_IF( m_diagonalPermeabilityTensor[0] > 0.0 && m_symmetricFullPermeabilityTensor[0] > 0.0,
"Only one between a diagonal permeability tensor and a full tensor permeability can be provided.");

if ( m_diagonalPermeabilityTensor[0] > 0.0 )
{
for( int i = 0; i < 3; i++ )
{
m_symmetricFullPermeabilityTensor[i] = m_diagonalPermeabilityTensor[i];
}
for( int j=4; j < 6; i++ )
{
m_symmetricFullPermeabilityTensor[i] = 0.0;
}
}
}

std::unique_ptr< ConstitutiveBase >
Expand All @@ -54,9 +82,8 @@ void ConstantPermeability::allocateConstitutiveData( dataRepository::Group & par
{
for( localIndex q = 0; q < numQuad; ++q )
{
m_permeability[ei][q][0] = m_permeabilityComponents[0];
m_permeability[ei][q][1] = m_permeabilityComponents[1];
m_permeability[ei][q][2] = m_permeabilityComponents[2];
for ( int c=0; c < 6; c++)
m_permeability[ei][q][c] = m_symmetricFullPermeabilityTensor[c];
Copy link
Contributor

Choose a reason for hiding this comment

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

Shouldn't we enforce positiveness of diagonal component here ?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

We should really throw an error earlier if someone inputs a negative value but, for now, I am using the default negative value to distinguish whether a diagonal or a full tensor was provided so I need to improve that logic fist.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

One option would be to simply provide 2 R1Tensors one for the main diagonal and one for the off diagonal terms and have the off diagonal components default to 0.0 if not provided.

Copy link
Contributor

Choose a reason for hiding this comment

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

Just noticed that the order would be imposed by LvArray::Ri_eq_symAijBj implicitly :/

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

yes, and it's different from the one used by vtk.

}
}
}
Expand All @@ -67,18 +94,21 @@ void ConstantPermeability::initializeState() const
integer constexpr numQuad = 1; // NOTE: enforcing 1 quadrature point

auto permView = m_permeability.toView();
real64 const permComponents[3] = { m_permeabilityComponents[0],
m_permeabilityComponents[1],
m_permeabilityComponents[2] };
real64 const permComponents[3] = { m_symmetricFullPermeabilityTensor[0],
m_symmetricFullPermeabilityTensor[1],
m_symmetricFullPermeabilityTensor[2],
m_symmetricFullPermeabilityTensor[3],
m_symmetricFullPermeabilityTensor[4],
m_symmetricFullPermeabilityTensor[5] };

forAll< parallelDevicePolicy<> >( numE, [=] GEOS_HOST_DEVICE ( localIndex const ei )
{
for( localIndex q = 0; q < numQuad; ++q )
{
for( integer dim=0; dim < 3; ++dim )
for( integer dim=0; dim < 6; ++dim )
{
// The default value is -1 so if it still -1 it needs to be set to something physical
if( permView[ei][q][0] < 0 )
if( permView[ei][q][dim] < 0 )
{
permView[ei][q][dim] = permComponents[dim];
}
Expand All @@ -87,9 +117,6 @@ void ConstantPermeability::initializeState() const
} );
}

void ConstantPermeability::postProcessInput()
{}

REGISTER_CATALOG_ENTRY( ConstitutiveBase, ConstantPermeability, string const &, Group * const )

}
Expand Down
8 changes: 6 additions & 2 deletions src/coreComponents/constitutive/permeability/ConstantPermeability.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,9 @@ class ConstantPermeability : public PermeabilityBase

struct viewKeyStruct : public PermeabilityBase::viewKeyStruct
{
static constexpr char const * permeabilityComponentsString() { return "permeabilityComponents"; }
static constexpr char const * diagonalPermeabilityTensorString() { return "diagonalPermeabilityTensor"; }
static constexpr char const * symmetricFullPermeabilityTensorString() { return "symmetricFullPermeabilityTensor"; }

} viewKeys;

virtual void initializeState() const override final;
Expand All @@ -84,7 +86,9 @@ class ConstantPermeability : public PermeabilityBase

private:

R1Tensor m_permeabilityComponents;
R1Tensor m_diagonalPermeabilityTensor;

R2SymTensor m_symmetricFullPermeabilityTensor;

};

Expand Down
4 changes: 2 additions & 2 deletions src/coreComponents/constitutive/permeability/PermeabilityBase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,8 +62,8 @@ void PermeabilityBase::allocateConstitutiveData( dataRepository::Group & parent,
localIndex const numConstitutivePointsPerParentIndex )
{
// NOTE: enforcing 1 quadrature point
m_permeability.resize( 0, 1, 3 );
m_dPerm_dPressure.resize( 0, 1, 3 );
m_permeability.resize( 0, 1, 6 );
m_dPerm_dPressure.resize( 0, 1, 6 );

ConstitutiveBase::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex );
}
Expand Down
3 changes: 2 additions & 1 deletion src/coreComponents/finiteVolume/BoundaryStencil.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -204,8 +204,9 @@ BoundaryStencilWrapper::
auto const compute = [&]
{
real64 faceConormal[3];
LvArray::tensorOps::hadamardProduct< 3 >( faceConormal, coef, faceNormal );
LvArray::tensorOps::Ri_eq_symAijBj< 3 >( faceConormal, coef, faceNormal )
CusiniM marked this conversation as resolved.
Show resolved Hide resolved
weight = LvArray::tensorOps::AiBi< 3 >( cellToFace, faceConormal );
// Matteo: this one also has to change!
LvArray::tensorOps::hadamardProduct< 3 >( dWeight_dCoef, cellToFace, faceNormal );
Copy link
Contributor

Choose a reason for hiding this comment

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

I am not sure about how this one should change.

};

Expand Down
8 changes: 3 additions & 5 deletions src/coreComponents/finiteVolume/CellElementStencilTPFA.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -242,23 +242,21 @@ CellElementStencilTPFAWrapper::
}

real64 faceConormal[3];
LvArray::tensorOps::hadamardProduct< 3 >( faceConormal, coefficient[er][esr][ei][0], faceNormal );
LvArray::tensorOps::Ri_eq_symAijBj< 3 >( faceConormal, coefficient[er][esr][ei][0], faceNormal )
halfWeight[i] *= LvArray::tensorOps::AiBi< 3 >( m_cellToFaceVec[iconn][i], faceConormal );

// correct negative weight issue arising from non-K-orthogonal grids
if( halfWeight[i] < 0.0 )
{
LvArray::tensorOps::hadamardProduct< 3 >( faceConormal,
coefficient[er][esr][ei][0],
m_cellToFaceVec[iconn][i] );
LvArray::tensorOps::Ri_eq_symAijBj< 3 >( faceConormal, coefficient[er][esr][ei][0], m_cellToFaceVec[iconn][i] )
halfWeight[i] = m_weights[iconn][i];
halfWeight[i] *= LvArray::tensorOps::AiBi< 3 >( m_cellToFaceVec[iconn][i], faceConormal );
}
}

// Do harmonic and arithmetic averaging
real64 const product = halfWeight[0]*halfWeight[1];
real64 const sum = halfWeight[0]+halfWeight[1];
real64 const sum = halfWeight[0]+halfWeight[1];

real64 const harmonicWeight = sum > 0 ? product / sum : 0.0;
real64 const arithmeticWeight = sum / 2;
Expand Down
Loading