From 121ddde0551b4f2c5b61d954d4789cf0b784b102 Mon Sep 17 00:00:00 2001 From: j507 Date: Thu, 11 Jul 2024 17:37:38 +0200 Subject: [PATCH] (A) WIP of bouncing algorithm and clean-up of monolithic algorithm --- include/gCP/fe_field.h | 5 +- include/gCP/gradient_crystal_plasticity.h | 53 +- source/fe_field.cc | 78 +-- source/gradient_crystal_plasticity/solve.cc | 563 ++++++++++++-------- 4 files changed, 408 insertions(+), 291 deletions(-) diff --git a/include/gCP/fe_field.h b/include/gCP/fe_field.h index 1d4d6d9..9cc887e 100644 --- a/include/gCP/fe_field.h +++ b/include/gCP/fe_field.h @@ -309,7 +309,10 @@ dealii::types::global_dof_index get_n_plastic_slip_dofs() const; */ unsigned int get_n_components() const; -std::vector get_l2_norms( +double get_l2_norm( + const dealii::LinearAlgebraTrilinos::MPI::BlockVector &vector) const; + +dealii::Vector get_sub_l2_norms( const dealii::LinearAlgebraTrilinos::MPI::BlockVector &vector) const; /*! diff --git a/include/gCP/gradient_crystal_plasticity.h b/include/gCP/gradient_crystal_plasticity.h index d5e48a2..61ba6e0 100644 --- a/include/gCP/gradient_crystal_plasticity.h +++ b/include/gCP/gradient_crystal_plasticity.h @@ -57,7 +57,7 @@ class GradientCrystalPlasticitySolver void set_macroscopic_strain( const dealii::SymmetricTensor<2,dim> macroscopic_strain); - std::tuple solve_nonlinear_system( + void solve_nonlinear_system( const bool flag_skip_extrapolation = false); std::shared_ptr> @@ -89,28 +89,6 @@ class GradientCrystalPlasticitySolver private: - struct SolverData - { - SolverData() - : - flag_successful_convergence(false), - flag_compute_active_set(true), - nonlinear_iteration(0), - residual_l2_norms(3,0.), - old_residual_l2_norms(3,0.) - {}; - - bool flag_successful_convergence; - - bool flag_compute_active_set; - - unsigned int nonlinear_iteration; - - std::vector residual_l2_norms; - - std::vector old_residual_l2_norms; - }; - const RunTimeParameters::SolverParameters ¶meters; const dealii::DiscreteTime &discrete_time; @@ -248,18 +226,20 @@ class GradientCrystalPlasticitySolver void copy_local_to_global_quadrature_point_history( const gCP::AssemblyData::QuadraturePointHistory::Copy &){}; - void monolithic_algorithm(SolverData &solver_data); + void monolithic_algorithm(); - void bouncing_algorithm(SolverData &solver_data); + void bouncing_algorithm(); - void embracing_algorihtm(SolverData &solver_data); + void embracing_algorihtm(); unsigned int solve_linearized_system( const RunTimeParameters::KrylovParameters &krylov_parameters, const unsigned int block_id, const double right_hand_side_l2_norm); - void update_trial_solution(const double relaxation_parameter); + void update_trial_solution( + const double relaxation_parameter, + const unsigned int block_id = 0); void update_trial_solution(const std::vector relaxation_parameter); @@ -269,6 +249,7 @@ class GradientCrystalPlasticitySolver double line_search_algorithm( dealii::Vector residual_l2_norms, const unsigned int block_id); + void store_trial_solution( const bool flag_store_initial_trial_solution = false); @@ -278,17 +259,19 @@ class GradientCrystalPlasticitySolver void extrapolate_initial_trial_solution( const bool flag_skip_extrapolation = false); + void distribute_affine_constraints_to_trial_solution(); + void update_and_output_nonlinear_solver_logger( - const std::vector residual_l2_norms); + const dealii::Vector residual_l2_norms); void update_and_output_nonlinear_solver_logger( - const unsigned int nonlinear_iteration, - const unsigned int n_krylov_iterations, - const unsigned int n_line_search_iterations, - const std::vector newton_update_l2_norms, - const std::vector residual_l2_norms, - const double order_of_convergence, - const double relaxation_parameter = 1.0); + const unsigned int nonlinear_iteration, + const unsigned int n_krylov_iterations, + const unsigned int n_line_search_iterations, + const dealii::Vector newton_update_l2_norms, + const dealii::Vector residual_l2_norms, + const double order_of_convergence, + const double relaxation_parameter = 1.0); /*! * @note Only for debugging purposes diff --git a/source/fe_field.cc b/source/fe_field.cc index d58ea8c..cb6cc0e 100644 --- a/source/fe_field.cc +++ b/source/fe_field.cc @@ -584,7 +584,26 @@ void FEField::reset_all_affine_constraints() template -std::vector FEField::get_l2_norms( +double FEField::get_l2_norm( + const dealii::LinearAlgebraTrilinos::MPI::BlockVector &vector) const +{ + if (vector.has_ghost_elements()) + { + dealii::LinearAlgebraTrilinos::MPI::BlockVector tmp = + get_distributed_vector_instance(vector); + + return (tmp.l2_norm()); + } + else + { + return (vector.l2_norm()); + } +} + + + +template +dealii::Vector FEField::get_sub_l2_norms( const dealii::LinearAlgebraTrilinos::MPI::BlockVector &vector) const { dealii::LinearAlgebraTrilinos::MPI::BlockVector distributed_vector_tmp; @@ -593,9 +612,9 @@ std::vector FEField::get_l2_norms( distributed_vector_tmp = vector; - std::vector l2_norms(3, 0.); + dealii::Vector l2_norms(2); - l2_norms[0] = distributed_vector_tmp.l2_norm(); + const double l2_norm = distributed_vector_tmp.l2_norm(); if (flag_use_single_block) { @@ -636,42 +655,37 @@ std::vector FEField::get_l2_norms( scalar_squared_entries, MPI_COMM_WORLD); - const double control_l2_norm = - std::sqrt(vector_squared_entries + scalar_squared_entries); - - auto to_string = - [](const double number) - { - std::ostringstream out; - out.precision(10); - out << std::scientific << number; - return std::move(out).str(); - }; - - const double factor = 1e4; - - Assert( - std::fabs(l2_norms[0] - control_l2_norm) < - std::numeric_limits::epsilon() * factor, - dealii::ExcMessage("The norms do not match (" + - to_string(l2_norms[0]) + ", " + to_string(control_l2_norm) + - ")")); + l2_norms[0] = std::sqrt(vector_squared_entries); - (void)control_l2_norm; - - (void)to_string; - - l2_norms[1] = std::sqrt(vector_squared_entries); - - l2_norms[2] = std::sqrt(scalar_squared_entries); + l2_norms[1] = std::sqrt(scalar_squared_entries); } else { - l2_norms[1] = distributed_vector_tmp.block(0).l2_norm(); + l2_norms[0] = distributed_vector_tmp.block(0).l2_norm(); - l2_norms[2] = distributed_vector_tmp.block(1).l2_norm(); + l2_norms[1] = distributed_vector_tmp.block(1).l2_norm(); } + auto to_string = + [](const double number) + { + std::ostringstream out; + out.precision(10); + out << std::scientific << number; + return std::move(out).str(); + }; + + const double factor = 1e4; + + Assert( + std::fabs(l2_norm - l2_norms.l2_norm()) < + std::numeric_limits::epsilon() * factor, + dealii::ExcMessage("The norms do not match (" + + to_string(l2_norm) + ", " + to_string(l2_norms.l2_norm()) + + ")")); + + (void)to_string; + return l2_norms; } diff --git a/source/gradient_crystal_plasticity/solve.cc b/source/gradient_crystal_plasticity/solve.cc index 91ac97d..ad30ab7 100644 --- a/source/gradient_crystal_plasticity/solve.cc +++ b/source/gradient_crystal_plasticity/solve.cc @@ -63,20 +63,35 @@ namespace gCP template - std::tuple GradientCrystalPlasticitySolver:: + void GradientCrystalPlasticitySolver:: + distribute_affine_constraints_to_trial_solution() + { + dealii::LinearAlgebraTrilinos::MPI::BlockVector + distributed_trial_solution = + fe_field->get_distributed_vector_instance(trial_solution); + + fe_field->get_affine_constraints().distribute( + distributed_trial_solution); + + trial_solution = distributed_trial_solution; + } + + + + template + void GradientCrystalPlasticitySolver:: solve_nonlinear_system(const bool flag_skip_extrapolation) { // Terminal and log output nonlinear_solver_logger.add_break( "Step " + std::to_string(discrete_time.get_step_number() + 1) + - ": Solving for t = " + std::to_string(discrete_time.get_next_time()) + - " with dt = " + std::to_string(discrete_time.get_next_step_size())); + ": Solving for t = " + + std::to_string(discrete_time.get_next_time()) + + " with dt = " + + std::to_string(discrete_time.get_next_step_size())); nonlinear_solver_logger.log_headers_to_terminal(); - // Initialize local variables - gCP::GradientCrystalPlasticitySolver::SolverData solver_data; - // Internal variables' values at the previous step are stored prepare_quadrature_point_history(); @@ -88,23 +103,24 @@ namespace gCP // Active set reset_internal_newton_method_constraints(); + // Solution algorithm switch (parameters.solution_algorithm) { case RunTimeParameters::SolutionAlgorithm::Monolithic: { - monolithic_algorithm(solver_data); + monolithic_algorithm(); } break; case RunTimeParameters::SolutionAlgorithm::Bouncing: { - bouncing_algorithm(solver_data); + bouncing_algorithm(); } break; case RunTimeParameters::SolutionAlgorithm::Embracing: { - embracing_algorihtm(solver_data); + embracing_algorihtm(); } break; @@ -116,23 +132,34 @@ namespace gCP break; } + // Compute and store the opening displacement at the quadrature + // points based on the converged solution store_effective_opening_displacement_in_quadrature_history(); + // Update the actual solution vector fe_field->solution = trial_solution; *pcout << std::endl; - - return (std::make_tuple( - solver_data.flag_successful_convergence, - solver_data.nonlinear_iteration)); } template - void GradientCrystalPlasticitySolver::monolithic_algorithm( - SolverData &solver_data) + void GradientCrystalPlasticitySolver::monolithic_algorithm() { + // Instance local variables and references + unsigned int nonlinear_iteration = 0; + + bool flag_successful_convergence = false, + flag_compute_active_set = true; + + dealii::Vector residual_l2_norms, old_residual_l2_norms; + + line_search = + std::make_unique( + parameters.monolithic_algorithm_parameters. + monolithic_system_solver_parameters.line_search_parameters); + const RunTimeParameters::NewtonRaphsonParameters &newton_parameters = parameters.monolithic_algorithm_parameters. monolithic_system_solver_parameters.newton_parameters; @@ -141,20 +168,14 @@ namespace gCP &krylov_parameters = parameters.monolithic_algorithm_parameters. monolithic_system_solver_parameters.krylov_parameters; - line_search = - std::make_unique( - parameters.monolithic_algorithm_parameters. - monolithic_system_solver_parameters.line_search_parameters); - // Newton-Raphson loop do { - // Increment iteration count - (solver_data.nonlinear_iteration)++; - //nonlinear_iteration++; + // Increase iteration counter + nonlinear_iteration++; AssertThrow( - solver_data.nonlinear_iteration <= + nonlinear_iteration <= newton_parameters.n_max_iterations, dealii::ExcMessage( "The nonlinear solver has reach the given maximum number of " @@ -162,88 +183,53 @@ namespace gCP std::to_string( newton_parameters.n_max_iterations) + ").")); - if (parameters.constitutive_laws_parameters. - scalar_microstress_law_parameters.flag_rate_independent) - { - if (solver_data.flag_compute_active_set) - { - determine_active_set(); - - determine_inactive_set(); - - reset_inactive_set_values(); - - solver_data.flag_compute_active_set = false; - } - } - else - { - locally_owned_active_set = - fe_field->get_locally_owned_plastic_slip_dofs(); - } + // Determine the active (and also inactive) set + active_set_algorithm(flag_compute_active_set); + // Output for debugging purposes debug_output(); // Constraints distribution (Done later to obtain a continous // start value for the active set determination) Temporary code - { - dealii::LinearAlgebraTrilinos::MPI::BlockVector - distributed_trial_solution = - fe_field->get_distributed_vector_instance(trial_solution); - - fe_field->get_affine_constraints().distribute( - distributed_trial_solution); - - trial_solution = distributed_trial_solution; - } + distribute_affine_constraints_to_trial_solution(); // Preparations for the Newton-Update and Line-Search - { - store_trial_solution(); + store_trial_solution(); - reset_and_update_quadrature_point_history(); - } + reset_and_update_quadrature_point_history(); // Assemble linear system - assemble_residual(); - assemble_jacobian(); - // Residuals - solver_data.residual_l2_norms = fe_field->get_l2_norms(residual); + assemble_residual(); - const double initial_objective_function_value = - solver_data.residual_l2_norms[0]; + // Store current l2-norm values and initial objective + // function + residual_l2_norms = fe_field->get_sub_l2_norms(residual); - const std::vector initial_objective_function_values = - LineSearch::get_objective_function_values( - std::vector( - solver_data.residual_l2_norms.begin() + 1, - solver_data.residual_l2_norms.end())); + old_residual_l2_norms = residual_l2_norms; - solver_data.old_residual_l2_norms = - solver_data.residual_l2_norms; + line_search->reinit( + LineSearch::get_objective_function_value( + residual_l2_norms.l2_norm()), + discrete_time.get_step_number() + 1, + nonlinear_iteration); // Terminal and log output - if (solver_data.nonlinear_iteration == 1) + if (nonlinear_iteration == 1) { update_and_output_nonlinear_solver_logger( - solver_data.residual_l2_norms); + residual_l2_norms); } - // Newton-Update + // Newton-Raphson update const unsigned int n_krylov_iterations = solve_linearized_system( krylov_parameters, 0, - solver_data.residual_l2_norms[0]); + residual_l2_norms.l2_norm()); - double relaxation_parameter = - newton_parameters.relaxation_parameter; - - std::vector relaxation_parameters( - 2, - newton_parameters.relaxation_parameter); + double relaxation_parameter = 1.0; update_trial_solution(relaxation_parameter); @@ -253,172 +239,291 @@ namespace gCP // (and possibly Line-Search) assemble_residual(); - solver_data.residual_l2_norms = fe_field->get_l2_norms(residual); - - double objective_function_value = - LineSearch::get_objective_function_value( - solver_data.residual_l2_norms[0]); - - std::vector objective_function_values = - LineSearch::get_objective_function_values( - std::vector( - solver_data.residual_l2_norms.begin() + 1, - solver_data.residual_l2_norms.end())); + residual_l2_norms = fe_field->get_sub_l2_norms(residual); // Line search algorithm if (newton_parameters.flag_line_search) { - line_search->reinit( - initial_objective_function_value, - discrete_time.get_step_number() + 1, - solver_data.nonlinear_iteration); + relaxation_parameter = line_search_algorithm(); - while (!line_search->suficient_descent_condition( - objective_function_value, relaxation_parameter)) - { - relaxation_parameter = - line_search->get_lambda(objective_function_value, - relaxation_parameter); + residual_l2_norms = fe_field->get_sub_l2_norms(residual); + } - reset_trial_solution(); + // Terminal and log output + { + const dealii::Vector newton_update_l2_norms = + fe_field->get_sub_l2_norms(newton_update); - update_trial_solution(relaxation_parameter); + const double order_of_convergence = + std::log(residual_l2_norms.l2_norm()) / + std::log(old_residual_l2_norms.l2_norm()); - reset_and_update_quadrature_point_history(); + update_and_output_nonlinear_solver_logger( + nonlinear_iteration, + n_krylov_iterations, + line_search->get_n_iterations(), + newton_update_l2_norms, + residual_l2_norms, + order_of_convergence, + relaxation_parameter); + } - objective_function_value = assemble_residual(); - } + //slip_rate_output(false); - /* - line_search->reinit( - initial_objective_function_values, - discrete_time.get_step_number() + 1, - nonlinear_iteration); + // Convergence check + flag_successful_convergence = + residual_l2_norms.l2_norm() < + newton_parameters.absolute_tolerance; + + // Check if the active set changed + if (flag_successful_convergence && + parameters.constitutive_laws_parameters. + scalar_microstress_law_parameters.flag_rate_independent) + { + const dealii::IndexSet original_active_set = + locally_owned_active_set; + + flag_compute_active_set = true; + + active_set_algorithm(flag_compute_active_set); - while (!line_search->suficient_descent_condition( - objective_function_values, - relaxation_parameters)) + if (original_active_set != locally_owned_active_set) { - std::cout << "Before: " - << relaxation_parameters[0] << ", " - << relaxation_parameters[1] << std::endl; + flag_successful_convergence = false; - relaxation_parameters = - line_search->get_lambdas( - objective_function_values, - relaxation_parameters); + nonlinear_iteration = 0; - std::cout << "After: " - << relaxation_parameters[0] << ", " - << relaxation_parameters[1] << std::endl; + reset_trial_solution(true); - reset_trial_solution(); + *pcout + << std::endl + << "Active set mismatch. Restarting Newton-loop with the " + << "new active set..." << std::endl << std::endl; + } + } + } while (!flag_successful_convergence); + } - update_trial_solution(std::min(relaxation_parameters[0], - relaxation_parameters[1])); - reset_and_update_quadrature_point_history(); - assemble_residual(); + template + void GradientCrystalPlasticitySolver::bouncing_algorithm() + { + // Declare and initilize local variables and references + unsigned int macro_nonlinear_iteration = 0, + micro_nonlinear_iteration = 0; - residual_l2_norms = - compute_residual_l2_norms(); + bool flag_successful_convergence = false, + flag_successful_macro_convergence = false, + flag_successful_micro_convergence = false, + flag_compute_active_set = true; - objective_function_values = - LineSearch::get_objective_function_values( - std::vector( - residual_l2_norms.begin() + 1, - residual_l2_norms.end())); - } */ - } + dealii::Vector residual_l2_norms, old_residual_l2_norms; - // Terminal and log output + line_search = + std::make_unique( + parameters.staggered_algorithm_parameters. + linear_momentum_solver_parameters.line_search_parameters); + + const RunTimeParameters::NewtonRaphsonParameters + ¯o_newton_parameters = parameters. + staggered_algorithm_parameters. + linear_momentum_solver_parameters.newton_parameters, + µ_newton_parameters = parameters. + staggered_algorithm_parameters.pseudo_balance_solver_parameters. + newton_parameters; + + const RunTimeParameters::KrylovParameters + ¯o_krylov_parameters = parameters. + staggered_algorithm_parameters. + linear_momentum_solver_parameters.krylov_parameters, + µ_krylov_parameters = parameters. + staggered_algorithm_parameters.pseudo_balance_solver_parameters. + krylov_parameters; + + // The distribution corresponds to another method for the sake of + // the monolithic algorithm + distribute_affine_constraints_to_trial_solution(); + + // Macro-Newton-Raphson loop + do + { + // Convergence check + residual_l2_norms = fe_field->get_sub_l2_norms(residual); + + flag_successful_macro_convergence = + residual_l2_norms[0] < + macro_newton_parameters.absolute_tolerance; + + if (flag_successful_macro_convergence) { - const std::vector newton_update_l2_norms = - fe_field->get_l2_norms(newton_update); + // If the linear momentum balance converges with the plastic + // slips values of the lattest converged micro Newton-Raphson + // loop, the solution is considered as found. + if (flag_successful_micro_convergence) + { + flag_successful_convergence = true; - solver_data.residual_l2_norms = - fe_field->get_l2_norms(residual); + continue; + } + // Otherwise compute a new trial solution for the plastic slips + else + { + // Determine the active (and also inactive) set + active_set_algorithm(flag_compute_active_set); - const double order_of_convergence = - std::log(solver_data.residual_l2_norms[0]) / - std::log(solver_data.old_residual_l2_norms[0]); + // Micro-Newton-Raphson loop + do + { + // Increase iteration counter + micro_nonlinear_iteration++; - update_and_output_nonlinear_solver_logger( - solver_data.nonlinear_iteration, - n_krylov_iterations, - line_search->get_n_iterations(), - newton_update_l2_norms, - solver_data.residual_l2_norms, - order_of_convergence, - relaxation_parameter); - } + // Preparations for the Newton-Update and Line-Search + store_trial_solution(); - solver_data.residual_l2_norms = fe_field->get_l2_norms(residual); + reset_and_update_quadrature_point_history(); - //slip_rate_output(false); + // Assemble linear system; + assemble_jacobian(); - solver_data.flag_successful_convergence = - solver_data.residual_l2_norms[0] < - newton_parameters.absolute_tolerance; + assemble_residual(); - if (solver_data.flag_successful_convergence && - parameters.constitutive_laws_parameters. - scalar_microstress_law_parameters.flag_rate_independent) - { - const dealii::IndexSet original_active_set = - locally_owned_active_set; + // Store current l2-norm values and initial objective + // function + residual_l2_norms = fe_field->get_sub_l2_norms(residual); - // Active set - reset_internal_newton_method_constraints(); + old_residual_l2_norms = residual_l2_norms; - if (parameters.constitutive_laws_parameters. - scalar_microstress_law_parameters.flag_rate_independent) - { - determine_active_set(); + line_search->reinit( + LineSearch::get_objective_function_value( + residual_l2_norms[1]), + discrete_time.get_step_number() + 1, + micro_nonlinear_iteration); + + // Terminal and log output + { + + } + + // Newton-Raphson update + const unsigned int n_krylov_iterations = + solve_linearized_system( + micro_krylov_parameters, + 1, + residual_l2_norms[1]); + + double relaxation_parameter = 1.0; + + update_trial_solution(relaxation_parameter, 1); + + reset_and_update_quadrature_point_history(); + + // Evaluate new trial solution + assemble_residual(); + + residual_l2_norms = + fe_field->get_sub_l2_norms(residual); - determine_inactive_set(); + // Line search + if (micro_newton_parameters.flag_line_search) + { - reset_inactive_set_values(); + } + + // Terminal and log output + { + (void)n_krylov_iterations; + } + + // Convergence check + flag_successful_micro_convergence = + residual_l2_norms[1] < + micro_newton_parameters.absolute_tolerance; + + } while (!flag_successful_micro_convergence); + } } else { - locally_owned_active_set = - fe_field->get_locally_owned_plastic_slip_dofs(); - } + // Increase iteration counter + macro_nonlinear_iteration++; - if (original_active_set != locally_owned_active_set) - { - solver_data.flag_successful_convergence = false; + // Preparations for the Newton-Update and Line-Search + store_trial_solution(); - solver_data.nonlinear_iteration = 0; + reset_and_update_quadrature_point_history(); - reset_trial_solution(true); + // Assemble linear system + if (macro_nonlinear_iteration == 1) + { + // The submatrix at (0,0) is a constant, therefore it only + // needs to be assembled once + assemble_jacobian(); + } - *pcout - << std::endl - << "Active set mismatch. Restarting Newton-loop with the " - << "new active set..." << std::endl << std::endl; - } - } + assemble_residual(); - } while (!solver_data.flag_successful_convergence); - } + // Store current l2-norm values and initial objective function + residual_l2_norms = fe_field->get_sub_l2_norms(residual); + old_residual_l2_norms = residual_l2_norms; + line_search->reinit( + LineSearch::get_objective_function_value( + residual_l2_norms[0]), + discrete_time.get_step_number() + 1, + macro_nonlinear_iteration); - template - void GradientCrystalPlasticitySolver::bouncing_algorithm( - SolverData &solver_data) - { + // Terminal and log output () + { + + } + + // Newton-Raphson update + const unsigned int n_krylov_iterations = + solve_linearized_system( + macro_krylov_parameters, + 0, + residual_l2_norms[0]); + + double relaxation_parameter = 1.0; + + update_trial_solution(relaxation_parameter, 0); + + reset_and_update_quadrature_point_history(); + + // Evaluate new trial solution + assemble_residual(); + + residual_l2_norms = + fe_field->get_sub_l2_norms(residual); + + // Line search algorithm + if (macro_newton_parameters.flag_line_search) + { + + } + + // Terminal and log output + { + (void)n_krylov_iterations; + } + // Convergence check + + flag_successful_macro_convergence = + residual_l2_norms[0] < + macro_newton_parameters.absolute_tolerance; + + flag_successful_micro_convergence = false; + } + } while (!flag_successful_convergence); } template - void GradientCrystalPlasticitySolver::embracing_algorihtm( - SolverData &solver_data) + void GradientCrystalPlasticitySolver::embracing_algorihtm() { } @@ -564,7 +669,8 @@ namespace gCP template void GradientCrystalPlasticitySolver::update_trial_solution( - const double relaxation_parameter) + const double relaxation_parameter, + const unsigned int block_id) { dealii::LinearAlgebraTrilinos::MPI::BlockVector distributed_trial_solution = @@ -574,13 +680,14 @@ namespace gCP distributed_newton_update = fe_field->get_distributed_vector_instance(newton_update); - distributed_trial_solution.add( - relaxation_parameter, distributed_newton_update); + distributed_trial_solution.block(block_id).add( + relaxation_parameter, distributed_newton_update.block(block_id)); fe_field->get_affine_constraints().distribute( distributed_trial_solution); - trial_solution = distributed_trial_solution; + trial_solution.block(block_id) = + distributed_trial_solution.block(block_id); } @@ -750,7 +857,7 @@ namespace gCP template void GradientCrystalPlasticitySolver:: update_and_output_nonlinear_solver_logger( - const std::vector residual_l2_norms) + const dealii::Vector residual_l2_norms) { nonlinear_solver_logger.update_value("N-Itr", 0); @@ -765,11 +872,11 @@ namespace gCP nonlinear_solver_logger.update_value("(NS_G)_L2", 0); nonlinear_solver_logger.update_value("(R)_L2", - residual_l2_norms[0]); + residual_l2_norms.l2_norm()); nonlinear_solver_logger.update_value("(R_U)_L2", - residual_l2_norms[1]); + residual_l2_norms[0]); nonlinear_solver_logger.update_value("(R_G)_L2", - residual_l2_norms[2]); + residual_l2_norms[1]); nonlinear_solver_logger.update_value("C-Rate", 0); @@ -783,13 +890,13 @@ namespace gCP template void GradientCrystalPlasticitySolver:: update_and_output_nonlinear_solver_logger( - const unsigned int nonlinear_iteration, - const unsigned int n_krylov_iterations, - const unsigned int n_line_search_iterations, - const std::vector newton_update_l2_norms, - const std::vector residual_l2_norms, - const double order_of_convergence, - const double relaxation_parameter) + const unsigned int nonlinear_iteration, + const unsigned int n_krylov_iterations, + const unsigned int n_line_search_iterations, + const dealii::Vector newton_update_l2_norms, + const dealii::Vector residual_l2_norms, + const double order_of_convergence, + const double relaxation_parameter) { nonlinear_solver_logger.update_value("N-Itr", nonlinear_iteration); @@ -799,19 +906,19 @@ namespace gCP n_line_search_iterations); nonlinear_solver_logger.update_value("(NS)_L2", relaxation_parameter * - newton_update_l2_norms[0]); + newton_update_l2_norms.l2_norm()); nonlinear_solver_logger.update_value("(NS_U)_L2", relaxation_parameter * - newton_update_l2_norms[1]); + newton_update_l2_norms[0]); nonlinear_solver_logger.update_value("(NS_G)_L2", relaxation_parameter * - newton_update_l2_norms[2]); + newton_update_l2_norms[1]); nonlinear_solver_logger.update_value("(R)_L2", - residual_l2_norms[0]); + residual_l2_norms.l2_norm()); nonlinear_solver_logger.update_value("(R_U)_L2", - residual_l2_norms[1]); + residual_l2_norms[0]); nonlinear_solver_logger.update_value("(R_G)_L2", - residual_l2_norms[2]); + residual_l2_norms[1]); nonlinear_solver_logger.update_value("C-Rate", order_of_convergence); @@ -826,11 +933,15 @@ namespace gCP -template void gCP::GradientCrystalPlasticitySolver<2>::extrapolate_initial_trial_solution(const bool); -template void gCP::GradientCrystalPlasticitySolver<3>::extrapolate_initial_trial_solution(const bool); +template void gCP::GradientCrystalPlasticitySolver<2>:: +extrapolate_initial_trial_solution(const bool); +template void gCP::GradientCrystalPlasticitySolver<3>:: +extrapolate_initial_trial_solution(const bool); -template std::tuple gCP::GradientCrystalPlasticitySolver<2>::solve_nonlinear_system(const bool); -template std::tuple gCP::GradientCrystalPlasticitySolver<3>::solve_nonlinear_system(const bool); +template void gCP::GradientCrystalPlasticitySolver<2>:: +solve_nonlinear_system(const bool); +template void gCP::GradientCrystalPlasticitySolver<3>:: +solve_nonlinear_system(const bool); template unsigned int gCP::GradientCrystalPlasticitySolver<2>:: solve_linearized_system( @@ -843,8 +954,14 @@ solve_linearized_system( const unsigned int, const double); -template void gCP::GradientCrystalPlasticitySolver<2>::update_trial_solution(const double); -template void gCP::GradientCrystalPlasticitySolver<3>::update_trial_solution(const double); +template void gCP::GradientCrystalPlasticitySolver<2>:: +update_trial_solution( + const double, + const unsigned int); +template void gCP::GradientCrystalPlasticitySolver<3>:: +update_trial_solution( + const double, + const unsigned int); template void gCP::GradientCrystalPlasticitySolver<2>:: update_trial_solution(const std::vector);