Skip to content

Commit

Permalink
Several corrections and improvements for the settlement workflow
Browse files Browse the repository at this point in the history 
- Make sure that the GiD output writer's lifetime matches the time step loop, to avoid that output results are being overwritten after each step.
- Make sure that nodal variables are reset and restored as needed.
- Added a missing process that finds neighboring elements.
- Added a missing process that deactivates conditions on inactive elements.
- Extracted several convenience functions to make the code more readable.
- Reformatted some code using clang-format.
- When creating an output writer only pass the relevant settings rather than all settings.
- Moved finalizing the writing of results to a more appropriate location (i.e. at the end of the time step loop).
- Initialize the strategy wrapper at the start of the time step loop.
- The test for the settlement workflow now also runs the fourth stage. Accepted the results of the corresponding Python workflow as the expected results.
  • Loading branch information
@rfaasse
rfaasse authored Oct 28, 2023
1 parent c914524 commit 7988d63
Show file tree
Hide file tree
Showing 15 changed files with 57,769 additions and 164 deletions.
8 changes: 6 additions & 2 deletions applications/GeoMechanicsApplication/custom_workflows/dgeoflow.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -388,7 +388,9 @@ namespace Kratos
if (!hasPiping)
{
MainExecution(model_part, processes, p_solving_strategy, 0.0, 1.0, 1);
GeoOutputWriter::WriteGiDOutput(model_part, projectfile, rWorkingDirectory);
const auto gid_output_settings = projectfile["output_processes"]["gid_output"][0]["Parameters"];
GeoOutputWriter writer{gid_output_settings, rWorkingDirectory, model_part};
writer.WriteGiDOutput(model_part, gid_output_settings);
}
else
{
Expand Down Expand Up @@ -473,7 +475,9 @@ namespace Kratos
break;
}

GeoOutputWriter::WriteGiDOutput(model_part, projectfile, rWorkingDirectory);
const auto gid_output_settings = projectfile["output_processes"]["gid_output"][0]["Parameters"];
GeoOutputWriter writer{gid_output_settings, rWorkingDirectory, model_part};
writer.WriteGiDOutput(model_part, gid_output_settings);

// Update boundary conditions for next search head.
if (RiverBoundary->Info() == "ApplyConstantScalarValueProcess")
Expand Down
70 changes: 52 additions & 18 deletions applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,6 @@
#include "input_output/logger.h"
#include "time_loop_executor_interface.h"

#include "utilities/variable_utils.h"

#include "custom_processes/apply_scalar_constraint_table_process.h"
#include "custom_processes/apply_normal_load_table_process.h"
#include "custom_processes/apply_vector_constraint_table_process.h"
Expand All @@ -29,6 +27,8 @@
#include "spaces/ublas_space.h"
#include "solving_strategy_wrapper.hpp"
#include "adaptive_time_incrementor.h"
#include "custom_processes/find_neighbour_elements_of_conditions_process.hpp"
#include "custom_processes/deactivate_conditions_on_inactive_elements_process.hpp"

namespace
{
Expand Down Expand Up @@ -88,7 +88,6 @@ bool GetResetDisplacementsFrom(const Parameters& rProjectParameters)
}



namespace Kratos
{

Expand All @@ -99,7 +98,7 @@ KratosGeoSettlement::KratosGeoSettlement(std::unique_ptr<InputUtility> pInputUti
mpProcessInfoParser{std::move(pProcessInfoParser)},
mpTimeLoopExecutor{std::move(pTimeLoopExecutorInterface)}
{
mKernel.GetApplicationsList().clear();
Kernel::GetApplicationsList().clear();
mModel.Reset();
KRATOS_INFO("KratosGeoSettlement") << "Setting up Kratos" << std::endl;
KRATOS_ERROR_IF_NOT(mpInputUtility) << "Invalid Input Utility";
Expand Down Expand Up @@ -135,7 +134,7 @@ void KratosGeoSettlement::InitializeProcessFactory()
mProcessFactory->AddCreator("ApplyNormalLoadTableProcess",
[this](const Parameters& rParameters)
{
return std::make_unique<ApplyNormalLoadTableProcess>(mModel.GetModelPart(mModelPartName),
return std::make_unique<ApplyNormalLoadTableProcess>(GetMainModelPart(),
rParameters);
});

Expand Down Expand Up @@ -216,6 +215,7 @@ int KratosGeoSettlement::RunStage(const std::filesystem::path& rWorki

std::vector<std::shared_ptr<Process>> processes = GetProcesses(project_parameters);
std::vector<std::weak_ptr<Process>> process_observables(processes.begin(), processes.end());

for (const auto& process : processes) {
process->ExecuteInitialize();
}
Expand All @@ -231,9 +231,14 @@ int KratosGeoSettlement::RunStage(const std::filesystem::path& rWorki
mpTimeLoopExecutor->SetSolverStrategyWrapper(MakeStrategyWrapper(project_parameters,
rWorkingDirectory));
// For now, pass a dummy state. THIS PROBABLY NEEDS TO BE REFINED AT SOME POINT!
TimeStepEndState dummy_state;
dummy_state.convergence_state = TimeStepEndState::ConvergenceState::converged;
mpTimeLoopExecutor->Run(dummy_state);
TimeStepEndState start_of_loop_state;
start_of_loop_state.convergence_state = TimeStepEndState::ConvergenceState::converged;
start_of_loop_state.time = GetStartTimeFrom(project_parameters);
mpTimeLoopExecutor->Run(start_of_loop_state);
}

for (const auto& process : processes) {
process->ExecuteFinalize();
}

FlushLoggingOutput(rLogCallback, logger_output, kratos_log_buffer);
Expand Down Expand Up @@ -348,11 +353,25 @@ std::unique_ptr<TimeIncrementor> KratosGeoSettlement::MakeTimeIncrementor(const
std::shared_ptr<StrategyWrapper> KratosGeoSettlement::MakeStrategyWrapper(const Parameters& rProjectParameters,
const std::filesystem::path& rWorkingDirectory)
{
auto& main_model_part = mModel.GetModelPart(mModelPartName);
auto solving_strategy = SolvingStrategyFactoryType::Create(rProjectParameters["solver_settings"],
main_model_part.GetSubModelPart(mComputationalSubModelPartName));
GetComputationalModelPart());
KRATOS_ERROR_IF_NOT(solving_strategy) << "No solving strategy was created!" << std::endl;

solving_strategy->Initialize();
GetMainModelPart().CloneTimeStep();

if (rProjectParameters["solver_settings"]["reset_displacements"].GetBool()) {
constexpr auto source_index = std::size_t{0};
constexpr auto destination_index = std::size_t{1};
RestoreValuesOfNodalVariable(DISPLACEMENT, source_index, destination_index);
RestoreValuesOfNodalVariable(ROTATION, source_index, destination_index);

VariableUtils{}.UpdateCurrentToInitialConfiguration(GetComputationalModelPart().Nodes());
}

FindNeighbourElementsOfConditionsProcess{GetComputationalModelPart()}.Execute();
DeactivateConditionsOnInactiveElements{GetComputationalModelPart()}.Execute();

// For now, we can create solving strategy wrappers only
using SolvingStrategyWrapperType = SolvingStrategyWrapper<SparseSpaceType, DenseSpaceType>;
return std::make_shared<SolvingStrategyWrapperType>(std::move(solving_strategy),
Expand All @@ -363,23 +382,28 @@ std::shared_ptr<StrategyWrapper> KratosGeoSettlement::MakeStrategyWrapper(const

void KratosGeoSettlement::PrepareModelPart(const Parameters& rSolverSettings)
{
auto& main_model_part = mModel.GetModelPart(mModelPartName);
auto& main_model_part = GetMainModelPart();
if (!main_model_part.HasSubModelPart(mComputationalSubModelPartName)) {
main_model_part.CreateSubModelPart(mComputationalSubModelPartName);
}
auto& computing_model_part = main_model_part.GetSubModelPart(mComputationalSubModelPartName);

if (rSolverSettings.Has("nodal_smoothing"))
{
main_model_part.GetProcessInfo().SetValue(NODAL_SMOOTHING, rSolverSettings["nodal_smoothing"].GetBool());
}

// Note that the computing part and the main model part _share_ their process info and properties
computing_model_part.SetProcessInfo(main_model_part.GetProcessInfo());
GetComputationalModelPart().SetProcessInfo(main_model_part.GetProcessInfo());
for (auto i = ModelPart::SizeType{0}; i < main_model_part.NumberOfMeshes(); ++i)
{
auto& mesh = main_model_part.GetMesh(i);
for (const auto& property : mesh.Properties())
{
computing_model_part.AddProperties( mesh.pGetProperties(property.GetId()));
GetComputationalModelPart().AddProperties( mesh.pGetProperties(property.GetId()));
}
}

computing_model_part.Set(ACTIVE);
GetComputationalModelPart().Set(ACTIVE);

const auto problem_domain_sub_model_part_list = rSolverSettings["problem_domain_sub_model_part_list"];
std::vector<std::string> domain_part_names;
Expand All @@ -395,7 +419,7 @@ void KratosGeoSettlement::PrepareModelPart(const Parameters& rSolverSettings)
node_id_set.insert(node.Id());
}
}
computing_model_part.AddNodes(std::vector<Node::IndexType>{node_id_set.begin(), node_id_set.end()});
GetComputationalModelPart().AddNodes(std::vector<Node::IndexType>{node_id_set.begin(), node_id_set.end()});

std::set<IndexedObject::IndexType> element_id_set;
for (const auto& name : domain_part_names) {
Expand All @@ -404,7 +428,7 @@ void KratosGeoSettlement::PrepareModelPart(const Parameters& rSolverSettings)
element_id_set.insert(element.Id());
}
}
computing_model_part.AddElements(std::vector<IndexedObject::IndexType>{element_id_set.begin(), element_id_set.end()});
GetComputationalModelPart().AddElements(std::vector<IndexedObject::IndexType>{element_id_set.begin(), element_id_set.end()});

const auto processes_sub_model_part_list = rSolverSettings["processes_sub_model_part_list"];
std::vector<std::string> domain_condition_names;
Expand All @@ -419,12 +443,22 @@ void KratosGeoSettlement::PrepareModelPart(const Parameters& rSolverSettings)
condition_id_set.insert(condition.Id());
}
}
computing_model_part.AddConditions(std::vector<IndexedObject::IndexType>{condition_id_set.begin(), condition_id_set.end()});
GetComputationalModelPart().AddConditions(std::vector<IndexedObject::IndexType>{condition_id_set.begin(), condition_id_set.end()});

// NOTE TO SELF: here we don't yet "Adding Computing Sub Sub Model Parts" (see check_and_prepare_model_process_geo.py)
// We are not sure yet whether that piece of code is relevant or not.
}

ModelPart& KratosGeoSettlement::GetMainModelPart()
{
return mModel.GetModelPart(mModelPartName);
}

ModelPart& KratosGeoSettlement::GetComputationalModelPart()
{
return GetMainModelPart().GetSubModelPart(mComputationalSubModelPartName);
}

// This default destructor is added in the cpp to be able to forward member variables in a unique_ptr
KratosGeoSettlement::~KratosGeoSettlement() = default;

Expand Down
20 changes: 20 additions & 0 deletions applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@

#include "geo_mechanics_application.h"
#include "linear_solvers_application.h"
#include "utilities/variable_utils.h"

namespace Kratos
{
Expand Down Expand Up @@ -55,6 +56,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) KratosGeoSettlement
ModelPart& AddNewModelPart(const std::string& rModelPartName);
void PrepareModelPart(const Parameters& rSolverSettings);

ModelPart& GetMainModelPart();
ModelPart& GetComputationalModelPart();

static void AddNodalSolutionStepVariablesTo(ModelPart& rModelPart);
static void AddDegreesOfFreedomTo(ModelPart& rModelPart);
void InitializeProcessFactory();
Expand All @@ -65,6 +69,22 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) KratosGeoSettlement
LoggerOutput::Pointer CreateLoggingOutput(std::stringstream& rKratosLogBuffer) const;
void FlushLoggingOutput(const std::function<void(const char*)>& rLogCallback, LoggerOutput::Pointer pLoggerOutput, const std::stringstream& rKratosLogBuffer) const;

template <typename TVariableType>
void RestoreValuesOfNodalVariable(const TVariableType& rVariable,
Node::IndexType SourceIndex,
Node::IndexType DestinationIndex)
{
if (!GetComputationalModelPart().HasNodalSolutionStepVariable(rVariable))
return;

VariableUtils{}.SetHistoricalVariableToZero(rVariable, GetComputationalModelPart().Nodes());

block_for_each(GetComputationalModelPart().Nodes(),
[&rVariable, SourceIndex, DestinationIndex](auto& node) {
node.GetSolutionStepValue(rVariable, DestinationIndex) = node.GetSolutionStepValue(rVariable, SourceIndex);
});
}

Kernel mKernel;
Model mModel;
std::string mModelPartName;
Expand Down
Loading

0 comments on commit 7988d63

Please sign in to comment.