Skip to content

Commit

Permalink
Merge branch 'master' into core/mpi/specialized-spatial-search-mpi
Browse files Browse the repository at this point in the history
  • Loading branch information
@loumalouomega
loumalouomega committed Jun 15, 2023
2 parents 202fbff + 94b7b41 commit 4300bac
Show file tree
Hide file tree
Showing 8 changed files with 260 additions and 114 deletions.
87 changes: 25 additions & 62 deletions applications/ParticleMechanicsApplication/custom_utilities/mpm_boundary_rotation_utility.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,19 +157,19 @@ class MPMBoundaryRotationUtility: public CoordinateTransformationUtils<TLocalMat
array_1d<double,3> rN = rGeometry[itNode].FastGetSolutionStepValue(NORMAL);
this->Normalize(rN);

for( unsigned int i = 0; i < j; ++i)// Skip term (i,i)
{
rLocalMatrix(i,j) = 0.0;
rLocalMatrix(j,i) = 0.0;
}
for( unsigned int i = j+1; i < LocalSize; ++i)
{
rLocalMatrix(i,j) = 0.0;
rLocalMatrix(j,i) = 0.0;
}

// Zero out row/column corresponding to normal displacement DoF except diagonal term (set to 1)
// Applied IFF the local matrix passed is not empty [otherwise does nothing -- RHS only case]
if (rLocalMatrix.size1() != 0) {
for( unsigned int i = 0; i < LocalSize; ++i)
{
rLocalMatrix(i,j) = 0.0;
rLocalMatrix(j,i) = 0.0;
}
rLocalMatrix(j, j) = 1.0; // set diagonal term to 1.0
}

// Set value of normal displacement at node directly to the normal displacement of the boundary mesh
rLocalVector[j] = inner_prod(rN,displacement);
rLocalMatrix(j,j) = 1.0;
}
}
}
Expand All @@ -179,26 +179,10 @@ class MPMBoundaryRotationUtility: public CoordinateTransformationUtils<TLocalMat
void ApplySlipCondition(TLocalVectorType& rLocalVector,
GeometryType& rGeometry) const override
{
if (rLocalVector.size() > 0)
{
for(unsigned int itNode = 0; itNode < rGeometry.PointsNumber(); ++itNode)
{
if( this->IsSlip(rGeometry[itNode]) )
{
// We fix the first momentum dof (normal component) for each rotated block
unsigned int j = itNode * this->GetBlockSize(); // +1

// Get the displacement of the boundary mesh, this does not assume that the mesh is moving.
// If the mesh is moving, need to consider the displacement of the moving mesh into account.
const array_1d<double,3> & displacement = rGeometry[itNode].FastGetSolutionStepValue(DISPLACEMENT);
array_1d<double,3> rN = rGeometry[itNode].FastGetSolutionStepValue(NORMAL);
this->Normalize(rN);

rLocalVector[j] = inner_prod(rN,displacement);

}
}
}
// creates an empty dummy matrix to pass into the 'full' ApplySlipCondition -- this dummy matrix is
// ignored, effectively only updating the RHS
TLocalMatrixType dummyMatrix;
this->ApplySlipCondition(dummyMatrix, rLocalVector, rGeometry);
}

// An extra function to distinguish the application of slip in element considering penalty imposition
Expand All @@ -207,7 +191,7 @@ class MPMBoundaryRotationUtility: public CoordinateTransformationUtils<TLocalMat
GeometryType& rGeometry) const
{
// If it is not a penalty element, do as standard
// Otherwise, if it is a penalty element, dont do anything
// Otherwise, if it is a penalty element, don't do anything
if (!this->IsPenalty(rGeometry))
{
this->ApplySlipCondition(rLocalMatrix, rLocalVector, rGeometry);
Expand All @@ -219,7 +203,7 @@ class MPMBoundaryRotationUtility: public CoordinateTransformationUtils<TLocalMat
GeometryType& rGeometry) const
{
// If it is not a penalty element, do as standard
// Otherwise, if it is a penalty element, dont do anything
// Otherwise, if it is a penalty element, don't do anything
if (!this->IsPenalty(rGeometry))
{
this->ApplySlipCondition(rLocalVector, rGeometry);
Expand Down Expand Up @@ -253,6 +237,7 @@ class MPMBoundaryRotationUtility: public CoordinateTransformationUtils<TLocalMat
unsigned int j = itNode * block_size;

// Copy all normal value in LHS to the temp_matrix
// [ does nothing for dummy rLocalMatrix (size1() == 0) -- RHS only case ]
for (unsigned int i = j; i < rLocalMatrix.size1(); i+= block_size)
{
temp_matrix(i,j) = rLocalMatrix(i,j);
Expand All @@ -266,6 +251,8 @@ class MPMBoundaryRotationUtility: public CoordinateTransformationUtils<TLocalMat
}
}
}
// All entries in penalty matrix zeroed out except for normal component
// [ no effect in case of empty dummy rLocalMatrix ]
rLocalMatrix = temp_matrix;
}
}
Expand All @@ -275,34 +262,10 @@ class MPMBoundaryRotationUtility: public CoordinateTransformationUtils<TLocalMat
void ConditionApplySlipCondition(TLocalVectorType& rLocalVector,
GeometryType& rGeometry) const
{
// If it is not a penalty condition, do as standard
if (!this->IsPenalty(rGeometry))
{
this->ApplySlipCondition(rLocalVector, rGeometry);
}
// Otherwise, if it is a penalty element, dont do anything
else
{
if (rLocalVector.size() > 0)
{
const unsigned int block_size = this->GetBlockSize();
for(unsigned int itNode = 0; itNode < rGeometry.PointsNumber(); ++itNode)
{
if( this->IsSlip(rGeometry[itNode]) )
{
// We fix the first momentum dof (normal component) for each rotated block
unsigned int j = itNode * block_size;

// Remove all other value than the normal component
for(unsigned int i = j; i < (j + block_size); ++i)
{
if (i!=j) rLocalVector[i] = 0.0;
}
}
}
}
}

// creates an empty dummy matrix to pass into the 'full' ConditionApplySlipCondition -- this dummy matrix is
// ignored, effectively only updating the RHS
TLocalMatrixType dummyMatrix;
this->ConditionApplySlipCondition(dummyMatrix, rLocalVector, rGeometry);
}

// Checking whether it is normal element or penalty element
Expand Down
1 change: 1 addition & 0 deletions applications/TrilinosApplication/custom_python/add_trilinos_strategies_to_python.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,7 @@ void AddStrategies(pybind11::module& m)
typedef TrilinosBlockBuilderAndSolver< TrilinosSparseSpaceType, TrilinosLocalSpaceType, TrilinosLinearSolverType > TrilinosBlockBuilderAndSolverType;
py::class_<TrilinosBlockBuilderAndSolverType, typename TrilinosBlockBuilderAndSolverType::Pointer, TrilinosBuilderAndSolverType>(m, "TrilinosBlockBuilderAndSolver")
.def(py::init<Epetra_MpiComm&, int, TrilinosLinearSolverType::Pointer > () )
.def(py::init<Epetra_MpiComm&, TrilinosLinearSolverType::Pointer, Parameters > () )
;

typedef TrilinosBlockBuilderAndSolverPeriodic< TrilinosSparseSpaceType, TrilinosLocalSpaceType, TrilinosLinearSolverType > TrilinosBlockBuilderAndSolverPeriodicType;
Expand Down
92 changes: 89 additions & 3 deletions ...inosApplication/custom_strategies/builder_and_solvers/trilinos_block_builder_and_solver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,9 @@ class TrilinosBlockBuilderAndSolver
///@name Type Definitions
///@{

/// Definition of the flags
KRATOS_DEFINE_LOCAL_FLAG( SILENT_WARNINGS );

/// Definition of the pointer
KRATOS_CLASS_POINTER_DEFINITION(TrilinosBlockBuilderAndSolver);

Expand Down Expand Up @@ -140,9 +143,19 @@ class TrilinosBlockBuilderAndSolver
}

/**
* @brief Default destructor.
* @brief Default constructor. (with parameters)
*/
~TrilinosBlockBuilderAndSolver() override = default;
explicit TrilinosBlockBuilderAndSolver(
EpetraCommunicatorType& rComm,
typename TLinearSolver::Pointer pNewLinearSystemSolver,
Parameters ThisParameters
) : BaseType(pNewLinearSystemSolver),
mrComm(rComm)
{
// Validate and assign defaults
ThisParameters = this->ValidateAndAssignParameters(ThisParameters, this->GetDefaultParameters());
this->AssignSettings(ThisParameters);
}

/**
* Copy constructor
Expand All @@ -154,6 +167,8 @@ class TrilinosBlockBuilderAndSolver
*/
TrilinosBlockBuilderAndSolver& operator=(const TrilinosBlockBuilderAndSolver& rOther) = delete;

// TODO: In order to create a Create method, the name of the DataCommunicator that is used needs to be passed in the settings (see DistributedImportModelPartUtility). Then an EpetraComm can be constructed from the MPI_Comm in the DataCommunicator

///@}
///@name Operators
///@{
Expand Down Expand Up @@ -831,6 +846,10 @@ class TrilinosBlockBuilderAndSolver
// defining a temporary vector to gather all of the values needed
Epetra_IntVector fixed(rA.ColMap());

// Detect if there is a line of all zeros and set the diagonal to a 1 if this happens
const auto& r_process_info = rModelPart.GetProcessInfo();
mScaleFactor = TSparseSpace::CheckAndCorrectZeroDiagonalValues(r_process_info, rA, rb, mScalingDiagonal);

// Importing in the new temp vector the values
int ierr = fixed.Import(fixed_local, dirichlet_importer, Insert);
KRATOS_ERROR_IF(ierr != 0) << "Epetra failure found" << std::endl;
Expand Down Expand Up @@ -889,6 +908,27 @@ class TrilinosBlockBuilderAndSolver
KRATOS_CATCH("");
}

/**
* @brief This method provides the defaults parameters to avoid conflicts between the different constructors
* @return The default parameters
*/
Parameters GetDefaultParameters() const override
{
Parameters default_parameters = Parameters(R"(
{
"name" : "trilinos_block_builder_and_solver",
"guess_row_size" : 45,
"block_builder" : true,
"diagonal_values_for_dirichlet_dofs" : "use_max_diagonal",
"silent_warnings" : false
})");

// Getting base class default parameters
const Parameters base_default_parameters = BaseType::GetDefaultParameters();
default_parameters.RecursivelyAddMissingParameters(base_default_parameters);
return default_parameters;
}

/**
* @brief Returns the name of the class as used in the settings (snake_case format)
* @return The name of the class
Expand Down Expand Up @@ -946,7 +986,13 @@ class TrilinosBlockBuilderAndSolver
int mGuessRowSize; /// The guess row size
IndexType mLocalSystemSize; /// The local system size
int mFirstMyId; /// Auxiliary Id (I)
int mLastMyId; /// Auxiliary Id (II) /// Some flags used internally
int mLastMyId; /// Auxiliary Id (II)

double mScaleFactor = 1.0; /// The manually set scale factor

/* Flags */
SCALING_DIAGONAL mScalingDiagonal = SCALING_DIAGONAL::CONSIDER_MAX_DIAGONAL; /// We identify the scaling considered for the dirichlet dofs
Flags mOptions; /// Some flags used internally

///@}
///@name Protected Operators
Expand All @@ -956,6 +1002,42 @@ class TrilinosBlockBuilderAndSolver
///@name Protected Operations
///@{

/**
* @brief This method assigns settings to member variables
* @param ThisParameters Parameters that are assigned to the member variables
*/
void AssignSettings(const Parameters ThisParameters) override
{
BaseType::AssignSettings(ThisParameters);

// Get guess row size
mGuessRowSize = ThisParameters["guess_row_size"].GetInt();

// Setting flags
const std::string& r_diagonal_values_for_dirichlet_dofs = ThisParameters["diagonal_values_for_dirichlet_dofs"].GetString();

const std::set<std::string> available_options_for_diagonal = {"no_scaling","use_max_diagonal","use_diagonal_norm","defined_in_process_info"};

if (available_options_for_diagonal.find(r_diagonal_values_for_dirichlet_dofs) == available_options_for_diagonal.end()) {
std::stringstream msg;
msg << "Currently prescribed diagonal values for dirichlet dofs : " << r_diagonal_values_for_dirichlet_dofs << "\n";
msg << "Admissible values for the diagonal scaling are : 'no_scaling', 'use_max_diagonal', 'use_diagonal_norm', or 'defined_in_process_info'" << "\n";
KRATOS_ERROR << msg.str() << std::endl;
}

// The first option will not consider any scaling (the diagonal values will be replaced with 1)
if (r_diagonal_values_for_dirichlet_dofs == "no_scaling") {
mScalingDiagonal = SCALING_DIAGONAL::NO_SCALING;
} else if (r_diagonal_values_for_dirichlet_dofs == "use_max_diagonal") {
mScalingDiagonal = SCALING_DIAGONAL::CONSIDER_MAX_DIAGONAL;
} else if (r_diagonal_values_for_dirichlet_dofs == "use_diagonal_norm") { // On this case the norm of the diagonal will be considered
mScalingDiagonal = SCALING_DIAGONAL::CONSIDER_NORM_DIAGONAL;
} else { // Otherwise we will assume we impose a numerical value
mScalingDiagonal = SCALING_DIAGONAL::CONSIDER_PRESCRIBED_DIAGONAL;
}
mOptions.Set(SILENT_WARNINGS, ThisParameters["silent_warnings"].GetBool());
}

///@}
///@name Protected Access
///@{
Expand Down Expand Up @@ -1012,6 +1094,10 @@ class TrilinosBlockBuilderAndSolver
///@name Type Definitions
///@{

// Here one should use the KRATOS_CREATE_LOCAL_FLAG, but it does not play nice with template parameters
template<class TSparseSpace, class TDenseSpace, class TLinearSolver>
const Kratos::Flags TrilinosBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver>::SILENT_WARNINGS(Kratos::Flags::Create(0));

///@}

} /* namespace Kratos.*/
Loading

0 comments on commit 4300bac

Please sign in to comment.