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

Staggered pf #2052

Merged
merged 19 commits into from Mar 7, 2018

Conversation

Projects
None yet
3 participants
@KeitaYoshioka
Copy link
Contributor

KeitaYoshioka commented Jan 22, 2018

Staggered phase-field first implementation. Two test cases have been added to test.

@@ -57,6 +57,7 @@ struct IntegrationPointData final
double integration_weight;
double history_variable;
double history_variable_prev;
double pressure = 0.1;

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

0.1 is a magic number; please add a comment.

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

It's been changed to 1.0 and expiations why have been added.

void PhaseFieldProcess<DisplacementDim>::postNonLinearSolverConcreteProcess(
GlobalVector const& x, const double t, const int process_id)
{
if (!_use_monolithic_scheme && process_id == 0)

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

@wenqing Can this decision be made before calling this function?

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

postNonLinearSolver function is yet to be used in this commit. Thus, the function has been removed.

private:
PhaseFieldProcessData<DisplacementDim> _process_data;

std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;

std::unique_ptr<NumLib::LocalToGlobalIndexMap>
_local_to_global_index_map_single_component;

/// Sparsity pattern for the phase field equation, and it is initialized
// only if the staggered scheme is used.

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

slash missing

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

fixed.

REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
DIFF_DATA
expected_2D_StaticCrack_pcs_1_ts_1_t_1.000000.vtu 2D_StaticCrack_pcs_1_ts_1_t_1.000000.vtu displacement displacement 1e-6 1e-6
expected_2D_StaticCrack_pcs_1_ts_1_t_1.000000.vtu 2D_StaticCrack_pcs_1_ts_1_t_1.000000.vtu phasefield phasefield 1e-6 1e-6

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

Is tighter tolerance possible? Something around 1e-12?
The first tolerance is the absolute one, the second is relative, are both applicable to this test?

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

They have been made tighter (1e-15).

<name>linear_solver_u</name>
<eigen>
<solver_type>BiCGSTAB</solver_type>
<precon_type>ILUT</precon_type>

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

tolerances and number of steps are missing.

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

tolerances and number of steps have been added.

<linear_solver>linear_solver_d</linear_solver>
</nonlinear_solver>
<nonlinear_solver>
<name>basic_newton_u</name>

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

Not sure here. Both nonlinear solvers are configured the same way. Is it possible to have only one?

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

It is certainly possible. They are deliberately defined separately so that one can easily adjust each using this template.

<repeat>10</repeat>
<each_steps>1</each_steps>
</pair>
<!--<pair><repeat>50</repeat><each_steps>100</each_steps></pair>-->

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

drop unused comment

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

dropped.

<convergence_criterion>
<type>DeltaX</type>
<norm_type>NORM2</norm_type>
<reltol>1.e-4</reltol>

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

Would a tighter tolerance be possible? If this is a simple test I'd think something around 1e-14 is easily achievable. The tighter the tolerances, here and in the ctest, the higher chance there is to spot a change in the process' behaviour.

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

The tolerances have been made tighter (1e-14).

/// process has process_id == 0 in the staggered scheme.
bool hasMechanicalProcess(int const process_id) const
{
return _use_monolithic_scheme || process_id == 0;

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

Move implementation into the implementation file.

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

It is the consistent structure with the HydroMechanics process.

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 24, 2018

Author Contributor

This function, hasMechanicalProcess has been removed.

/// Check whether the process represented by \c process_id is/has
/// mechanical process. In the present implementation, the mechanical
/// process has process_id == 0 in the staggered scheme.
bool hasMechanicalProcess(int const process_id) const

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

From the comment, I'd say this function answers the question isMechanicalProcess, not has.
Where is this function used?

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 24, 2018

Author Contributor

It was used by postNonlienar but this function has been removed as above. Thus, hasMechanicalProcess has been removed as well.

@@ -46,7 +46,7 @@ struct IntegrationPointData final

typename BMatricesType::KelvinVectorType sigma_tensile, sigma_compressive,
sigma_real_prev, sigma_real;
double strain_energy_tensile;
double strain_energy_tensile = 0.0;

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

Because none of the members is initialized here, I'd recommend to move the initializations together.

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

It has been moved to the initialisations.

Kud.noalias() += Kud_ip_contribution;

if (history_variable_prev < strain_energy_tensile) {
history_variable = strain_energy_tensile;

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

Formatting is wrong.

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

clang-format has been re-run.


namespace ProcessLib
{
/// TODO docu

This comment has been minimized.

@endJunction

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

Removed the comment.

auto const& nodes = _mesh.getNodes();
for (auto const* n : nodes)
{
std::size_t id = n->getID();

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

node_id would be a better name

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

changed to node_id

const auto g_idx =
_dof_table.getGlobalIndex(l, _variable_id, _component_id);

if (x[id] <= irrevDamage)

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

mostly taste, but full words are easier to understand. What is easier on the brain:
irreversibleDamage or
irvDmg ?

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

changed to irreversibleDamage

}

// update new values and corresponding indices.
void SolutionDependentDirichletBoundaryCondition::preTimestep(

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

This is very specific for the phase field process. It would be better either to generalize this boundary condition, such that the process dependent code resides in the process implementation, or at least rename it something like PhaseFieldIrreversibleDamageOracleBoundaryCondition.

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

the class name has been changed to PhaseFieldIrreversibleDamageOracleBoundaryCondition

}

virtual ~BoundaryCondition() = default;
//! Applies natural BCs (i.e. non-Dirichlet BCs) to the stiffness matrix

This comment has been minimized.

@endJunction

endJunction Jan 22, 2018

Member

Please revert the formatting to the original. There is .clang-format file in the source code tree which must be used with clang-format. Otherwise some other default formatting is used. There must be an option in qtcreator to set the formatting style to "file".

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Jan 23, 2018

Author Contributor

clang-format has been re-run.

This comment has been minimized.

@wenqing

wenqing Feb 9, 2018

Member

It looks that the clang-format in this file has been changed.

@endJunction

This comment has been minimized.

Copy link
Member

endJunction commented Jan 22, 2018

Only superficial comments, didn't look into the code, yet.

@endJunction endJunction referenced this pull request Jan 23, 2018

Closed

Staggered pf #2051

@endJunction

This comment has been minimized.

Copy link
Member

endJunction commented Feb 8, 2018

@KeitaYoshioka Please have a look on the rebased branch in my repo. There are fixes for formatting and the fixed size matrices included.
Still open is the docu, and there is a TODO commit---the reason why the ctests do not pass.

I was not sure if you want to include the hydro_cracking option here; it's in the rebased branch.

@endJunction endJunction added the WIP 🏗 label Feb 8, 2018

}

virtual ~BoundaryCondition() = default;
//! Applies natural BCs (i.e. non-Dirichlet BCs) to the stiffness matrix

This comment has been minimized.

@wenqing

wenqing Feb 9, 2018

Member

It looks that the clang-format in this file has been changed.

void PhaseFieldIrreversibleDamageOracleBoundaryCondition::preTimestep(
const double /*t*/, const GlobalVector& x)
{
double irreversibleDamage = 0.05;

This comment has been minimized.

@wenqing

wenqing Feb 9, 2018

Member

A comment is needed for the number of 0.05.

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Feb 14, 2018

Author Contributor

comments were added.

DBUG(
"Constructing PhaseFieldIrreversibleDamageOracleBoundaryCondition from "
"config.");
//! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}

This comment has been minimized.

@wenqing

wenqing Feb 9, 2018

Member

Parameter doc files are not generated.

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Feb 14, 2018

Author Contributor

done.

else
{
// For the staggered scheme
if (_coupled_solutions->process_id == 0)

This comment has been minimized.

@wenqing

wenqing Feb 9, 2018

Member

Change the process orders as that in assembleWithJacobianConcreteProcess, thought the Picard scheme is not used in the current implementation.

This comment has been minimized.

@KeitaYoshioka

KeitaYoshioka Feb 14, 2018

Author Contributor

The process order is now consistent.

@KeitaYoshioka KeitaYoshioka force-pushed the KeitaYoshioka:StaggeredPF branch from 38dab64 to ef1cd2f Feb 14, 2018

@@ -38,7 +38,7 @@ set(OGS_CPU_ARCHITECTURE "native" CACHE STRING "Processor architecture, \
option(OGS_ENABLE_AVX2 "Enable the use of AVX2 instructions" OFF)
option(OGS_BUILD_TESTS "Should the test executables be built?" ON)
option(OGS_USE_PCH "Should pre-compiled headers be used?" ON)
option(OGS_USE_CONAN "Should Conan package manager be used?" OFF)
option(OGS_USE_CONAN "Should Conan package manager be used?" ON)

This comment has been minimized.

@endJunction

endJunction Feb 15, 2018

Member

I'd think you didn't want commit this change.

@wenqing
Copy link
Member

wenqing left a comment

This PR could be separated into 3 PRs about

  1. staggeredPF, the original one, can be merged with the recent improvement.
  2. PETSc SNES.
  3. staggeredHydroPF
@@ -336,6 +339,16 @@ createNonlinearSolver(GlobalLinearSolver& linear_solver,
std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
tag);
}
OGS_FATAL("Unsupported nonlinear solver type");
#ifdef USE_PETSC

This comment has been minimized.

@wenqing

wenqing Feb 15, 2018

Member

Introducing SNES solver is nice. Could you please separate the SNES stuff from this PR? For example, a PR for SNES only with a unit test, which can be merged before the original staggeredPF.

: _linear_solver(linear_solver)
{
SNESCreate(PETSC_COMM_WORLD, &_snes_solver);
// SNESSetType(_snes_solver, "vi");

This comment has been minimized.

@wenqing

wenqing Feb 15, 2018

Member

Is this comment needed?


// TODO Maybe need to overwrite the ogs x with petsc_x.
/*
DBUG("BEFORE ASSEMBLY")

This comment has been minimized.

@wenqing

wenqing Feb 15, 2018

Member

These lines can be removed.

context->system->applyKnownSolutionsNewton(*context->J, *context->r,
*context->x);
/*
DBUG("AFTER ASSEMBLY");

This comment has been minimized.

@wenqing

wenqing Feb 15, 2018

Member

The same.

VecCopy(context->r->getRawVector(), petsc_r);

/*
DBUG("The petsc-x vector.")

This comment has been minimized.

@wenqing

wenqing Feb 15, 2018

Member

The same.

DBUG("PETScNonlinearSolver: set constraints");
// Constraints
Vec xl, xu;
VecDuplicate(x.getRawVector(), &xl);

This comment has been minimized.

@wenqing

wenqing Feb 15, 2018

Member

Frequently allocating/releasing memory for xl, xu, J, r, x may slow down the computation. It could be better to set these variables as the members of this class, and allocate memory for them only once.

_equation_system->assemble(x);
}

bool solve(GlobalVector& x,

This comment has been minimized.

@wenqing

wenqing Feb 15, 2018

Member

It is better to move the definition of this function to cpp file.

MeshLib::Mesh const& mesh, int const variable_id,
int const component_id)
:

This comment has been minimized.

@wenqing

wenqing Feb 15, 2018

Member

too much space here. Need run clang-format.

@@ -0,0 +1 @@
---

This comment has been minimized.

@wenqing

wenqing Feb 15, 2018

Member

This file has to be removed.

@@ -0,0 +1,3 @@
Start testing: Jan 22 16:49 CET
----------------------------------------------------------
End testing: Jan 22 16:49 CET

This comment has been minimized.

@wenqing

wenqing Feb 15, 2018

Member

This file has to be removed.

@endJunction endJunction force-pushed the KeitaYoshioka:StaggeredPF branch from ef1cd2f to dc0abf9 Mar 1, 2018

@endJunction

This comment has been minimized.

Copy link
Member

endJunction commented Mar 1, 2018

Staggerd PhaseField only. The first four commits are independent PR #2091.
Doxygen/web docu and ctests might need some attention, of which I'll take care today, and when ready set the label to "please review".

@endJunction endJunction force-pushed the KeitaYoshioka:StaggeredPF branch from dc0abf9 to 439f0fc Mar 5, 2018

@endJunction endJunction added please review and removed WIP 🏗 labels Mar 5, 2018

@endJunction

This comment has been minimized.

Copy link
Member

endJunction commented Mar 5, 2018

The test documentation will be provided when the corresponding paper is accepted.

@endJunction endJunction force-pushed the KeitaYoshioka:StaggeredPF branch from 439f0fc to f3b810f Mar 6, 2018

@endJunction endJunction force-pushed the KeitaYoshioka:StaggeredPF branch from f3b810f to ced5d85 Mar 6, 2018

@wenqing

wenqing approved these changes Mar 7, 2018

Copy link
Member

wenqing left a comment

Only minor changes needed.

void preTimestep(const double t, GlobalVector const& x)
{
auto const n_bcs = _boundary_conditions.size();
for (std::size_t i = 0; i < n_bcs; ++i)

This comment has been minimized.

@wenqing

wenqing Mar 7, 2018

Member

Why not use for range loop? ✔️

@@ -56,24 +55,24 @@ std::unique_ptr<Process> createPhaseFieldProcess(
auto per_process_variables = findProcessVariables(
variables, pv_config,
{//! \ogs_file_param_special{prj__processes__process__PHASE_FIELD__process_variables__phasefield}

This comment has been minimized.

@wenqing

wenqing Mar 7, 2018

Member

The parameter documentation file is missing.

This comment has been minimized.

@endJunction

endJunction Mar 7, 2018

Member

Maybe I'm missing something, but Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/process_variables/t_phasefield.md is there....

This comment has been minimized.

@wenqing

wenqing Mar 7, 2018

Member

Sorry. It already exists.


for (auto const& ip_data : _ip_data)
for (unsigned ip = 0; ip < num_intpts; ++ip)

This comment has been minimized.

@wenqing

wenqing Mar 7, 2018

Member

unsigned --> int

@@ -191,10 +192,29 @@ std::unique_ptr<Process> createPhaseFieldProcess(
std::copy_n(b.data(), b.size(), specific_body_force.data());
}

auto const crack_scheme =
//! \ogs_file_param{prj__processes__process__PHASE_FIELD__hydro_crack_scheme}

This comment has been minimized.

@wenqing

wenqing Mar 7, 2018

Member

This keyword seems not be for this PR, but It is OK to keep it.

This comment has been minimized.

@endJunction

endJunction Mar 7, 2018

Member

Yes, you are right. This is already for crack propagation but it is used in Tests/Data/PhaseField/2D_static.prj.

@endJunction endJunction force-pushed the KeitaYoshioka:StaggeredPF branch from ced5d85 to bb9b871 Mar 7, 2018

@endJunction endJunction merged commit 3351c3e into ufz:master Mar 7, 2018

2 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/jenkins/pr-merge This commit looks good
Details

@KeitaYoshioka KeitaYoshioka deleted the KeitaYoshioka:StaggeredPF branch Mar 7, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.