From 78fe922ed8bbde387697d4f0acf0415e2d67cfc7 Mon Sep 17 00:00:00 2001 From: Menno Fraters Date: Sat, 14 Jun 2025 00:15:16 +0200 Subject: [PATCH 1/8] Add direction dilation. --- include/aspect/material_model/interface.h | 4 +- source/material_model/interface.cc | 15 ++++++-- source/simulator/assemblers/newton_stokes.cc | 40 ++++++++++++++++---- source/simulator/assemblers/stokes.cc | 36 +++++++++++++++--- tests/prescribed_dilation.cc | 5 ++- 5 files changed, 80 insertions(+), 20 deletions(-) diff --git a/include/aspect/material_model/interface.h b/include/aspect/material_model/interface.h index 9e9974585a3..6a86a999c9b 100644 --- a/include/aspect/material_model/interface.h +++ b/include/aspect/material_model/interface.h @@ -1247,10 +1247,10 @@ namespace aspect std::vector dilation_lhs_term; /** - * A scalar value per evaluation point corresponding to the RHS term + * A scalar value per direction per evaluation point corresponding to the RHS term * due to plastic dilation. */ - std::vector dilation_rhs_term; + std::vector> dilation_rhs_term; }; diff --git a/source/material_model/interface.cc b/source/material_model/interface.cc index 7c7cd1b4455..7971a76717a 100644 --- a/source/material_model/interface.cc +++ b/source/material_model/interface.cc @@ -1004,6 +1004,9 @@ namespace aspect std::vector names; names.emplace_back("dilation_lhs_term"); names.emplace_back("dilation_rhs_term"); + names.emplace_back("dilation_rhs_term"); + if (dim == 3) + names.emplace_back("dilation_rhs_term"); return names; } } @@ -1014,7 +1017,7 @@ namespace aspect PrescribedPlasticDilation::PrescribedPlasticDilation (const unsigned int n_points) : NamedAdditionalMaterialOutputs(make_prescribed_dilation_outputs_names()), dilation_lhs_term(n_points, numbers::signaling_nan()), - dilation_rhs_term(n_points, numbers::signaling_nan()) + dilation(dim, std::vector(n_points, numbers::signaling_nan())) {} @@ -1022,14 +1025,20 @@ namespace aspect template std::vector PrescribedPlasticDilation::get_nth_output(const unsigned int idx) const { - AssertIndexRange (idx, 2); + AssertIndexRange (idx, 4); switch (idx) { case 0: return dilation_lhs_term; case 1: - return dilation_rhs_term; + return dilation_rhs_term[0]; + + case 2: + return dilation_rhs_term[1]; + + case 3: + return dilation_rhs_term[2]; default: AssertThrow(false, ExcInternalError()); diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc index 516116489fd..0e2870575f7 100644 --- a/source/simulator/assemblers/newton_stokes.cc +++ b/source/simulator/assemblers/newton_stokes.cc @@ -461,13 +461,34 @@ namespace aspect // dilation term, but also the LHS dilation term should be included in the // system residual if (enable_prescribed_dilation) - data.local_rhs(i) += ( - - pressure_scaling - * (prescribed_dilation->dilation_rhs_term[q] - - prescribed_dilation->dilation_lhs_term[q] * - scratch.material_model_inputs.pressure[q]) - * scratch.phi_p[i] - ) * JxW; + { + const unsigned int index_direction=fe.system_to_component_index(i).first; + // Only want the velocity components and not the pressure one (which is the last one), so add 1 + if (introspection.is_stokes_component(index_direction+1)) + data.local_rhs(i) += ( + // RHS of - (div u,q) = - (R,q) + - pressure_scaling + * (prescribed_dilation->dilation_rhs_term[index_direction][q] - + prescribed_dilation->dilation_lhs_term[q] * + scrasch.material_model_inputs.pressure[q]) + * scratch.phi_p[i] + ) * JxW; + } + + // Only assemble this term if we are running incompressible, otherwise this term + // is already included on the LHS of the equation. + if (enable_prescribed_dilation && !material_model_is_compressible) + { + const unsigned int index_direction=fe.system_to_component_index(i).first; + // Only want the velocity components and not the pressure one (which is the last one), so add 1 + if (introspection.is_stokes_component(index_direction+1)) + data.local_rhs(i) += ( + // RHS of momentum eqn: - \int 2/3 eta R, div v + - 2.0 / 3.0 * eta + * prescribed_dilation->dilation_rhs_term[index_direction][q] + * scratch.div_phi_u[i] + ) * JxW; + } } // and then the matrix, if necessary @@ -665,7 +686,10 @@ namespace aspect Assert(!this->get_parameters().enable_prescribed_dilation || (outputs.template get_additional_output_object>()->dilation_lhs_term.size() == n_points && - outputs.template get_additional_output_object>()->dilation_rhs_term.size() == n_points), + outputs.template get_additional_output_object>()->dilation_rhs_term.size() == dim && + outputs.template get_additional_output_object>()->dilation_rhs_term[0].size() == n_points && + outputs.template get_additional_output_object>()->dilation_rhs_term[1].size() == n_points && + (dim == 2 || (dim == 3 && outputs.template get_additional_output_object>()->dilation_rhs_term[2].size() == n_points))), ExcInternalError()); if (this->get_newton_handler().parameters.newton_derivative_scaling_factor != 0) diff --git a/source/simulator/assemblers/stokes.cc b/source/simulator/assemblers/stokes.cc index 7fd884b0394..7ab92d12ace 100644 --- a/source/simulator/assemblers/stokes.cc +++ b/source/simulator/assemblers/stokes.cc @@ -451,11 +451,32 @@ namespace aspect * JxW; if (prescribed_dilation != nullptr) - data.local_rhs(i) += ( - - pressure_scaling - * prescribed_dilation->dilation_rhs_term[q] - * scratch.phi_p[i] - ) * JxW; + { + const unsigned int index_direction=fe.system_to_component_index(i).first; + // Only want the velocity components and not the pressure one (which is the last one), so add 1 + if (introspection.is_stokes_component(index_direction+1)) + data.local_rhs(i) += ( + // RHS of - (div u,q) = - (R,q) + - pressure_scaling + * prescribed_dilation->dilation_rhs_term[index_direction][q] + * scratch.phi_p[i] + ) * JxW; + } + + // Only assemble this term if we are running incompressible, otherwise this term + // is already included on the LHS of the equation. + if (prescribed_dilation != nullptr && !material_model_is_compressible) + { + const unsigned int index_direction=fe.system_to_component_index(i).first; + // Only want the velocity components and not the pressure one (which is the last one), so add 1 + if (introspection.is_stokes_component(index_direction+1)) + data.local_rhs(i) += ( + // RHS of momentum eqn: - \int 2/3 eta R, div v + - 2.0 / 3.0 * eta + * prescribed_dilation->dilation_rhs_term[index_direction][q] + * scratch.div_phi_u[i] + ) * JxW; + } if (scratch.rebuild_stokes_matrix) for (unsigned int j=0; jget_parameters().enable_prescribed_dilation || (outputs.template get_additional_output_object>()->dilation_lhs_term.size() == n_points && - outputs.template get_additional_output_object>()->dilation_rhs_term.size() == n_points), + outputs.template get_additional_output_object>()->dilation_rhs_term.size() == dim && + outputs.template get_additional_output_object>()->dilation_rhs_term[0].size() == n_points && + outputs.template get_additional_output_object>()->dilation_rhs_term[1].size() == n_points && + (dim == 2 || (dim == 3 && outputs.template get_additional_output_object>()->dilation_rhs_term[2].size() == n_points))), ExcInternalError()); // Elasticity: diff --git a/tests/prescribed_dilation.cc b/tests/prescribed_dilation.cc index cb4afcf89b8..2a1be20a0c8 100644 --- a/tests/prescribed_dilation.cc +++ b/tests/prescribed_dilation.cc @@ -208,7 +208,10 @@ namespace aspect if (prescribed_dilation) { prescribed_dilation->dilation_rhs_term[i] = x; - prescribed_dilation->dilation_lhs_term[i] = 0.; + prescribed_dilation->dilation_lhs_term[0][i] = 0.; + prescribed_dilation->dilation_lhs_term[1][i] = 0.; + if (dim == 3) + prescribed_dilation->dilation_lhs_term[2][i] = 0.; } } From 771661a0fcfacc4db5b58937249d14f0d3cc165a Mon Sep 17 00:00:00 2001 From: Menno Fraters Date: Wed, 25 Jun 2025 12:50:20 +0200 Subject: [PATCH 2/8] fixup --- source/material_model/interface.cc | 6 +++--- source/material_model/visco_plastic.cc | 10 ++++++++-- source/simulator/assemblers/newton_stokes.cc | 3 ++- source/simulator/assemblers/stokes.cc | 1 + 4 files changed, 14 insertions(+), 6 deletions(-) diff --git a/source/material_model/interface.cc b/source/material_model/interface.cc index 7971a76717a..4ef6eb0bbb1 100644 --- a/source/material_model/interface.cc +++ b/source/material_model/interface.cc @@ -999,7 +999,7 @@ namespace aspect namespace { - std::vector make_prescribed_dilation_outputs_names() + std::vector make_prescribed_dilation_outputs_names(int dim) { std::vector names; names.emplace_back("dilation_lhs_term"); @@ -1015,9 +1015,9 @@ namespace aspect template PrescribedPlasticDilation::PrescribedPlasticDilation (const unsigned int n_points) - : NamedAdditionalMaterialOutputs(make_prescribed_dilation_outputs_names()), + : NamedAdditionalMaterialOutputs(make_prescribed_dilation_outputs_names(dim)), dilation_lhs_term(n_points, numbers::signaling_nan()), - dilation(dim, std::vector(n_points, numbers::signaling_nan())) + dilation_rhs_term(dim, std::vector(n_points, numbers::signaling_nan())) {} diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc index 0e54fad8dee..fd267df8ef7 100644 --- a/source/material_model/visco_plastic.cc +++ b/source/material_model/visco_plastic.cc @@ -299,9 +299,15 @@ namespace aspect // should cancel out (RHS = LHS * p) so as to satisfy the loading-unloading conditions. plastic_dilation->dilation_lhs_term[i] = dilation_lhs_term; if (dilation_rhs_term - dilation_lhs_term * in.pressure[i] > 0) - plastic_dilation->dilation_rhs_term[i] = dilation_rhs_term; + { + for (unsigned int dim_i = 0; dim_i < dim; ++dim_i) + plastic_dilation->dilation_rhs_term[dim_i][i] = dilation_rhs_term; + } else - plastic_dilation->dilation_rhs_term[i] = dilation_lhs_term * in.pressure[i]; + { + for (unsigned int dim_i = 0; dim_i < dim; ++dim_i) + plastic_dilation->dilation_rhs_term[dim_i][i] = dilation_lhs_term * in.pressure[i]; + } } } diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc index 0e2870575f7..c24030a4eec 100644 --- a/source/simulator/assemblers/newton_stokes.cc +++ b/source/simulator/assemblers/newton_stokes.cc @@ -436,6 +436,7 @@ namespace aspect const double JxW = scratch.finite_element_values.JxW(q); const double pressure_scaling = this->get_pressure_scaling(); + bool material_model_is_compressible = (this->get_material_model().is_compressible()); // first assemble the rhs for (unsigned int i=0; idilation_rhs_term[index_direction][q] - prescribed_dilation->dilation_lhs_term[q] * - scrasch.material_model_inputs.pressure[q]) + scratch.material_model_inputs.pressure[q]) * scratch.phi_p[i] ) * JxW; } diff --git a/source/simulator/assemblers/stokes.cc b/source/simulator/assemblers/stokes.cc index 7ab92d12ace..934d00b1f67 100644 --- a/source/simulator/assemblers/stokes.cc +++ b/source/simulator/assemblers/stokes.cc @@ -436,6 +436,7 @@ namespace aspect const double density = scratch.material_model_outputs.densities[q]; const double JxW = scratch.finite_element_values.JxW(q); + bool material_model_is_compressible = (this->get_material_model().is_compressible()); for (unsigned int i=0; i Date: Mon, 28 Jul 2025 12:46:18 +0200 Subject: [PATCH 3/8] Fixes to assemblers and test. --- source/simulator/assemblers/newton_stokes.cc | 32 +++++--------------- source/simulator/assemblers/stokes.cc | 16 ++++++---- tests/prescribed_dilation.cc | 8 ++--- 3 files changed, 21 insertions(+), 35 deletions(-) diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc index c24030a4eec..f035d5efeaf 100644 --- a/source/simulator/assemblers/newton_stokes.cc +++ b/source/simulator/assemblers/newton_stokes.cc @@ -464,31 +464,13 @@ namespace aspect if (enable_prescribed_dilation) { const unsigned int index_direction=fe.system_to_component_index(i).first; - // Only want the velocity components and not the pressure one (which is the last one), so add 1 - if (introspection.is_stokes_component(index_direction+1)) - data.local_rhs(i) += ( - // RHS of - (div u,q) = - (R,q) - - pressure_scaling - * (prescribed_dilation->dilation_rhs_term[index_direction][q] - - prescribed_dilation->dilation_lhs_term[q] * - scratch.material_model_inputs.pressure[q]) - * scratch.phi_p[i] - ) * JxW; - } - - // Only assemble this term if we are running incompressible, otherwise this term - // is already included on the LHS of the equation. - if (enable_prescribed_dilation && !material_model_is_compressible) - { - const unsigned int index_direction=fe.system_to_component_index(i).first; - // Only want the velocity components and not the pressure one (which is the last one), so add 1 - if (introspection.is_stokes_component(index_direction+1)) - data.local_rhs(i) += ( - // RHS of momentum eqn: - \int 2/3 eta R, div v - - 2.0 / 3.0 * eta - * prescribed_dilation->dilation_rhs_term[index_direction][q] - * scratch.div_phi_u[i] - ) * JxW; + data.local_rhs(i) += ( + - pressure_scaling + * (prescribed_dilation->dilation_rhs_term[index_direction][q] - + prescribed_dilation->dilation_lhs_term[q] * + scratch.material_model_inputs.pressure[q]) + * scratch.phi_p[i] + ) * JxW; } } diff --git a/source/simulator/assemblers/stokes.cc b/source/simulator/assemblers/stokes.cc index 934d00b1f67..9f09f2879e5 100644 --- a/source/simulator/assemblers/stokes.cc +++ b/source/simulator/assemblers/stokes.cc @@ -456,12 +456,16 @@ namespace aspect const unsigned int index_direction=fe.system_to_component_index(i).first; // Only want the velocity components and not the pressure one (which is the last one), so add 1 if (introspection.is_stokes_component(index_direction+1)) - data.local_rhs(i) += ( - // RHS of - (div u,q) = - (R,q) - - pressure_scaling - * prescribed_dilation->dilation_rhs_term[index_direction][q] - * scratch.phi_p[i] - ) * JxW; + { + const unsigned int index_direction=fe.system_to_component_index(i).first; + data.local_rhs(i) += ( + - pressure_scaling + * (prescribed_dilation->dilation_rhs_term[index_direction][q] - + prescribed_dilation->dilation_lhs_term[q] * + scratch.material_model_inputs.pressure[q]) + * scratch.phi_p[i] + ) * JxW; + } } // Only assemble this term if we are running incompressible, otherwise this term diff --git a/tests/prescribed_dilation.cc b/tests/prescribed_dilation.cc index 2a1be20a0c8..6679c19033a 100644 --- a/tests/prescribed_dilation.cc +++ b/tests/prescribed_dilation.cc @@ -207,11 +207,11 @@ namespace aspect } if (prescribed_dilation) { - prescribed_dilation->dilation_rhs_term[i] = x; - prescribed_dilation->dilation_lhs_term[0][i] = 0.; - prescribed_dilation->dilation_lhs_term[1][i] = 0.; + prescribed_dilation->dilation_lhs_term[i] = 0; + prescribed_dilation->dilation_rhs_term[0][i] = x; + prescribed_dilation->dilation_rhs_term[1][i] = x; if (dim == 3) - prescribed_dilation->dilation_lhs_term[2][i] = 0.; + prescribed_dilation->dilation_rhs_term[2][i] = x; } } From 9b868dae5c8a2416c6aafdc3b2cd441a3aaf848e Mon Sep 17 00:00:00 2001 From: Menno Fraters Date: Mon, 28 Jul 2025 13:29:19 +0200 Subject: [PATCH 4/8] more fixes. --- source/simulator/assemblers/newton_stokes.cc | 15 +++++----- source/simulator/assemblers/stokes.cc | 30 ++++++-------------- 2 files changed, 16 insertions(+), 29 deletions(-) diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc index f035d5efeaf..aa3b99acbeb 100644 --- a/source/simulator/assemblers/newton_stokes.cc +++ b/source/simulator/assemblers/newton_stokes.cc @@ -464,13 +464,14 @@ namespace aspect if (enable_prescribed_dilation) { const unsigned int index_direction=fe.system_to_component_index(i).first; - data.local_rhs(i) += ( - - pressure_scaling - * (prescribed_dilation->dilation_rhs_term[index_direction][q] - - prescribed_dilation->dilation_lhs_term[q] * - scratch.material_model_inputs.pressure[q]) - * scratch.phi_p[i] - ) * JxW; + if (index_direction < dim) + data.local_rhs(i) += ( + - pressure_scaling + * (prescribed_dilation->dilation_rhs_term[index_direction][q] - + prescribed_dilation->dilation_lhs_term[q] * + scratch.material_model_inputs.pressure[q]) + * scratch.phi_p[i] + ) * JxW; } } diff --git a/source/simulator/assemblers/stokes.cc b/source/simulator/assemblers/stokes.cc index 9f09f2879e5..7d29fa8cc4c 100644 --- a/source/simulator/assemblers/stokes.cc +++ b/source/simulator/assemblers/stokes.cc @@ -458,31 +458,17 @@ namespace aspect if (introspection.is_stokes_component(index_direction+1)) { const unsigned int index_direction=fe.system_to_component_index(i).first; - data.local_rhs(i) += ( - - pressure_scaling - * (prescribed_dilation->dilation_rhs_term[index_direction][q] - - prescribed_dilation->dilation_lhs_term[q] * - scratch.material_model_inputs.pressure[q]) - * scratch.phi_p[i] - ) * JxW; + if (index_direction < dim) + data.local_rhs(i) += ( + - pressure_scaling + * (prescribed_dilation->dilation_rhs_term[index_direction][q] - + prescribed_dilation->dilation_lhs_term[q] * + scratch.material_model_inputs.pressure[q]) + * scratch.phi_p[i] + ) * JxW; } } - // Only assemble this term if we are running incompressible, otherwise this term - // is already included on the LHS of the equation. - if (prescribed_dilation != nullptr && !material_model_is_compressible) - { - const unsigned int index_direction=fe.system_to_component_index(i).first; - // Only want the velocity components and not the pressure one (which is the last one), so add 1 - if (introspection.is_stokes_component(index_direction+1)) - data.local_rhs(i) += ( - // RHS of momentum eqn: - \int 2/3 eta R, div v - - 2.0 / 3.0 * eta - * prescribed_dilation->dilation_rhs_term[index_direction][q] - * scratch.div_phi_u[i] - ) * JxW; - } - if (scratch.rebuild_stokes_matrix) for (unsigned int j=0; j Date: Mon, 28 Jul 2025 17:09:27 +0200 Subject: [PATCH 5/8] Remove compressible variable. --- source/simulator/assemblers/newton_stokes.cc | 1 - source/simulator/assemblers/stokes.cc | 1 - 2 files changed, 2 deletions(-) diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc index aa3b99acbeb..c0445427228 100644 --- a/source/simulator/assemblers/newton_stokes.cc +++ b/source/simulator/assemblers/newton_stokes.cc @@ -436,7 +436,6 @@ namespace aspect const double JxW = scratch.finite_element_values.JxW(q); const double pressure_scaling = this->get_pressure_scaling(); - bool material_model_is_compressible = (this->get_material_model().is_compressible()); // first assemble the rhs for (unsigned int i=0; iget_material_model().is_compressible()); for (unsigned int i=0; i Date: Wed, 30 Jul 2025 15:33:14 +0200 Subject: [PATCH 6/8] Change order of vectors. --- source/simulator/assemblers/newton_stokes.cc | 10 ++++------ source/simulator/assemblers/stokes.cc | 1 - 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc index c0445427228..56aa15b2a6a 100644 --- a/source/simulator/assemblers/newton_stokes.cc +++ b/source/simulator/assemblers/newton_stokes.cc @@ -463,10 +463,10 @@ namespace aspect if (enable_prescribed_dilation) { const unsigned int index_direction=fe.system_to_component_index(i).first; - if (index_direction < dim) + if (introspection.is_stokes_component(index_direction+1)) data.local_rhs(i) += ( - pressure_scaling - * (prescribed_dilation->dilation_rhs_term[index_direction][q] - + * (prescribed_dilation->dilation_rhs_term[q][index_direction] - prescribed_dilation->dilation_lhs_term[q] * scratch.material_model_inputs.pressure[q]) * scratch.phi_p[i] @@ -669,10 +669,8 @@ namespace aspect Assert(!this->get_parameters().enable_prescribed_dilation || (outputs.template get_additional_output_object>()->dilation_lhs_term.size() == n_points && - outputs.template get_additional_output_object>()->dilation_rhs_term.size() == dim && - outputs.template get_additional_output_object>()->dilation_rhs_term[0].size() == n_points && - outputs.template get_additional_output_object>()->dilation_rhs_term[1].size() == n_points && - (dim == 2 || (dim == 3 && outputs.template get_additional_output_object>()->dilation_rhs_term[2].size() == n_points))), + outputs.template get_additional_output_object>()->dilation_rhs_term.size() == n_points + ), ExcInternalError()); if (this->get_newton_handler().parameters.newton_derivative_scaling_factor != 0) diff --git a/source/simulator/assemblers/stokes.cc b/source/simulator/assemblers/stokes.cc index c9f447a1232..04447264ff3 100644 --- a/source/simulator/assemblers/stokes.cc +++ b/source/simulator/assemblers/stokes.cc @@ -456,7 +456,6 @@ namespace aspect // Only want the velocity components and not the pressure one (which is the last one), so add 1 if (introspection.is_stokes_component(index_direction+1)) { - const unsigned int index_direction=fe.system_to_component_index(i).first; if (index_direction < dim) data.local_rhs(i) += ( - pressure_scaling From 2513e7f8ff925473ea0713bc8a7b3dd5bcf04aeb Mon Sep 17 00:00:00 2001 From: Menno Fraters Date: Wed, 30 Jul 2025 15:36:02 +0200 Subject: [PATCH 7/8] fixup --- source/simulator/assemblers/stokes.cc | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/source/simulator/assemblers/stokes.cc b/source/simulator/assemblers/stokes.cc index 04447264ff3..fd75f07db6c 100644 --- a/source/simulator/assemblers/stokes.cc +++ b/source/simulator/assemblers/stokes.cc @@ -456,14 +456,13 @@ namespace aspect // Only want the velocity components and not the pressure one (which is the last one), so add 1 if (introspection.is_stokes_component(index_direction+1)) { - if (index_direction < dim) - data.local_rhs(i) += ( - - pressure_scaling - * (prescribed_dilation->dilation_rhs_term[index_direction][q] - - prescribed_dilation->dilation_lhs_term[q] * - scratch.material_model_inputs.pressure[q]) - * scratch.phi_p[i] - ) * JxW; + data.local_rhs(i) += ( + - pressure_scaling + * (prescribed_dilation->dilation_rhs_term[index_direction][q] - + prescribed_dilation->dilation_lhs_term[q] * + scratch.material_model_inputs.pressure[q]) + * scratch.phi_p[i] + ) * JxW; } } From c4eae98e45087cc8f3834c5d348f9de16a3701dd Mon Sep 17 00:00:00 2001 From: Menno Fraters Date: Wed, 30 Jul 2025 15:56:47 +0200 Subject: [PATCH 8/8] revert some changes. --- source/simulator/assemblers/newton_stokes.cc | 8 +++++--- tests/prescribed_dilation.cc | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc index 56aa15b2a6a..489e1e91a27 100644 --- a/source/simulator/assemblers/newton_stokes.cc +++ b/source/simulator/assemblers/newton_stokes.cc @@ -466,7 +466,7 @@ namespace aspect if (introspection.is_stokes_component(index_direction+1)) data.local_rhs(i) += ( - pressure_scaling - * (prescribed_dilation->dilation_rhs_term[q][index_direction] - + * (prescribed_dilation->dilation_rhs_term[index_direction][q] - prescribed_dilation->dilation_lhs_term[q] * scratch.material_model_inputs.pressure[q]) * scratch.phi_p[i] @@ -669,8 +669,10 @@ namespace aspect Assert(!this->get_parameters().enable_prescribed_dilation || (outputs.template get_additional_output_object>()->dilation_lhs_term.size() == n_points && - outputs.template get_additional_output_object>()->dilation_rhs_term.size() == n_points - ), + outputs.template get_additional_output_object>()->dilation_rhs_term.size() == dim && + outputs.template get_additional_output_object>()->dilation_rhs_term[0].size() == n_points && + outputs.template get_additional_output_object>()->dilation_rhs_term[1].size() == n_points && + (dim == 2 || (dim == 3 && outputs.template get_additional_output_object>()->dilation_rhs_term[2].size() == n_points))), ExcInternalError()); if (this->get_newton_handler().parameters.newton_derivative_scaling_factor != 0) diff --git a/tests/prescribed_dilation.cc b/tests/prescribed_dilation.cc index 6679c19033a..465feede6ae 100644 --- a/tests/prescribed_dilation.cc +++ b/tests/prescribed_dilation.cc @@ -207,11 +207,11 @@ namespace aspect } if (prescribed_dilation) { - prescribed_dilation->dilation_lhs_term[i] = 0; prescribed_dilation->dilation_rhs_term[0][i] = x; prescribed_dilation->dilation_rhs_term[1][i] = x; if (dim == 3) prescribed_dilation->dilation_rhs_term[2][i] = x; + prescribed_dilation->dilation_lhs_term[i] = 0; } }