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 thermo-mechanical phase-field #2102

Merged
merged 59 commits into from May 22, 2018

Conversation

Projects
None yet
3 participants
@xingyuanmiao
Copy link
Contributor

xingyuanmiao commented Mar 15, 2018

The thermo-mechanical phase-field process is implemented. The implementation includes the following tests: a 3D beam undergoing thermal shrinkage, a 3D cube for homogeneous damage; two large tests, one for iglu project, one for thermal shock testing.

@xingyuanmiao xingyuanmiao force-pushed the xingyuanmiao:staggeredThmermoMechanicalPhaseField branch 3 times, most recently from d48f361 to 1c44568 Mar 15, 2018

@wenqing

This comment has been minimized.

Copy link
Member

wenqing commented Mar 16, 2018

@xingyuanmiao You can merge my PR to your branch for this PR.

int mechanics_related_process_id = 0;
int phase_field_process_id = 0;
int heat_conduction_process_id = 0;
if (use_monolithic_scheme) // monolithic scheme.

This comment has been minimized.

@endJunction

endJunction Mar 17, 2018

Member

Since you did not implement the monolithic scheme, I'd replace these code parts with OGS_FATAL calls with explanatory messages.
Update: Since you actually have implemented the monolithic scheme, it might be worth keep it, but you have to add ctests for that case.

This comment has been minimized.

@endJunction

endJunction Apr 10, 2018

Member

This is still not addressed.

This comment has been minimized.

@wenqing

wenqing Apr 11, 2018

Member

An OGS_FATAL is already in line 43 for the monolithic scheme. If the monolithic scheme is not interesting for solving of coupled TMPh problems , all monolithic scheme related stuff in this PR can be removed.

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 12, 2018

Author Contributor

I agree with Wenqing, judging by my experience on testing specific cases for model verification and applications, the monolithic scheme is not an efficient method for solving complex cases.

MathLib::KelvinVector::KelvinVectorDimensions<
DisplacementDim>::value;
ip_data.eps.setZero(kelvin_vector_size);
ip_data.eps_prev.setZero(kelvin_vector_size);

This comment has been minimized.

@endJunction

endJunction Mar 17, 2018

Member

Initialization of the _prev quantities is not necessary. Better to use resize like it is done in PhaseFieldProcess. Initialization of these quantities conseils improper "push_back_state" implementation or other initial conditions implementation. ✔️

typename BMatricesType::KelvinVectorType eps_m, eps_m_prev;

typename BMatricesType::KelvinVectorType sigma_tensile, sigma_compressive,
sigma_real_prev, sigma_real;

This comment has been minimized.

@endJunction

endJunction Mar 17, 2018

Member

I didn't see where the sigma_real_prev is used; could you show me its usage, please?

This comment has been minimized.

@endJunction

endJunction Mar 17, 2018

Member

Another question, just for understanding: What does the _real postfix mean? Probably it is not about complex numbers, but then there might be a better (more descriptive) name.

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Mar 17, 2018

Author Contributor

if the _real need to be changed, then it should also be renamed in the PhaseField process for consistency.

This comment has been minimized.

@endJunction

endJunction Mar 20, 2018

Member

See #2103 ✔️

typename BMatricesType::KelvinVectorType sigma_tensile, sigma_compressive,
sigma_real_prev, sigma_real;
double strain_energy_tensile, elastic_energy;
typename ShapeMatrixType::GlobalDimVectorType heatflux, heatflux_prev;

This comment has been minimized.

@endJunction

endJunction Mar 17, 2018

Member

I also could not find the usage of heatflux_prev. Am I missing something obvious?

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Mar 17, 2018

Author Contributor

check the pushBackState.

This comment has been minimized.

@endJunction

endJunction Mar 20, 2018

Member

I have seen, that there is an assignment to that variable, but where is it read?
Same question goes on for the history_variable_prev.

This comment has been minimized.

@endJunction

endJunction Apr 10, 2018

Member

The history_variable_prev question is not addressed. ✔️

{
if (history_variable_prev < history_variable)
{
history_variable_prev = history_variable;

This comment has been minimized.

@endJunction

endJunction Mar 17, 2018

Member

minor: could also be

history_variable_prev = std::max(history_variable_prev, history_variable);

For the curious minds, this indeed produces no branches and one less assembly instruction, compared to the if-clause implementation. ✔️

using Invariants = MathLib::KelvinVector::Invariants<
MathLib::KelvinVector::KelvinVectorDimensions<
DisplacementDim>::value>;
double const epsm_trace = Invariants::trace(eps_m);

This comment has been minimized.

@endJunction

endJunction Mar 17, 2018

Member

eps_m_trace

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Mar 17, 2018

Author Contributor

the name has been changed.

typename ShapeMatrixType::GlobalDimNodalMatrixType dNdx;

typename BMatricesType::KelvinVectorType eps, eps_prev;
typename BMatricesType::KelvinVectorType eps_m, eps_m_prev;

This comment has been minimized.

@endJunction

endJunction Mar 17, 2018

Member

Also I don't see usage of eps_m_prev aside of copying and initialization.

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Mar 17, 2018

Author Contributor

the unused parameters are removed. BTW: it was used for updating the mechanical strain.


#include "CreateThermoMechanicalPhaseFieldProcess.h"

#include <cassert>

This comment has been minimized.

@endJunction

endJunction Mar 19, 2018

Member

Please check, which includes are needed (in all files), and remove the unnecessary ones.

This comment has been minimized.

@endJunction

endJunction Mar 20, 2018

Member

I see you have found some other includes, but missed this cassert include..... ✔️


auto const& sigma_real = _ip_data[ip].sigma_real;

// auto const [&](member){ return _process_data.member(t,

This comment has been minimized.

@endJunction

endJunction Mar 19, 2018

Member

Please remove all commented code parts or add a debugging description if something is needed for such a case. ✔️

@xingyuanmiao xingyuanmiao force-pushed the xingyuanmiao:staggeredThmermoMechanicalPhaseField branch 2 times, most recently from adb9d1e to b916c7b Mar 21, 2018


if (history_variable_prev < strain_energy_tensile)
{
// INFO("History variable %g:", history_variable);

This comment has been minimized.

@wenqing

wenqing Mar 21, 2018

Member

If these two INFOs are not needed, please remove it.

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 10, 2018

Author Contributor

removed unnecessary INFOs.


//
// temperature equation, temperature part;
// temperature equation, phasefield part;

This comment has been minimized.

@wenqing

wenqing Mar 21, 2018

Member

Are these comments needed?

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Mar 22, 2018

Author Contributor

These comments are used to distinguish the components of matrix.

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 10, 2018

Author Contributor

removed unnecessary comments.

double const T_ip = N.dot(T);
double const delta_T = T_ip - T0;
// calculate real density
double const rho_s = rho_sr * (1 - 3 * alpha * delta_T);

This comment has been minimized.

@wenqing

wenqing Mar 21, 2018

Member

Where does this density formula come from?

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Mar 22, 2018

Author Contributor

rho_s * (V_0 + dV) = rho_sr * V_0; dV = 3 * alpha * dT * V_0; rho_s = rho_sr / (1 + 3 * alpha * dT )

This comment has been minimized.

@wenqing

wenqing Mar 23, 2018

Member

You can either use this corrected formula by adding the above description as comment, or use rho_sr as rho_s directly.


#include "ThermoMechanicalPhaseFieldFEM.h"
#include "ThermoMechanicalPhaseFieldProcessData.h"
#include "ThermoMechanicalPhaseFieldProcess.h"

This comment has been minimized.

@wenqing

wenqing Mar 21, 2018

Member

Move it to the line before #include <cassert>. In other words, it should be the first include.

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 10, 2018

Author Contributor

done.

*/

#include "ThermoMechanicalPhaseFieldProcess-impl.h"
#include "ThermoMechanicalPhaseFieldProcess.h"

This comment has been minimized.

@wenqing

wenqing Mar 21, 2018

Member

Move #include "ThermoMechanicalPhaseFieldProcess.h" to the line before #include "ThermoMechanicalPhaseFieldProcess-impl.h"

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 10, 2018

Author Contributor

done.

std::vector<double>& local_Jac_data)
{
using RhsVector =
typename ShapeMatricesType::template VectorType<phasefield_size +

This comment has been minimized.

@wenqing

wenqing Mar 21, 2018

Member

If this function is for the monolithic scheme of phase-field thermomechanics, the size of RhsVector or JacobianMatrix should be temperature_size + phasefield_size + displacement_size.

This comment has been minimized.

@wenqing

wenqing Mar 21, 2018

Member

If only the staggered scheme is needed, you can remove this function.

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 10, 2018

Author Contributor

this function is removed as well as the monolithic scheme.


double const d_dot_ip = N.dot(d_dot);

local_rhs.template block<phasefield_size, 1>(phasefield_index, 0)

This comment has been minimized.

@endJunction

endJunction Mar 21, 2018

Member

blocks with 1 columns are usually created using segment

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 10, 2018

Author Contributor

fixed.

i, i * displacement_size / DisplacementDim)
.noalias() = N;

using Invariants = MathLib::KelvinVector::Invariants<

This comment has been minimized.

@endJunction

endJunction Mar 21, 2018

Member

Unused typedef.

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 10, 2018

Author Contributor

removed all unused functions and definitions.

int DisplacementDim>
void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
DisplacementDim>::
assembleWithJacobian(double const t, std::vector<double> const& local_x,

This comment has been minimized.

@endJunction

endJunction Mar 21, 2018

Member

This function is not tested with the ctest, but should be.

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 10, 2018

Author Contributor

the monolithic scheme is removed in current process.

@endJunction

This comment has been minimized.

Copy link
Member

endJunction commented Mar 22, 2018

There also are some doxygen warnings. see Parameter QA page https://jenkins.opengeosys.org/job/ufz/job/ogs/job/PR-2102/12/Doxygen/
✔️

@@ -0,0 +1,27 @@
AddTest(

This comment has been minimized.

@endJunction

endJunction Mar 22, 2018

Member

In the PR description you mention four tests, but there are only two ctests: 3D_non-isothermal and 3D_beam_LARGE. Something is missing....

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 10, 2018

Author Contributor

other two tests including the IGLU TES and thermal shock slab are added in ctests.

* the latest published benchmark book
* "Thermo-Hydro-Mechanical-Chemical Processes in Fractured Porous Media:
* Modelling and Benchmarking"
* (Chapter 9.7 -- A Phase-Field Model for Brittle Fracturing of Thermo-Elastic Solids)

This comment has been minimized.

@endJunction

endJunction Mar 22, 2018

Member

Please add a proper bibtex reference into Documentation/bibliography.bib and \cite it here. ✔️

@@ -439,7 +439,7 @@ void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
double const T_ip = N.dot(T);
double const delta_T = T_ip - T0;
// calculate real density
double const rho_s = rho_sr * (1 - 3 * alpha * delta_T);
double const rho_s = rho_sr / (1 + 3 * alpha * delta_T);

This comment has been minimized.

@endJunction

endJunction Mar 22, 2018

Member

Why is this change not reflected in the ctests? There must be a significant difference, isn't it?

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Mar 22, 2018

Author Contributor
  1. body force is not considered in current testing cases.
  2. under the present values of parameters, the density does not have large changes.
@@ -0,0 +1 @@
An optional input to select the coupling scheme. The monolithic scheme is currently not verified for THERMO_MECHANICAL_PHASE_FIELD.

This comment has been minimized.

@endJunction

endJunction Mar 22, 2018

Member

The monolithic scheme must be validated and verified in this PR. Then the second comment can be removed.

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 10, 2018

Author Contributor

removed the monolithic scheme.

This comment has been minimized.

@endJunction

endJunction Apr 10, 2018

Member

If the monolithic scheme is removed, this is then also obsolete...
Update: Maybe not clear; by "this" I mean the tag <coupling_scheme>

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 12, 2018

Author Contributor

the coupling_scheme option is removed.

@@ -0,0 +1 @@
Residual thermal conductivity represents the property of left medium in cracked region.

This comment has been minimized.

@endJunction

endJunction Mar 22, 2018

Member

Probably "left" is not right, but "remaining".

This comment has been minimized.

@endJunction

endJunction Apr 10, 2018

Member

not addressed

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 12, 2018

Author Contributor

it was rephrased.

@wenqing
Copy link
Member

wenqing left a comment

I am fine with this PR as long as my last four comments are addressed.

DIFF_DATA
expected_beam3d_pcs_2_ts_1_t_10.000000.vtu beam3d_pcs_2_ts_1_t_10.000000.vtu displacement displacement 1e-5 0
expected_beam3d_pcs_2_ts_1_t_10.000000.vtu beam3d_pcs_2_ts_1_t_10.000000.vtu phasefield phasefield 1e-6 0
expected_beam3d_pcs_2_ts_1_t_10.000000.vtu beam3d_pcs_2_ts_1_t_10.000000.vtu temperature temperature 1e-6 0

This comment has been minimized.

@wenqing

wenqing Mar 23, 2018

Member

Heat flux is output but its components are all zero.

This comment has been minimized.

@wenqing

wenqing Mar 23, 2018

Member

It could be nice to take one benchmark from PhaseField and test it with homogeneous temperature field by this implementation.

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 10, 2018

Author Contributor

the benchmark for comparison of pure mechanical phasefield is added.

This comment has been minimized.

@endJunction

endJunction Apr 10, 2018

Member

Could you tell the name of the benchmark, please.

This comment has been minimized.

@xingyuanmiao

xingyuanmiao Apr 12, 2018

Author Contributor

ThermoMechanicalPhaseField_3D_non-isothermal

This comment has been minimized.

@wenqing

wenqing May 18, 2018

Member

ThermoMechanicalPhaseField_3D_non-isothermal --> ThermoMechanicalPhaseField_3D_isothermal

This comment has been minimized.

@xingyuanmiao

xingyuanmiao May 19, 2018

Author Contributor

changed.

double const T_ip = N.dot(T);
double const delta_T = T_ip - T0;
// calculate real density
double const rho_s = rho_sr / (1 + 3 * alpha * delta_T);

This comment has been minimized.

@wenqing

wenqing Mar 23, 2018

Member

See my comment for another rho_s

double const T_dot_ip = N.dot(T_dot);
double const delta_T = T_ip - T0;
// calculate real density
double const rho_s = rho_sr * (1 - 3 * alpha * delta_T);

This comment has been minimized.

@wenqing

wenqing Mar 23, 2018

Member

Correct it or use rho_sr as rho_s directly. For solid of geomaterials or metal, the factor 1 / (1 +3 * alpha * dT) to their density are tiny.

@xingyuanmiao xingyuanmiao force-pushed the xingyuanmiao:staggeredThmermoMechanicalPhaseField branch 2 times, most recently from f71feac to c332b0a Mar 24, 2018

@endJunction endJunction force-pushed the xingyuanmiao:staggeredThmermoMechanicalPhaseField branch from fbdfc94 to e56c47e May 21, 2018

@endJunction

This comment has been minimized.

Copy link
Member

endJunction commented May 22, 2018

Technically the PR is OK. All code parts are exercised by the tests and the complex tests look reasonable.
There is also some documentation, although the notation is not consistent there.
The basic tests need careful checking.

@endJunction endJunction merged commit 63c2ef2 into ufz:master May 22, 2018

1 of 2 checks passed

continuous-integration/jenkins/pr-merge This commit is being built
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
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.