From 6463b7fa4b440e6c971f0629ea30eef1b12f0eab Mon Sep 17 00:00:00 2001 From: Yimin Jin Date: Sun, 30 Jun 2024 22:37:16 +0800 Subject: [PATCH 1/4] Add dilation relevant terms to the Newton solver --- .../material_model/rheology/drucker_prager.h | 17 +++ .../material_model/rheology/visco_plastic.h | 9 +- include/aspect/material_model/visco_plastic.h | 4 +- include/aspect/newton.h | 10 ++ include/aspect/stokes_matrix_free.h | 23 ++- .../material_model/rheology/drucker_prager.cc | 50 ++++++- .../material_model/rheology/visco_plastic.cc | 136 ++++++++++++------ source/material_model/visco_plastic.cc | 15 +- source/simulator/assemblers/newton_stokes.cc | 30 ++++ source/simulator/core.cc | 4 +- source/simulator/newton.cc | 2 + source/simulator/stokes_matrix_free.cc | 82 ++++++++--- 12 files changed, 306 insertions(+), 76 deletions(-) diff --git a/include/aspect/material_model/rheology/drucker_prager.h b/include/aspect/material_model/rheology/drucker_prager.h index d787f070357..904aa86e160 100644 --- a/include/aspect/material_model/rheology/drucker_prager.h +++ b/include/aspect/material_model/rheology/drucker_prager.h @@ -45,6 +45,11 @@ namespace aspect */ double angle_internal_friction; + /** + * Dilation angle (psi) for the current composition and phase + */ + double angle_dilation; + /** * Cohesion for the current composition and phase */ @@ -138,9 +143,21 @@ namespace aspect compute_derivative (const double angle_internal_friction, const double effective_strain_rate) const; + /** + * Compute the plastic dilation term, which is defined as the product of + * the plastic multiplier and the derivative of plastic potential with + * respect to minus pressure. + */ + double + compute_plastic_dilation (const double angle_dilation, + const double non_yielding_viscosity, + const double effective_viscosity, + const double effective_strain_rate) const; + private: std::vector angles_internal_friction; + std::vector angles_dilation; std::vector cohesions; double max_yield_stress; diff --git a/include/aspect/material_model/rheology/visco_plastic.h b/include/aspect/material_model/rheology/visco_plastic.h index e4f6b171102..9879c0e2939 100644 --- a/include/aspect/material_model/rheology/visco_plastic.h +++ b/include/aspect/material_model/rheology/visco_plastic.h @@ -108,6 +108,13 @@ namespace aspect * The current cohesion. */ std::vector current_cohesions; + + /** + * The plastic dilation terms, defined as the product of the + * plastic multiplier and the derivative of plastic potential + * with respect to negative pressure. + */ + std::vector plastic_dilation; }; namespace Rheology @@ -149,7 +156,7 @@ namespace aspect */ void compute_viscosity_derivatives(const unsigned int point_index, const std::vector &volume_fractions, - const std::vector &composition_viscosities, + const IsostrainViscosities &isostrain_values, const MaterialModel::MaterialModelInputs &in, MaterialModel::MaterialModelOutputs &out, const std::vector &phase_function_values = std::vector(), diff --git a/include/aspect/material_model/visco_plastic.h b/include/aspect/material_model/visco_plastic.h index 9ba9bf9fa67..62637a82048 100644 --- a/include/aspect/material_model/visco_plastic.h +++ b/include/aspect/material_model/visco_plastic.h @@ -23,7 +23,7 @@ #include #include -#include +#include #include #include @@ -268,7 +268,7 @@ namespace aspect /** * Object for computing the equation of state. */ - EquationOfState::MulticomponentIncompressible equation_of_state; + EquationOfState::MulticomponentCompressible equation_of_state; /** * Object that handles phase transitions. diff --git a/include/aspect/newton.h b/include/aspect/newton.h index f7e4fc70c46..591c24a73f8 100644 --- a/include/aspect/newton.h +++ b/include/aspect/newton.h @@ -63,6 +63,16 @@ namespace aspect * derivatives when material averaging is applied. */ std::vector viscosity_derivative_averaging_weights; + + /** + * The derivatives of the prescribed dilation with respect to strain rate. + */ + std::vector> dilation_derivative_wrt_strain_rate; + + /** + * The derivatives of the prescribed dilation term with respect to pressure. + */ + std::vector dilation_derivative_wrt_pressure; }; } diff --git a/include/aspect/stokes_matrix_free.h b/include/aspect/stokes_matrix_free.h index e90f857cc4a..8e2be8c0380 100644 --- a/include/aspect/stokes_matrix_free.h +++ b/include/aspect/stokes_matrix_free.h @@ -125,6 +125,12 @@ namespace aspect */ bool apply_stabilization_free_surface_faces; + /** + * If true, derivatives of the prescribed dilation are + * part of the operator. + */ + bool enable_prescribed_dilation; + /** * Table which stores viscosity values for each cell. * @@ -144,7 +150,7 @@ namespace aspect * variables: viscosity derivative with respect to pressure, * the Newton derivative scaling factor, and the averaging weight. */ - Table<2, VectorizedArray> newton_factor_wrt_pressure_table; + Table<2, VectorizedArray> viscosity_derivative_wrt_pressure_table; /** * Table which stores the product of the following four @@ -154,7 +160,20 @@ namespace aspect * is PD or SPD, otherwise, it is 1. */ Table<2, SymmetricTensor<2, dim, VectorizedArray>> - newton_factor_wrt_strain_rate_table; + viscosity_derivative_wrt_strain_rate_table; + + /** + * Table which stores the product of the dilation derivative + * with respect to pressure and the Newton derivative scaling factor. + */ + Table<2, VectorizedArray> dilation_derivative_wrt_pressure_table; + + /** + * Table which stores the product of the dilation derivative + * with respect to strain rate and the Newton derivative scaling factor. + */ + Table<2, SymmetricTensor<2, dim, VectorizedArray>> + dilation_derivative_wrt_strain_rate_table; /** * Table which stores the product of the pressure perturbation diff --git a/source/material_model/rheology/drucker_prager.cc b/source/material_model/rheology/drucker_prager.cc index aa5856a7a5e..67436cc4a8e 100644 --- a/source/material_model/rheology/drucker_prager.cc +++ b/source/material_model/rheology/drucker_prager.cc @@ -34,6 +34,7 @@ namespace aspect { DruckerPragerParameters::DruckerPragerParameters() : angle_internal_friction (numbers::signaling_nan()), + angle_dilation (numbers::signaling_nan()), cohesion (numbers::signaling_nan()), max_yield_stress (numbers::signaling_nan()) {} @@ -60,6 +61,7 @@ namespace aspect { // no phases drucker_prager_parameters.angle_internal_friction = angles_internal_friction[composition]; + drucker_prager_parameters.angle_dilation = angles_dilation[composition]; drucker_prager_parameters.cohesion = cohesions[composition]; } else @@ -67,6 +69,8 @@ namespace aspect // Average among phases drucker_prager_parameters.angle_internal_friction = MaterialModel::MaterialUtilities::phase_average_value(phase_function_values, n_phase_transitions_per_composition, angles_internal_friction, composition); + drucker_prager_parameters.angle_dilation = MaterialModel::MaterialUtilities::phase_average_value(phase_function_values, n_phase_transitions_per_composition, + angles_dilation, composition); drucker_prager_parameters.cohesion = MaterialModel::MaterialUtilities::phase_average_value(phase_function_values, n_phase_transitions_per_composition, cohesions, composition); } @@ -181,6 +185,26 @@ namespace aspect + template + double + DruckerPrager::compute_plastic_dilation(const double angle_dilation, + const double non_yielding_viscosity, + const double effective_viscosity, + const double effective_strain_rate) const + { + const double sin_psi = std::sin(angle_dilation); + const double minus_dQ_dp = (dim == 3 ? + 6.0 * sin_psi / (std::sqrt(3.0) * (3.0 + sin_psi)) : + sin_psi); + + const double plastic_multiplier = + std::max(0.0, (1.0 - effective_viscosity / non_yielding_viscosity) * 2.0 * effective_strain_rate); + + return minus_dQ_dp * plastic_multiplier; + } + + + template void DruckerPrager::declare_parameters (ParameterHandler &prm) @@ -192,6 +216,13 @@ namespace aspect "those corresponding to chemical compositions. " "For a value of zero, in 2d the von Mises criterion is retrieved. " "Angles higher than 30 degrees are harder to solve numerically. Units: degrees."); + prm.declare_entry ("Angles of dilation", "0.", + Patterns::Anything(), + "List of angles of plastic dilation, $\\psi$, for background material and compositional fields, " + "for a total of N+1 values, where N is the number of all compositional fields or only " + "those corresponding to chemical compositions. " + "For a value of zero, the von Mises flow rule is retrieved. " + "The dilation angle should never exceed the internal friction angle."); prm.declare_entry ("Cohesions", "1e20", Patterns::Anything(), "List of cohesions, $C$, for background material and compositional fields, " @@ -248,10 +279,25 @@ namespace aspect angles_internal_friction = Utilities::MapParsing::parse_map_to_double_array(prm.get("Angles of internal friction"), options); + angles_dilation = Utilities::MapParsing::parse_map_to_double_array(prm.get("Angles of dilation"), + options); // Convert angles from degrees to radians - for (double &angle : angles_internal_friction) - angle *= constants::degree_to_radians; + for (unsigned int i = 0; i < angles_internal_friction.size(); ++i) + { + angles_internal_friction[i] *= constants::degree_to_radians; + angles_dilation[i] *= constants::degree_to_radians; + + AssertThrow(angles_dilation[i] <= angles_internal_friction[i], + ExcMessage("The dilation angle is greater than the internal friction angle for " + "composition " + Utilities::to_string(i) + ".")); + + if (angles_dilation[i] > 0.0) + AssertThrow(this->get_parameters().enable_prescribed_dilation == true, + ExcMessage("ASPECT detected a nonzero dilation angle, but dilation is " + "not enabled. Please set parameter entry 'Enable prescribed " + "dilation' to 'true'.")); + } options.property_name = "Cohesions"; cohesions = Utilities::MapParsing::parse_map_to_double_array(prm.get("Cohesions"), diff --git a/source/material_model/rheology/visco_plastic.cc b/source/material_model/rheology/visco_plastic.cc index 2fe6259e7d2..55fe2cab300 100644 --- a/source/material_model/rheology/visco_plastic.cc +++ b/source/material_model/rheology/visco_plastic.cc @@ -109,6 +109,7 @@ namespace aspect output_parameters.composition_viscosities.resize(volume_fractions.size(), numbers::signaling_nan()); output_parameters.current_friction_angles.resize(volume_fractions.size(), numbers::signaling_nan()); output_parameters.current_cohesions.resize(volume_fractions.size(), numbers::signaling_nan()); + output_parameters.plastic_dilation.resize(volume_fractions.size(), numbers::signaling_nan()); // Assemble stress tensor if elastic behavior is enabled SymmetricTensor<2,dim> stress_old = numbers::signaling_nan>(); @@ -383,7 +384,23 @@ namespace aspect MaterialModel::MaterialUtilities::PhaseUtilities::logarithmic ); output_parameters.composition_viscosities[j] = std::min(std::max(effective_viscosity, minimum_viscosity_for_composition), maximum_viscosity_for_composition); + + // Compute the plastic dilation if necessary. + if (this->get_parameters().enable_prescribed_dilation == true) + { + if (output_parameters.composition_yielding[j] == true) + { + const double current_dilation = drucker_prager_parameters.angle_dilation * weakening_factors[1]; + output_parameters.plastic_dilation[j] = drucker_prager_plasticity.compute_plastic_dilation (current_dilation, + non_yielding_viscosity, + output_parameters.composition_viscosities[j], + effective_edot_ii); + } + else + output_parameters.plastic_dilation[j] = 0; + } } + return output_parameters; } @@ -394,7 +411,7 @@ namespace aspect ViscoPlastic:: compute_viscosity_derivatives(const unsigned int i, const std::vector &volume_fractions, - const std::vector &composition_viscosities, + const IsostrainViscosities ¤t_isostrain_values, const MaterialModel::MaterialModelInputs &in, MaterialModel::MaterialModelOutputs &out, const std::vector &phase_function_values, @@ -403,16 +420,22 @@ namespace aspect MaterialModel::MaterialModelDerivatives *derivatives = out.template get_additional_output>(); + const bool enable_dilation = this->get_parameters().enable_prescribed_dilation; + if (derivatives != nullptr) { // compute derivatives if necessary - std::vector> composition_viscosities_derivatives(volume_fractions.size()); - std::vector composition_dviscosities_dpressure(volume_fractions.size()); + const unsigned int n_compositions = volume_fractions.size(); + std::vector> composition_viscosity_derivatives_wrt_strain_rate(n_compositions); + std::vector composition_viscosity_derivatives_wrt_pressure(n_compositions); + + std::vector> composition_dilation_derivatives_wrt_strain_rate(enable_dilation ? n_compositions : 0); + std::vector composition_dilation_derivatives_wrt_pressure(enable_dilation ? n_compositions : 0); const double finite_difference_accuracy = 1e-7; // A new material model inputs variable that uses the strain rate and pressure difference. - MaterialModel::MaterialModelInputs in_derivatives = in; + MaterialModel::MaterialModelInputs in_forward = in; Assert(std::isfinite(in.strain_rate[i].norm()), ExcMessage("Invalid strain_rate in the MaterialModelInputs. This is likely because it was " @@ -425,59 +448,62 @@ namespace aspect { const TableIndices<2> strain_rate_indices = SymmetricTensor<2,dim>::unrolled_to_component_indices (component); - 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(component); + const double strain_rate_difference = std::max(std::fabs(deviatoric_strain_rate[strain_rate_indices]), min_strain_rate) + * finite_difference_accuracy; + const SymmetricTensor<2,dim> forward_strain_rate = deviatoric_strain_rate + + strain_rate_difference * Utilities::nth_basis_for_symmetric_tensors(component); - in_derivatives.strain_rate[i] = strain_rate_difference; + in_forward.strain_rate[i] = forward_strain_rate; - std::vector eta_component = - calculate_isostrain_viscosities(in_derivatives, i, volume_fractions, - phase_function_values, n_phase_transitions_per_composition).composition_viscosities; + const IsostrainViscosities forward_isostrain_values = + calculate_isostrain_viscosities(in_forward, i, volume_fractions, + phase_function_values, n_phase_transitions_per_composition); // For each composition of the independent component, compute the derivative. - for (unsigned int composition_index = 0; composition_index < eta_component.size(); ++composition_index) + for (unsigned int composition_index = 0; composition_index < n_compositions; ++composition_index) { - // compute the difference between the viscosity with and without the strain-rate difference. - double viscosity_derivative = eta_component[composition_index] - composition_viscosities[composition_index]; - if (viscosity_derivative != 0) - { - // when the difference is non-zero, divide by the difference. - viscosity_derivative /= std::max(std::fabs(strain_rate_difference[strain_rate_indices]), min_strain_rate) - * finite_difference_accuracy; - } - composition_viscosities_derivatives[composition_index][strain_rate_indices] = viscosity_derivative; + composition_viscosity_derivatives_wrt_strain_rate[composition_index][strain_rate_indices] + = ( forward_isostrain_values.composition_viscosities[composition_index] - + current_isostrain_values.composition_viscosities[composition_index] + ) / strain_rate_difference; + + if (enable_dilation) + composition_dilation_derivatives_wrt_strain_rate[composition_index][strain_rate_indices] + = ( forward_isostrain_values.plastic_dilation[composition_index] - + current_isostrain_values.plastic_dilation[composition_index] + ) / strain_rate_difference; } } // Now compute the derivative of the viscosity to the pressure - const double pressure_difference = in.pressure[i] + (std::fabs(in.pressure[i]) * finite_difference_accuracy); + const double pressure_difference = std::fabs(in.pressure[i]) * finite_difference_accuracy; + const double forward_pressure = in.pressure[i] + pressure_difference; - in_derivatives.pressure[i] = pressure_difference; + in_forward.pressure[i] = forward_pressure; // Modify the in_derivatives object again to take the original strain rate. - in_derivatives.strain_rate[i] = in.strain_rate[i]; + in_forward.strain_rate[i] = in.strain_rate[i]; - const std::vector viscosity_difference = - calculate_isostrain_viscosities(in_derivatives, i, volume_fractions, - phase_function_values, n_phase_transitions_per_composition).composition_viscosities; + const IsostrainViscosities forward_isostrain_values = + calculate_isostrain_viscosities(in_forward, i, volume_fractions, + phase_function_values, n_phase_transitions_per_composition); - for (unsigned int composition_index = 0; composition_index < viscosity_difference.size(); ++composition_index) + for (unsigned int composition_index = 0; composition_index < n_compositions; ++composition_index) { - double viscosity_derivative = viscosity_difference[composition_index] - composition_viscosities[composition_index]; - if (viscosity_difference[composition_index] != 0) - { - if (in.pressure[i] != 0) - { - viscosity_derivative /= std::fabs(in.pressure[i]) * finite_difference_accuracy; - } - else - { - viscosity_derivative = 0; - } - } - composition_dviscosities_dpressure[composition_index] = viscosity_derivative; + composition_viscosity_derivatives_wrt_pressure[composition_index] + = pressure_difference > 0.0 ? + ( forward_isostrain_values.composition_viscosities[composition_index] - + current_isostrain_values.composition_viscosities[composition_index] + ) / pressure_difference : + 0.0; + + if (enable_dilation) + composition_dilation_derivatives_wrt_pressure[composition_index] + = pressure_difference > 0.0 ? + ( forward_isostrain_values.plastic_dilation[composition_index] - + current_isostrain_values.plastic_dilation[composition_index] + ) / pressure_difference : + 0.0; } double viscosity_averaging_p = 0; // Geometric @@ -491,15 +517,33 @@ namespace aspect derivatives->viscosity_derivative_wrt_strain_rate[i] = Utilities::derivative_of_weighted_p_norm_average(out.viscosities[i], volume_fractions, - composition_viscosities, - composition_viscosities_derivatives, + current_isostrain_values.composition_viscosities, + composition_viscosity_derivatives_wrt_strain_rate, viscosity_averaging_p); derivatives->viscosity_derivative_wrt_pressure[i] = Utilities::derivative_of_weighted_p_norm_average(out.viscosities[i], volume_fractions, - composition_viscosities, - composition_dviscosities_dpressure, + current_isostrain_values.composition_viscosities, + composition_viscosity_derivatives_wrt_pressure, viscosity_averaging_p); + + if (enable_dilation) + { + double dummy_plastic_dilation = numbers::signaling_nan(); + // always use arithmetic average for plastic dilation terms + derivatives->dilation_derivative_wrt_strain_rate[i] = + Utilities::derivative_of_weighted_p_norm_average(dummy_plastic_dilation, + volume_fractions, + current_isostrain_values.plastic_dilation, + composition_dilation_derivatives_wrt_strain_rate, + 1); + derivatives->dilation_derivative_wrt_pressure[i] = + Utilities::derivative_of_weighted_p_norm_average(dummy_plastic_dilation, + volume_fractions, + current_isostrain_values.plastic_dilation, + composition_dilation_derivatives_wrt_pressure, + 1); + } } } diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc index cef81e360e1..c0a2b33638d 100644 --- a/source/material_model/visco_plastic.cc +++ b/source/material_model/visco_plastic.cc @@ -247,7 +247,7 @@ namespace aspect out.template get_additional_output>()) rheology->compute_viscosity_derivatives(i, volume_fractions, - isostrain_viscosities.composition_viscosities, + isostrain_viscosities, in, out, phase_function_values, n_phase_transitions_for_each_chemical_composition); } @@ -297,6 +297,14 @@ namespace aspect elastic_out->elastic_shear_moduli[i] = average_elastic_shear_moduli[i]; } } + + // Compute the arithmetic average of plastic dilation term if necessary + if (PrescribedPlasticDilation *plastic_dilation = + out.template get_additional_output>()) + plastic_dilation->dilation[i] + = MaterialUtilities::average_value (volume_fractions, + isostrain_viscosities.plastic_dilation, + MaterialUtilities::arithmetic); } // If we use the full strain tensor, compute the change in the individual tensor components. @@ -340,7 +348,7 @@ namespace aspect { MaterialUtilities::PhaseFunction::declare_parameters(prm); - EquationOfState::MulticomponentIncompressible::declare_parameters (prm); + EquationOfState::MulticomponentCompressible::declare_parameters (prm); Rheology::ViscoPlastic::declare_parameters(prm); @@ -402,8 +410,7 @@ namespace aspect // Equation of state parameters equation_of_state.initialize_simulator (this->get_simulator()); - equation_of_state.parse_parameters (prm, - std::make_unique>(n_phases_for_each_chemical_composition)); + equation_of_state.parse_parameters(prm); // Make options file for parsing maps to double arrays std::vector chemical_field_names = this->introspection().chemical_composition_field_names(); diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc index 174bab9cbc2..6a9c5cbab65 100644 --- a/source/simulator/assemblers/newton_stokes.cc +++ b/source/simulator/assemblers/newton_stokes.cc @@ -67,6 +67,8 @@ namespace aspect AssertThrow(elastic_out != nullptr, ExcMessage("Error: The Newton method requires ElasticOutputs when elasticity is enabled.")); + const bool enable_prescribed_dilation = this->get_parameters().enable_prescribed_dilation; + // First loop over all dofs and find those that are in the Stokes system // save the component (pressure and dim velocities) each belongs to. for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/) @@ -229,6 +231,17 @@ namespace aspect ) * JxW; } + + if (enable_prescribed_dilation) + { + for (unsigned int i=0; idilation_derivative_wrt_pressure[q] + * JxW; + + } } } #if DEBUG @@ -524,6 +537,23 @@ namespace aspect " = " + Utilities::to_string(eta))); } } + + if (enable_prescribed_dilation) + { + std::vector dilation_differentiations(stokes_dofs_per_cell); + for (unsigned int k=0; kdilation_derivative_wrt_strain_rate[q] * scratch.grads_phi_u[k] + + derivatives->dilation_derivative_wrt_pressure[q] * pressure_scaling * scratch.phi_p[k] + ) * derivative_scaling_factor; + + for (unsigned int i=0; i()) , viscosity_derivative_wrt_strain_rate(n_points, numbers::signaling_nan>()) , viscosity_derivative_averaging_weights(n_points, numbers::signaling_nan()) + , dilation_derivative_wrt_strain_rate(n_points, numbers::signaling_nan>()) + , dilation_derivative_wrt_pressure(n_points, numbers::signaling_nan()) {} } diff --git a/source/simulator/stokes_matrix_free.cc b/source/simulator/stokes_matrix_free.cc index 646e3d8437b..6858d601f8d 100644 --- a/source/simulator/stokes_matrix_free.cc +++ b/source/simulator/stokes_matrix_free.cc @@ -342,9 +342,11 @@ namespace aspect OperatorCellData::memory_consumption() const { return viscosity.memory_consumption() - + newton_factor_wrt_pressure_table.memory_consumption() + + viscosity_derivative_wrt_pressure_table.memory_consumption() + strain_rate_table.memory_consumption() - + newton_factor_wrt_strain_rate_table.memory_consumption(); + + viscosity_derivative_wrt_strain_rate_table.memory_consumption() + + dilation_derivative_wrt_pressure_table.memory_consumption() + + dilation_derivative_wrt_strain_rate_table.memory_consumption(); } @@ -355,9 +357,11 @@ namespace aspect { enable_newton_derivatives = false; viscosity.clear(); - newton_factor_wrt_pressure_table.clear(); + viscosity_derivative_wrt_pressure_table.clear(); strain_rate_table.clear(); - newton_factor_wrt_strain_rate_table.clear(); + viscosity_derivative_wrt_strain_rate_table.clear(); + dilation_derivative_wrt_pressure_table.clear(); + dilation_derivative_wrt_strain_rate_table.clear(); } } @@ -455,9 +459,19 @@ namespace aspect const VectorizedArray val_p_q = p_eval.get_value(q); // Terms to be tested by phi_p: - const VectorizedArray pressure_terms = + VectorizedArray pressure_terms = -cell_data->pressure_scaling * div_u_q; + if (cell_data->enable_newton_derivatives && + cell_data->enable_prescribed_dilation) + { + pressure_terms += ( cell_data->dilation_derivative_wrt_strain_rate_table(cell,q) + * sym_grad_u_q ) + + + ( cell_data->dilation_derivative_wrt_pressure_table(cell,q) + * cell_data->pressure_scaling * val_p_q ); + } + // Terms to be tested by the symmetric gradients of phi_u: SymmetricTensor<2,dim,VectorizedArray> velocity_terms = viscosity_x_2 * sym_grad_u_q; @@ -477,9 +491,9 @@ namespace aspect VectorizedArray deta_dp_times_p(0.); for (const unsigned int r : u_eval.quadrature_point_indices()) { - deta_deps_times_sym_grad_u += cell_data->newton_factor_wrt_strain_rate_table(cell,r) + deta_deps_times_sym_grad_u += cell_data->viscosity_derivative_wrt_strain_rate_table(cell,r) * sym_grad_u[r]; - deta_dp_times_p += cell_data->newton_factor_wrt_pressure_table(cell,r) * val_p[r]; + deta_dp_times_p += cell_data->viscosity_derivative_wrt_pressure_table(cell,r) * val_p[r]; if (cell_data->symmetrize_newton_system) eps_times_sym_grad_u += cell_data->strain_rate_table(cell,r) * sym_grad_u[r]; } @@ -487,7 +501,7 @@ namespace aspect velocity_terms += ( cell_data->symmetrize_newton_system ? ( cell_data->strain_rate_table(cell,q) * deta_deps_times_sym_grad_u + - cell_data->newton_factor_wrt_strain_rate_table(cell,q) * eps_times_sym_grad_u ) : + cell_data->viscosity_derivative_wrt_strain_rate_table(cell,q) * eps_times_sym_grad_u ) : 2. * cell_data->strain_rate_table(cell,q) * deta_deps_times_sym_grad_u ) + 2. * cell_data->strain_rate_table(cell,q) * deta_dp_times_p; @@ -1375,8 +1389,15 @@ namespace aspect const unsigned int n_q_points = quadrature_formula.size(); active_cell_data.strain_rate_table.reinit(TableIndices<2>(n_cells, n_q_points)); - active_cell_data.newton_factor_wrt_pressure_table.reinit(TableIndices<2>(n_cells, n_q_points)); - active_cell_data.newton_factor_wrt_strain_rate_table.reinit(TableIndices<2>(n_cells, n_q_points)); + active_cell_data.viscosity_derivative_wrt_pressure_table.reinit(TableIndices<2>(n_cells, n_q_points)); + active_cell_data.viscosity_derivative_wrt_strain_rate_table.reinit(TableIndices<2>(n_cells, n_q_points)); + + if (sim.parameters.enable_prescribed_dilation) + { + active_cell_data.enable_prescribed_dilation = true; + active_cell_data.dilation_derivative_wrt_pressure_table.reinit(TableIndices<2>(n_cells, n_q_points)); + active_cell_data.dilation_derivative_wrt_strain_rate_table.reinit(TableIndices<2>(n_cells, n_q_points)); + } for (unsigned int cell=0; cellviscosity_derivative_wrt_pressure[q] * derivatives->viscosity_derivative_averaging_weights[q] * newton_derivative_scaling_factor; - Assert(std::isfinite(active_cell_data.newton_factor_wrt_pressure_table(cell,q)[i]), - ExcMessage("active_cell_data.newton_factor_wrt_pressure_table is not finite: " + std::to_string(active_cell_data.newton_factor_wrt_pressure_table(cell,q)[i]) + + Assert(std::isfinite(active_cell_data.viscosity_derivative_wrt_pressure_table(cell,q)[i]), + ExcMessage("active_cell_data.viscosity_derivative_wrt_pressure_table is not finite: " + std::to_string(active_cell_data.viscosity_derivative_wrt_pressure_table(cell,q)[i]) + ". Relevant variables are derivatives->viscosity_derivative_wrt_pressure[q] = " + std::to_string(derivatives->viscosity_derivative_wrt_pressure[q]) + ", derivatives->viscosity_derivative_averaging_weights[q] = " + std::to_string(derivatives->viscosity_derivative_averaging_weights[q]) + ", and newton_derivative_scaling_factor = " + std::to_string(newton_derivative_scaling_factor))); @@ -1451,16 +1472,39 @@ namespace aspect active_cell_data.strain_rate_table(cell, q)[m][n][i] = effective_strain_rate[m][n]; - active_cell_data.newton_factor_wrt_strain_rate_table(cell, q)[m][n][i] + active_cell_data.viscosity_derivative_wrt_strain_rate_table(cell, q)[m][n][i] = derivatives->viscosity_derivative_wrt_strain_rate[q][m][n] * derivatives->viscosity_derivative_averaging_weights[q] * newton_derivative_scaling_factor * alpha; Assert(std::isfinite(active_cell_data.strain_rate_table(cell, q)[m][n][i]), ExcMessage("active_cell_data.strain_rate_table has an element which is not finite: " + std::to_string(active_cell_data.strain_rate_table(cell, q)[m][n][i]))); - Assert(std::isfinite(active_cell_data.newton_factor_wrt_strain_rate_table(cell, q)[m][n][i]), - ExcMessage("active_cell_data.newton_factor_wrt_strain_rate_table has an element which is not finite: " + std::to_string(active_cell_data.newton_factor_wrt_strain_rate_table(cell, q)[m][n][i]))); + Assert(std::isfinite(active_cell_data.viscosity_derivative_wrt_strain_rate_table(cell, q)[m][n][i]), + ExcMessage("active_cell_data.viscosity_derivative_wrt_strain_rate_table has an element which is not finite: " + std::to_string(active_cell_data.viscosity_derivative_wrt_strain_rate_table(cell, q)[m][n][i]))); } + + if (active_cell_data.enable_prescribed_dilation == true) + { + active_cell_data.dilation_derivative_wrt_pressure_table(cell,q)[i] + = derivatives->dilation_derivative_wrt_pressure[q] * + newton_derivative_scaling_factor; + Assert(std::isfinite(active_cell_data.dilation_derivative_wrt_pressure_table(cell,q)[i]), + ExcMessage("active_cell_data.dilation_derivative_wrt_pressure_table is not finite: " + std::to_string(active_cell_data.dilation_derivative_wrt_pressure_table(cell,q)[i]) + + ". Relevant variables are derivatives->dilation_derivative_wrt_pressure[q] = " + std::to_string(derivatives->dilation_derivative_wrt_pressure[q]) + + " and newton_derivative_scaling_factor = " + std::to_string(newton_derivative_scaling_factor))); + + for (unsigned int m=0; mdilation_derivative_wrt_strain_rate[q][m][n] * + newton_derivative_scaling_factor; + Assert(std::isfinite(active_cell_data.dilation_derivative_wrt_pressure_table(cell,q)[i]), + ExcMessage("active_cell_data.dilation_derivative_wrt_strain_rate_table is not finite: " + std::to_string(active_cell_data.dilation_derivative_wrt_pressure_table(cell,q)[i]) + + ". Relevant variables are derivatives->dilation_derivative_wrt_strain_rate[q] = " + std::to_string(derivatives->dilation_derivative_wrt_pressure[q]) + + " and newton_derivative_scaling_factor = " + std::to_string(newton_derivative_scaling_factor))); + } + } } } } @@ -1476,9 +1520,11 @@ namespace aspect // delete data used for Newton derivatives if necessary // TODO: use Table::clear() once implemented in 10.0.pre active_cell_data.enable_newton_derivatives = false; - active_cell_data.newton_factor_wrt_pressure_table.reinit(TableIndices<2>(0,0)); + active_cell_data.viscosity_derivative_wrt_pressure_table.reinit(TableIndices<2>(0,0)); active_cell_data.strain_rate_table.reinit(TableIndices<2>(0,0)); - active_cell_data.newton_factor_wrt_strain_rate_table.reinit(TableIndices<2>(0,0)); + active_cell_data.viscosity_derivative_wrt_strain_rate_table.reinit(TableIndices<2>(0,0)); + active_cell_data.dilation_derivative_wrt_pressure_table.reinit(TableIndices<2>(0,0)); + active_cell_data.dilation_derivative_wrt_strain_rate_table.reinit(TableIndices<2>(0,0)); for (unsigned int level=0; level Date: Thu, 25 Jul 2024 16:54:36 +0800 Subject: [PATCH 2/4] fix bugs: missing pressure_scaling in local_apply() --- source/simulator/stokes_matrix_free.cc | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/source/simulator/stokes_matrix_free.cc b/source/simulator/stokes_matrix_free.cc index 6858d601f8d..0d6da231d95 100644 --- a/source/simulator/stokes_matrix_free.cc +++ b/source/simulator/stokes_matrix_free.cc @@ -465,11 +465,14 @@ namespace aspect if (cell_data->enable_newton_derivatives && cell_data->enable_prescribed_dilation) { - pressure_terms += ( cell_data->dilation_derivative_wrt_strain_rate_table(cell,q) - * sym_grad_u_q ) - + - ( cell_data->dilation_derivative_wrt_pressure_table(cell,q) - * cell_data->pressure_scaling * val_p_q ); + pressure_terms += ( + ( cell_data->dilation_derivative_wrt_strain_rate_table(cell,q) + * sym_grad_u_q ) + + + ( cell_data->dilation_derivative_wrt_pressure_table(cell,q) + * cell_data->pressure_scaling * val_p_q ) + ) + * cell_data->pressure_scaling; } // Terms to be tested by the symmetric gradients of phi_u: @@ -504,7 +507,7 @@ namespace aspect cell_data->viscosity_derivative_wrt_strain_rate_table(cell,q) * eps_times_sym_grad_u ) : 2. * cell_data->strain_rate_table(cell,q) * deta_deps_times_sym_grad_u ) + - 2. * cell_data->strain_rate_table(cell,q) * deta_dp_times_p; + 2. * cell_data->strain_rate_table(cell,q) * deta_dp_times_p * cell_data->pressure_scaling; } u_eval.submit_symmetric_gradient(velocity_terms, q); From 251f1441a079a558e0037f4fd6e87fe9e99c15eb Mon Sep 17 00:00:00 2001 From: Yimin Jin Date: Fri, 26 Jul 2024 00:48:29 +0800 Subject: [PATCH 3/4] change the eos back to incompressible model --- include/aspect/material_model/visco_plastic.h | 4 ++-- source/material_model/visco_plastic.cc | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/include/aspect/material_model/visco_plastic.h b/include/aspect/material_model/visco_plastic.h index 62637a82048..9ba9bf9fa67 100644 --- a/include/aspect/material_model/visco_plastic.h +++ b/include/aspect/material_model/visco_plastic.h @@ -23,7 +23,7 @@ #include #include -#include +#include #include #include @@ -268,7 +268,7 @@ namespace aspect /** * Object for computing the equation of state. */ - EquationOfState::MulticomponentCompressible equation_of_state; + EquationOfState::MulticomponentIncompressible equation_of_state; /** * Object that handles phase transitions. diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc index c0a2b33638d..bb23225e649 100644 --- a/source/material_model/visco_plastic.cc +++ b/source/material_model/visco_plastic.cc @@ -348,7 +348,7 @@ namespace aspect { MaterialUtilities::PhaseFunction::declare_parameters(prm); - EquationOfState::MulticomponentCompressible::declare_parameters (prm); + EquationOfState::MulticomponentIncompressible::declare_parameters (prm); Rheology::ViscoPlastic::declare_parameters(prm); From b1c748eb309a5f553833267f80c891e9921fbc4e Mon Sep 17 00:00:00 2001 From: Yimin Jin Date: Fri, 26 Jul 2024 00:49:29 +0800 Subject: [PATCH 4/4] treat plastic dilation in the same way as compressible model --- source/simulator/assemblers/newton_stokes.cc | 10 ---------- source/simulator/assemblers/stokes.cc | 14 -------------- source/simulator/assembly.cc | 2 +- source/simulator/newton.cc | 3 ++- source/simulator/stokes_matrix_free.cc | 9 ++++++--- 5 files changed, 9 insertions(+), 29 deletions(-) diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc index 6a9c5cbab65..2f23efb9f86 100644 --- a/source/simulator/assemblers/newton_stokes.cc +++ b/source/simulator/assemblers/newton_stokes.cc @@ -451,16 +451,6 @@ namespace aspect * prescribed_dilation->dilation[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) - data.local_rhs(i) += ( - // RHS of momentum eqn: - \int 2/3 eta R, div v - - 2.0 / 3.0 * eta - * prescribed_dilation->dilation[q] - * scratch.div_phi_u[i] - ) * JxW; } // and then the matrix, if necessary diff --git a/source/simulator/assemblers/stokes.cc b/source/simulator/assemblers/stokes.cc index 164a7a7913e..935889fca30 100644 --- a/source/simulator/assemblers/stokes.cc +++ b/source/simulator/assemblers/stokes.cc @@ -374,10 +374,6 @@ namespace aspect { scratch.grads_phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].symmetric_gradient(i,q); } - else if (prescribed_dilation && !material_model_is_compressible) - { - scratch.div_phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].divergence (i, q); - } ++i_stokes; } ++i; @@ -426,16 +422,6 @@ namespace aspect * 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) - data.local_rhs(i) += ( - // RHS of momentum eqn: - \int 2/3 eta R, div v - - 2.0 / 3.0 * eta - * prescribed_dilation->dilation[q] - * scratch.div_phi_u[i] - ) * JxW; - if (scratch.rebuild_stokes_matrix) for (unsigned int j=0; jstokes_preconditioner.push_back(std::make_unique>()); assemblers->stokes_system.push_back(std::make_unique>()); - if (material_model->is_compressible()) + if (material_model->is_compressible() || parameters.enable_prescribed_dilation) { // The compressible part of the preconditioner is only necessary if we use the simplified A block if (parameters.use_full_A_block_preconditioner == false) diff --git a/source/simulator/newton.cc b/source/simulator/newton.cc index d0fea7e73e8..49d2ea4d76d 100644 --- a/source/simulator/newton.cc +++ b/source/simulator/newton.cc @@ -52,7 +52,8 @@ namespace aspect assemblers.stokes_preconditioner.push_back(std::make_unique>()); assemblers.stokes_system.push_back(std::make_unique>()); - if (this->get_material_model().is_compressible()) + if (this->get_material_model().is_compressible() || + this->get_parameters().enable_prescribed_dilation) { // The compressible part of the preconditioner is only necessary if we use the simplified A block if (this->get_parameters().use_full_A_block_preconditioner == false) diff --git a/source/simulator/stokes_matrix_free.cc b/source/simulator/stokes_matrix_free.cc index 0d6da231d95..ac693e83920 100644 --- a/source/simulator/stokes_matrix_free.cc +++ b/source/simulator/stokes_matrix_free.cc @@ -482,7 +482,8 @@ namespace aspect for (unsigned int d=0; dpressure_scaling * val_p_q; - if (cell_data->is_compressible) + if (cell_data->is_compressible || + cell_data->enable_prescribed_dilation) for (unsigned int d=0; dis_compressible) + if (cell_data->is_compressible || + cell_data->enable_prescribed_dilation) { const VectorizedArray div = trace(sym_grad_u); for (unsigned int d=0; dis_compressible(); + const bool is_compressible = sim.material_model->is_compressible() || + sim.parameters.enable_prescribed_dilation; dealii::LinearAlgebra::distributed::BlockVector rhs_correction(2); dealii::LinearAlgebra::distributed::BlockVector u0(2);