Skip to content

Commit ed987d7

Browse files
MFratersMFraters
committed
Add direction dilation.
1 parent e0cfd33 commit ed987d7

File tree

5 files changed

+80
-20
lines changed

5 files changed

+80
-20
lines changed

‎include/aspect/material_model/interface.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1247,10 +1247,10 @@ namespace aspect
12471247
std::vector<double> dilation_lhs_term;
12481248

12491249
/**
1250-
* A scalar value per evaluation point corresponding to the RHS term
1250+
* A scalar value per direction per evaluation point corresponding to the RHS term
12511251
* due to plastic dilation.
12521252
*/
1253-
std::vector<double> dilation_rhs_term;
1253+
std::vector<std::vector<double>> dilation_rhs_term;
12541254
};
12551255

12561256

‎source/material_model/interface.cc‎

Lines changed: 12 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1005,6 +1005,9 @@ namespace aspect
10051005
std::vector<std::string> names;
10061006
names.emplace_back("dilation_lhs_term");
10071007
names.emplace_back("dilation_rhs_term");
1008+
names.emplace_back("dilation_rhs_term");
1009+
if (dim == 3)
1010+
names.emplace_back("dilation_rhs_term");
10081011
return names;
10091012
}
10101013
}
@@ -1015,22 +1018,28 @@ namespace aspect
10151018
PrescribedPlasticDilation<dim>::PrescribedPlasticDilation (const unsigned int n_points)
10161019
: NamedAdditionalMaterialOutputs<dim>(make_prescribed_dilation_outputs_names()),
10171020
dilation_lhs_term(n_points, numbers::signaling_nan<double>()),
1018-
dilation_rhs_term(n_points, numbers::signaling_nan<double>())
1021+
dilation(dim, std::vector<double>(n_points, numbers::signaling_nan<double>()))
10191022
{}
10201023

10211024

10221025

10231026
template <int dim>
10241027
std::vector<double> PrescribedPlasticDilation<dim>::get_nth_output(const unsigned int idx) const
10251028
{
1026-
AssertIndexRange (idx, 2);
1029+
AssertIndexRange (idx, 4);
10271030
switch (idx)
10281031
{
10291032
case 0:
10301033
return dilation_lhs_term;
10311034

10321035
case 1:
1033-
return dilation_rhs_term;
1036+
return dilation_rhs_term[0];
1037+
1038+
case 2:
1039+
return dilation_rhs_term[1];
1040+
1041+
case 3:
1042+
return dilation_rhs_term[2];
10341043

10351044
default:
10361045
AssertThrow(false, ExcInternalError());

‎source/simulator/assemblers/newton_stokes.cc‎

Lines changed: 32 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -461,13 +461,34 @@ namespace aspect
461461
// dilation term, but also the LHS dilation term should be included in the
462462
// system residual
463463
if (enable_prescribed_dilation)
464-
data.local_rhs(i) += (
465-
- pressure_scaling
466-
* (prescribed_dilation->dilation_rhs_term[q] -
467-
prescribed_dilation->dilation_lhs_term[q] *
468-
scratch.material_model_inputs.pressure[q])
469-
* scratch.phi_p[i]
470-
) * JxW;
464+
{
465+
const unsigned int index_direction=fe.system_to_component_index(i).first;
466+
// Only want the velocity components and not the pressure one (which is the last one), so add 1
467+
if (introspection.is_stokes_component(index_direction+1))
468+
data.local_rhs(i) += (
469+
// RHS of - (div u,q) = - (R,q)
470+
- pressure_scaling
471+
* (prescribed_dilation->dilation_rhs_term[index_direction][q] -
472+
prescribed_dilation->dilation_lhs_term[q] *
473+
scrasch.material_model_inputs.pressure[q])
474+
* scratch.phi_p[i]
475+
) * JxW;
476+
}
477+
478+
// Only assemble this term if we are running incompressible, otherwise this term
479+
// is already included on the LHS of the equation.
480+
if (enable_prescribed_dilation && !material_model_is_compressible)
481+
{
482+
const unsigned int index_direction=fe.system_to_component_index(i).first;
483+
// Only want the velocity components and not the pressure one (which is the last one), so add 1
484+
if (introspection.is_stokes_component(index_direction+1))
485+
data.local_rhs(i) += (
486+
// RHS of momentum eqn: - \int 2/3 eta R, div v
487+
- 2.0 / 3.0 * eta
488+
* prescribed_dilation->dilation_rhs_term[index_direction][q]
489+
* scratch.div_phi_u[i]
490+
) * JxW;
491+
}
471492
}
472493

473494
// and then the matrix, if necessary
@@ -665,7 +686,10 @@ namespace aspect
665686
Assert(!this->get_parameters().enable_prescribed_dilation
666687
||
667688
(outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation_lhs_term.size() == n_points &&
668-
outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation_rhs_term.size() == n_points),
689+
outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation_rhs_term.size() == dim &&
690+
outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation_rhs_term[0].size() == n_points &&
691+
outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation_rhs_term[1].size() == n_points &&
692+
(dim == 2 || (dim == 3 && outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation_rhs_term[2].size() == n_points))),
669693
ExcInternalError());
670694

671695
if (this->get_newton_handler().parameters.newton_derivative_scaling_factor != 0)

‎source/simulator/assemblers/stokes.cc‎

Lines changed: 30 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -451,11 +451,32 @@ namespace aspect
451451
* JxW;
452452

453453
if (prescribed_dilation != nullptr)
454-
data.local_rhs(i) += (
455-
- pressure_scaling
456-
* prescribed_dilation->dilation_rhs_term[q]
457-
* scratch.phi_p[i]
458-
) * JxW;
454+
{
455+
const unsigned int index_direction=fe.system_to_component_index(i).first;
456+
// Only want the velocity components and not the pressure one (which is the last one), so add 1
457+
if (introspection.is_stokes_component(index_direction+1))
458+
data.local_rhs(i) += (
459+
// RHS of - (div u,q) = - (R,q)
460+
- pressure_scaling
461+
* prescribed_dilation->dilation_rhs_term[index_direction][q]
462+
* scratch.phi_p[i]
463+
) * JxW;
464+
}
465+
466+
// Only assemble this term if we are running incompressible, otherwise this term
467+
// is already included on the LHS of the equation.
468+
if (prescribed_dilation != nullptr && !material_model_is_compressible)
469+
{
470+
const unsigned int index_direction=fe.system_to_component_index(i).first;
471+
// Only want the velocity components and not the pressure one (which is the last one), so add 1
472+
if (introspection.is_stokes_component(index_direction+1))
473+
data.local_rhs(i) += (
474+
// RHS of momentum eqn: - \int 2/3 eta R, div v
475+
- 2.0 / 3.0 * eta
476+
* prescribed_dilation->dilation_rhs_term[index_direction][q]
477+
* scratch.div_phi_u[i]
478+
) * JxW;
479+
}
459480

460481
if (scratch.rebuild_stokes_matrix)
461482
for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
@@ -531,7 +552,10 @@ namespace aspect
531552
Assert(!this->get_parameters().enable_prescribed_dilation
532553
||
533554
(outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation_lhs_term.size() == n_points &&
534-
outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation_rhs_term.size() == n_points),
555+
outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation_rhs_term.size() == dim &&
556+
outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation_rhs_term[0].size() == n_points &&
557+
outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation_rhs_term[1].size() == n_points &&
558+
(dim == 2 || (dim == 3 && outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation_rhs_term[2].size() == n_points))),
535559
ExcInternalError());
536560

537561
// Elasticity:

‎tests/prescribed_dilation.cc‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -208,7 +208,10 @@ namespace aspect
208208
if (prescribed_dilation)
209209
{
210210
prescribed_dilation->dilation_rhs_term[i] = x;
211-
prescribed_dilation->dilation_lhs_term[i] = 0.;
211+
prescribed_dilation->dilation_lhs_term[0][i] = 0.;
212+
prescribed_dilation->dilation_lhs_term[1][i] = 0.;
213+
if (dim == 3)
214+
prescribed_dilation->dilation_lhs_term[2][i] = 0.;
212215
}
213216

214217
}

0 commit comments

Comments
 (0)