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

Jump to bottom

Fixup newton solver and elastic rheology #5580

Merged
merged 21 commits into from
Apr 11, 2024
Merged

Fixup newton solver and elastic rheology #5580

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
8fa9034
modify the SPD factor
Feb 17, 2024
e767c30
change elastic_force from Tensor to SymmetricTensor and add viscoelas…
Feb 17, 2024
638630b
modify functions Elasticity::fill_elastic_force_outputs() and Elastic…
Feb 17, 2024
852b9a0
apply deviator() to the strain rate in system jacobian and include th…
Feb 17, 2024
24fe57b
restore the original form of Newton residual
Feb 27, 2024
ff6aa4f
Update include/aspect/material_model/interface.h
YiminJin Feb 29, 2024
4f30907
Update include/aspect/material_model/rheology/elasticity.h
YiminJin Feb 29, 2024
f49375b
make the Newton method compatible with material averaging
Mar 4, 2024
4c8a936
make the viscosity derivatives compatible with material averaging
Mar 6, 2024
19d0e36
make the elastic force compatible with material averaging
Mar 6, 2024
e265b12
update GMG solver with newton method
Mar 6, 2024
ff98dd6
fix bugs
Mar 7, 2024
cac47c4
modify the comments
Mar 7, 2024
0d5edfe
indent
Mar 7, 2024
bd3dc15
make function local_apply() more readable
Mar 8, 2024
8aaa42f
apply the changes of test results
Mar 8, 2024
70a53d6
indent
Mar 8, 2024
2d797c5
discard 'project to Q1' when filling elastic outputs and reaction ter…
Mar 9, 2024
b78e839
apply changes in test results
Mar 9, 2024
21e5eb4
Update source/simulator/assemblers/newton_stokes.cc
YiminJin Mar 20, 2024
fa2be51
Update source/simulator/assemblers/newton_stokes.cc
YiminJin Mar 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
11 changes: 9 additions & 2 deletions include/aspect/material_model/interface.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1151,7 +1151,8 @@ namespace aspect
{
public:
ElasticOutputs(const unsigned int n_points)
: elastic_force(n_points, numbers::signaling_nan<Tensor<2,dim>>() )
: elastic_force(n_points, numbers::signaling_nan<SymmetricTensor<2,dim>>())
, viscoelastic_strain_rate(n_points, numbers::signaling_nan<SymmetricTensor<2,dim>>())
{}

~ElasticOutputs() override
Expand All @@ -1170,7 +1171,13 @@ namespace aspect
* momentum equation (first part of the Stokes equation) in each
* quadrature point.
*/
std::vector<Tensor<2,dim>> elastic_force;
std::vector<SymmetricTensor<2,dim>> elastic_force;

/**
* Strain tensor that includes the contribution of elastic shear, which is
YiminJin marked this conversation as resolved.
Show resolved Hide resolved
* required by the Newton solver.
*/
std::vector<SymmetricTensor<2,dim>> viscoelastic_strain_rate;
};


Expand Down
10 changes: 5 additions & 5 deletions include/aspect/material_model/rheology/elasticity.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,9 +89,9 @@ namespace aspect
* fill an additional material model outputs objects with the elastic force terms.
*/
void
fill_elastic_force_outputs (const MaterialModel::MaterialModelInputs<dim> &in,
const std::vector<double> &average_elastic_shear_moduli,
MaterialModel::MaterialModelOutputs<dim> &out) const;
fill_elastic_outputs (const MaterialModel::MaterialModelInputs<dim> &in,
anne-glerum marked this conversation as resolved.
Show resolved Hide resolved
const std::vector<double> &average_elastic_shear_moduli,
MaterialModel::MaterialModelOutputs<dim> &out) const;

/**
* Given the stress of the previous time step in the material model inputs @p in,
Expand Down Expand Up @@ -128,7 +128,7 @@ namespace aspect
const double shear_modulus) const;

/**
* Calculate the norm of the effective deviatoric strain rate tensor,
* Calculate the effective deviatoric strain rate tensor,
* which equals the norm of the true deviatoric strain rate plus
YiminJin marked this conversation as resolved.
Show resolved Hide resolved
* a fictional strain rate which would arise from stored elastic stresses.
* In ASPECT, this additional strain rate is
Expand All @@ -137,7 +137,7 @@ namespace aspect
* by ensuring that the resulting strain rate tensor is equal to the
* total current stress tensor multiplied by a scalar.
*/
double
SymmetricTensor<2,dim>
calculate_viscoelastic_strain_rate (const SymmetricTensor<2,dim> &strain_rate,
const SymmetricTensor<2,dim> &stored_stress,
const double shear_modulus) const;
Expand Down
6 changes: 6 additions & 0 deletions include/aspect/newton.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,6 +234,12 @@ namespace aspect
void
execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base,
internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;

/**
* Create additional material models outputs for computing viscoelastic strain rate when
* elasticity is enabled.
*/
void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const override;
};

/**
Expand Down
31 changes: 13 additions & 18 deletions source/material_model/rheology/elasticity.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -242,16 +242,15 @@ namespace aspect

template <int dim>
void
Elasticity<dim>::fill_elastic_force_outputs (const MaterialModel::MaterialModelInputs<dim> &in,
const std::vector<double> &average_elastic_shear_moduli,
MaterialModel::MaterialModelOutputs<dim> &out) const
Elasticity<dim>::fill_elastic_outputs (const MaterialModel::MaterialModelInputs<dim> &in,
anne-glerum marked this conversation as resolved.
Show resolved Hide resolved
const std::vector<double> &average_elastic_shear_moduli,
MaterialModel::MaterialModelOutputs<dim> &out) const
{
// Create a reference to the structure for the elastic force terms that are needed to compute the
// right-hand side of the Stokes system
// Create a reference to the structure for the elastic outputs
MaterialModel::ElasticOutputs<dim>
*force_out = out.template get_additional_output<MaterialModel::ElasticOutputs<dim>>();
*elastic_out = out.template get_additional_output<MaterialModel::ElasticOutputs<dim>>();

if (force_out == nullptr)
if (elastic_out == nullptr)
return;

if (in.requests_property(MaterialProperties::additional_outputs))
Expand All @@ -262,12 +261,10 @@ namespace aspect
for (unsigned int j=0; j < SymmetricTensor<2,dim>::n_independent_components; ++j)
stress_old[SymmetricTensor<2,dim>::unrolled_to_component_indices(j)] = in.composition[i][j];

// Average viscoelastic viscosity
const double average_viscoelastic_viscosity = out.viscosities[i];

// Fill elastic force outputs (See equation 30 in Moresi et al., 2003, J. Comp. Phys.)
force_out->elastic_force[i] = -1. * ( average_viscoelastic_viscosity / calculate_elastic_viscosity(average_elastic_shear_moduli[i]) * stress_old );

// The viscoelastic strain rate is needed only when the Newton method or the DCP method is used for solving.
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
// The viscoelastic strain rate is needed only when the Newton method or the DCP method is used for solving.
// The viscoelastic strain rate is needed only when the Newton or defect correction picard nonlinear solvers are selected.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@naliboff Thanks very much for reviewing my codes!

I have tried the benchmark and the modified Newton solver performs poorly :-(. I think the problem is caused by the following two facts:

  1. In the first 4 time steps the solver should converge within 2 iterations since plastic yielding has not occurred and the system is linear. However, the modified Newton solver needs a lot more iterations to converge. This is because on the RHS of the Newton system I have changed
    $$2\eta\dot{\boldsymbol{\varepsilon}} - \frac{\eta}{G\Delta t}\boldsymbol{\tau}^{\text{old}}$$
    to
    $$2\eta\dot{\boldsymbol{\varepsilon}}_{ve},$$

which looks identical. However, the viscoelastic strain rate is calculated by
$$\dot{\boldsymbol{\varepsilon}}_{ve} = \text{dev}(\dot{\boldsymbol{\varepsilon}}) + \frac{\boldsymbol{\tau}^{\text{old}}}{2G\Delta t},$$
where deviator() is applied to the current strain rate (this is necessary, since the viscosity is calculated with the deviatoric strain rate). Consequently, the RHS is not precisely equal to the system residual. This problem has been fixed by restoring the original RHS, and the covergence rate is much better (see the log files below), but still not as fast as expected.

  1. The GMG solver has been selected in this benchmark, which requires the viscosity to be averaged in each cell. However, function MaterialModel::MaterialAveraging::average() is called AFTER MaterialModel::evaluate(), which means that $\eta$ is changed after $\partial\eta/\partial\dot{\boldsymbol{\varepsilon}}$ is calculated, hence the LHS can no longer represent the system tangent. Therefore, the modified Newton solver cannot be used simultaneously with GMG solver for now. To fix this problem we may have to do a lot more work (to calculate $\partial\eta/\partial\dot{\boldsymbol{\varepsilon}}$ after cell averaging).

The following log files are the results of the main branch, this PR (after fixing problem 1), and without cell averaging. The convergence rate without cell averaging is much faster than the previous two.
log_main.txt
log_pr5580.txt
log_no_avg.txt

elastic_out->elastic_force[i] = -out.viscosities[i] / calculate_elastic_viscosity(average_elastic_shear_moduli[i]) * stress_old;
if (Parameters<dim>::is_defect_correction(this->get_parameters().nonlinear_solver))
Copy link
Contributor

Choose a reason for hiding this comment

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

Will this condition also capture the case for the full Newton solver?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes. In fact, after restoring the RHS only the full Newton solver needs the viscoelastic strain rate now.

Copy link
Contributor

Choose a reason for hiding this comment

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

So does this condition need to be updated to only reflect the full Newton solver?

elastic_out->viscoelastic_strain_rate[i] = calculate_viscoelastic_strain_rate(in.strain_rate[i], stress_old, average_elastic_shear_moduli[i]);
}
}

Expand Down Expand Up @@ -412,7 +409,7 @@ namespace aspect


template <int dim>
double
SymmetricTensor<2,dim>
Elasticity<dim>::
calculate_viscoelastic_strain_rate(const SymmetricTensor<2,dim> &strain_rate,
const SymmetricTensor<2,dim> &stored_stress,
Expand All @@ -425,10 +422,8 @@ namespace aspect
// elastic stresses stored from the last time step.
// Note the parallels with the viscous part of the strain rate deviator,
// which is equal to 0.5 * stress / viscosity.
const SymmetricTensor<2,dim> edot_deviator = deviator(strain_rate) + 0.5*stored_stress /
calculate_elastic_viscosity(shear_modulus);
// Return the norm of the strain rate, or 0, whichever is larger.
return std::sqrt(std::max(-second_invariant(edot_deviator), 0.));
return deviator(strain_rate) + 0.5 * deviator(stored_stress) /
calculate_elastic_viscosity(shear_modulus);
}
}
}
Expand Down
19 changes: 9 additions & 10 deletions source/material_model/rheology/visco_plastic.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -261,11 +261,12 @@ namespace aspect
Assert(std::isfinite(in.strain_rate[i].norm()),
ExcMessage("Invalid strain_rate in the MaterialModelInputs. This is likely because it was "
"not filled by the caller."));
const double effective_strain_rate_invariant = elastic_rheology.calculate_viscoelastic_strain_rate(in.strain_rate[i],
stress_old,
elastic_shear_moduli[j]);
const SymmetricTensor<2,dim> effective_strain_rate =
elastic_rheology.calculate_viscoelastic_strain_rate(in.strain_rate[i],
stress_old,
elastic_shear_moduli[j]);

effective_edot_ii = std::max(effective_strain_rate_invariant,
effective_edot_ii = std::max(std::sqrt(std::max(-second_invariant(effective_strain_rate), 0.)),
min_strain_rate);
}

Expand Down Expand Up @@ -403,17 +404,15 @@ namespace aspect
ExcMessage("Invalid strain_rate in the MaterialModelInputs. This is likely because it was "
"not filled by the caller."));

const SymmetricTensor<2,dim> deviatoric_strain_rate = deviator(in.strain_rate[i]);

// For each independent component, compute the derivative.
for (unsigned int component = 0; component < SymmetricTensor<2,dim>::n_independent_components; ++component)
{
const TableIndices<2> strain_rate_indices = SymmetricTensor<2,dim>::unrolled_to_component_indices (component);

// components that are not on the diagonal are multiplied by 0.5, because the symmetric tensor
// is modified by 0.5 in both symmetric directions (xy/yx) simultaneously and we compute the combined
// derivative
const SymmetricTensor<2,dim> strain_rate_difference = in.strain_rate[i]
+ std::max(std::fabs(in.strain_rate[i][strain_rate_indices]), min_strain_rate)
* (component > dim-1 ? 0.5 : 1 )
anne-glerum marked this conversation as resolved.
Show resolved Hide resolved
const SymmetricTensor<2,dim> strain_rate_difference = deviatoric_strain_rate
+ std::max(std::fabs(deviatoric_strain_rate[strain_rate_indices]), min_strain_rate)
* finite_difference_accuracy
* Utilities::nth_basis_for_symmetric_tensors<dim>(component);

Expand Down
2 changes: 1 addition & 1 deletion source/material_model/visco_plastic.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -306,7 +306,7 @@ namespace aspect

if (this->get_parameters().enable_elasticity)
{
rheology->elastic_rheology.fill_elastic_force_outputs(in, average_elastic_shear_moduli, out);
rheology->elastic_rheology.fill_elastic_outputs(in, average_elastic_shear_moduli, out);
rheology->elastic_rheology.fill_reaction_outputs(in, average_elastic_shear_moduli, out);
}
}
Expand Down
2 changes: 1 addition & 1 deletion source/material_model/viscoelastic.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ namespace aspect
}
}

elastic_rheology.fill_elastic_force_outputs(in, average_elastic_shear_moduli, out);
elastic_rheology.fill_elastic_outputs(in, average_elastic_shear_moduli, out);
elastic_rheology.fill_reaction_outputs(in, average_elastic_shear_moduli, out);

}
Expand Down
Loading
Loading