From 9788e909767dd2431a8964e50182d66054dd259c Mon Sep 17 00:00:00 2001 From: Manav Bhatia Date: Thu, 26 Mar 2020 16:44:41 -0500 Subject: [PATCH] -- Added example for transient conduction analysis -- Updated tutorial page to list the new example -- Bug-fix for transient direct sensitivity analysis -- Modification to fluid example sensitivity analysis consistent with the bug-fix --- doc/tutorials.dox | 1 + examples/CMakeLists.txt | 1 + examples/conduction/CMakeLists.txt | 2 + examples/conduction/example_1/CMakeLists.txt | 12 + examples/conduction/example_1/example_1.cpp | 343 ++++++++++++++++++ examples/fluid/example_1/example_1.cpp | 11 +- .../heat_conduction_elem_base.cpp | 84 ++++- .../heat_conduction_transient_assembly.cpp | 21 +- .../first_order_newmark_transient_solver.cpp | 5 +- .../second_order_newmark_transient_solver.cpp | 3 + src/solver/transient_solver_base.cpp | 16 +- src/solver/transient_solver_base.h | 14 +- 12 files changed, 488 insertions(+), 25 deletions(-) create mode 100644 examples/conduction/CMakeLists.txt create mode 100644 examples/conduction/example_1/CMakeLists.txt create mode 100644 examples/conduction/example_1/example_1.cpp diff --git a/doc/tutorials.dox b/doc/tutorials.dox index 5a4756ad..13f2f020 100644 --- a/doc/tutorials.dox +++ b/doc/tutorials.dox @@ -16,5 +16,6 @@ - \subpage structural_example_7 - \subpage fluid_example_1 - \subpage fsi_example_1 + - \subpage conduction_example_1 */ diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 1f04258c..60a27fb4 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,3 +1,4 @@ +add_subdirectory(conduction) add_subdirectory(structural) add_subdirectory(fluid) add_subdirectory(fsi) diff --git a/examples/conduction/CMakeLists.txt b/examples/conduction/CMakeLists.txt new file mode 100644 index 00000000..09263b08 --- /dev/null +++ b/examples/conduction/CMakeLists.txt @@ -0,0 +1,2 @@ +add_subdirectory(example_1) # 1D transient conduction + diff --git a/examples/conduction/example_1/CMakeLists.txt b/examples/conduction/example_1/CMakeLists.txt new file mode 100644 index 00000000..f6068106 --- /dev/null +++ b/examples/conduction/example_1/CMakeLists.txt @@ -0,0 +1,12 @@ +add_executable(conduction_example_1 + example_1.cpp) + +target_include_directories(conduction_example_1 PRIVATE + ${CMAKE_CURRENT_LIST_DIR} + ${PROJECT_SOURCE_DIR} + ${PROJECT_SOURCE_DIR}/src) +target_link_libraries(conduction_example_1 mast) + +install(TARGETS conduction_example_1 + RUNTIME DESTINATION ${CMAKE_INSTALL_PREFIX}/examples) + diff --git a/examples/conduction/example_1/example_1.cpp b/examples/conduction/example_1/example_1.cpp new file mode 100644 index 00000000..1c878e38 --- /dev/null +++ b/examples/conduction/example_1/example_1.cpp @@ -0,0 +1,343 @@ +/* + * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit + * Copyright (C) 2013-2020 Manav Bhatia and MAST authors + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ + +// C/C++ includes. +#include + +// MAST includes. +#include "examples/base/input_wrapper.h" +#include "base/nonlinear_system.h" +#include "base/transient_assembly.h" +#include "base/boundary_condition_base.h" +#include "base/field_function_base.h" +#include "base/parameter.h" +#include "base/constant_field_function.h" +#include "base/physics_discipline_base.h" +#include "property_cards/isotropic_material_property_card.h" +#include "property_cards/solid_1d_section_element_property_card.h" +#include "boundary_condition/dirichlet_boundary_condition.h" +#include "heat_conduction/heat_conduction_system_initialization.h" +#include "heat_conduction/heat_conduction_nonlinear_assembly.h" +#include "heat_conduction/heat_conduction_transient_assembly.h" +#include "solver/first_order_newmark_transient_solver.h" + +// libMesh includes +#include "libmesh/libmesh.h" +#include "libmesh/string_to_enum.h" +#include "libmesh/enum_order.h" +#include "libmesh/enum_fe_family.h" +#include "libmesh/enum_elem_type.h" // ElemType +#include "libmesh/fe_type.h" // FEFamily, Order +#include "libmesh/serial_mesh.h" +#include "libmesh/equation_systems.h" +#include "libmesh/dof_map.h" +#include "libmesh/numeric_vector.h" +#include "libmesh/exodusII_io.h" +#include "libmesh/nonlinear_solver.h" +#include "libmesh/mesh_generation.h" + +// Forward declerations +void compute_transient_solution(MAST::PhysicsDisciplineBase& discipline, + MAST::SystemInitialization& sys_init, + MAST::Examples::GetPotWrapper& input); + +void compute_transient_sensitivity(MAST::PhysicsDisciplineBase& discipline, + MAST::SystemInitialization& sys_init, + MAST::Examples::GetPotWrapper& input, + MAST::Parameter& p); + + +int main(int argc, const char** argv) +{ + // BEGIN_TRANSLATE Extension of bar + // + // \tableofcontents + // + // This example solves an axial bar extension problem. + // + // Initialize libMesh library. + libMesh::LibMeshInit init(argc, argv); + + MAST::Examples::GetPotWrapper + input(argc, argv, "input"); + + // Create Mesh object on default MPI communicator and generate a line mesh (5 elements, 10 units long). + // Note that in libMesh, all meshes are parallel by default in the sense that the equations on the mesh are solved in parallel by PETSc. + // A "ReplicatedMesh" is one where all MPI processes have the full mesh in memory, as compared to a "DistributedMesh" where the mesh is + // "chunked" up and distributed across processes, each having their own piece. + libMesh::ReplicatedMesh mesh(init.comm()); + libMesh::MeshTools::Generation::build_line(mesh, 20, 0.0, 10.0); + mesh.print_info(); + mesh.boundary_info->print_info(); + + // Create EquationSystems object, which is a container for multiple systems of equations that are defined on a given mesh. + libMesh::EquationSystems equation_systems(mesh); + + // Add system of type MAST::NonlinearSystem (which wraps libMesh::NonlinearImplicitSystem) to the EquationSystems container. + MAST::NonlinearSystem & system = equation_systems.add_system("conduction"); + + // Create a finite element type for the system. Here we use first order + // Lagrangian-type finite elements. + libMesh::FEType fetype(libMesh::FIRST, libMesh::LAGRANGE); + + // Initialize the system to the correct set of variables for a conduction + // analysis. In libMesh this is analogous to adding variables (each with + // specific finite element type/order to the system for a particular + // system of equations. + MAST::HeatConductionSystemInitialization conduction_system(system, + system.name(), + fetype); + + // Initialize a new conduction discipline using equation_systems. + MAST::PhysicsDisciplineBase discipline(equation_systems); + + // Initialize the equation system since we now know the size of our + // system matrices (based on mesh, element type, variables in the + // conduction_system) as well as the setup of dirichlet boundary conditions. + // This initialization process is basically a pre-processing step to + // preallocate storage and spread it across processors. + equation_systems.init(); + equation_systems.print_info(); + + // Create parameters. + MAST::Parameter zero("zero", 0.0); + MAST::Parameter kappa_yy("kappa_yy", 5./6.); + MAST::Parameter kappa_zz("kappa_zz", 5./6.); + MAST::Parameter thickness_y("thy", input("thy", "section y-thickness", 0.06)); + MAST::Parameter thickness_z("thz", input("thz", "section z-thickness", 0.02)); + MAST::Parameter k ("k", input("k", "thermal conductivity", 190.)); + MAST::Parameter rho ("rho", input("rho", "material density", 2700.)); + MAST::Parameter cp ("cp", input("cp", "thermal capacitance", 864.)); + MAST::Parameter q ("flux",input("flux", "heat flux", -2.e3)); + + // Create ConstantFieldFunctions used to spread parameters throughout the model. + MAST::ConstantFieldFunction thy_f( "hy", thickness_y); + MAST::ConstantFieldFunction thz_f( "hz", thickness_z); + MAST::ConstantFieldFunction hyoff_f( "hy_off", zero); + MAST::ConstantFieldFunction hzoff_f( "hz_off", zero); + MAST::ConstantFieldFunction k_f( "k_th", k); + MAST::ConstantFieldFunction rho_f( "rho", rho); + MAST::ConstantFieldFunction cp_f( "cp", cp); + MAST::ConstantFieldFunction q_f( "heat_flux", q); + MAST::ConstantFieldFunction kappa_yy_f("Kappayy", kappa_yy); + MAST::ConstantFieldFunction kappa_zz_f("Kappazz", kappa_zz); + + // Initialize load. + // TODO - Switch this to a concentrated/point load on the right end of the bar. + MAST::BoundaryConditionBase right_end_flux(MAST::HEAT_FLUX); + right_end_flux.add(q_f); + discipline.add_side_load(1, right_end_flux); + + // Create the material property card ("card" is NASTRAN lingo) and the + // relevant parameters to it. An isotropic material needs thermal + // conductivity (k), density (rho) and thermal capacitance (cp) to + // describe its behavior. + MAST::IsotropicMaterialPropertyCard material; + material.add(k_f); + material.add(rho_f); + material.add(cp_f); + + // Create the section property card. Attach all property values. + MAST::Solid1DSectionElementPropertyCard section; + section.add(thy_f); + section.add(thz_f); + section.add(hyoff_f); + section.add(hzoff_f); + section.add(kappa_yy_f); + section.add(kappa_zz_f); + + // Specify a section orientation point and add it to the section. + RealVectorX orientation = RealVectorX::Zero(3); + orientation(1) = 1.0; + section.y_vector() = orientation; + + // Attach material to the card. + section.set_material(material); + section.set_diagonal_mass_matrix(true); + + // Initialize the section and specify the subdomain in the mesh that it applies to. + section.init(); + discipline.set_property_for_subdomain(0, section); + + // transient solution + compute_transient_solution(discipline, conduction_system, input); + + // transient solution sensitivity with respect to section y-thickness + compute_transient_sensitivity(discipline, conduction_system, input, thickness_y); + + // END_TRANSLATE + return 0; +} + + +// \section conduction_computation Computation +// \subsection conduction_transient_analysis Transient analysis +void compute_transient_solution(MAST::PhysicsDisciplineBase& discipline, + MAST::SystemInitialization& sys_init, + MAST::Examples::GetPotWrapper& input) { + + std::string + output_name = input("output_file_root", "prefix of output file names", "output"), + transient_output_name = output_name + "_transient.exo"; + + MAST::NonlinearSystem& sys = sys_init.system(); + libMesh::EquationSystems& eq_sys = sys.get_equation_systems(); + + // create the nonlinear assembly object + MAST::TransientAssembly assembly; + MAST::HeatConductionTransientAssemblyElemOperations elem_ops; + MAST::FirstOrderNewmarkTransientSolver solver; + + assembly.set_discipline_and_system(discipline, sys_init); + elem_ops.set_discipline_and_system(discipline, sys_init); + solver.set_discipline_and_system(discipline, sys_init); + solver.set_elem_operation_object(elem_ops); + + // initialize the solution to zero, or to something that the + // user may have provided + sys.solution->zero(); + sys.solution->close(); + sys.update(); + + // file to write the solution for visualization + libMesh::ExodusII_IO transient_output(sys.get_mesh()); + + unsigned int + t_step = 0, + n_steps = input("n_transient_steps", "number of transient time-steps", 100); + solver.dt = input("dt", "time-step size", 1.e+3); + solver.beta = input("beta", "Newmark solver beta parameter ", 0.5); + + // ask the solver to update the initial condition for d2(X)/dt2 + // This is recommended only for the initial time step, since the time + // integration scheme updates the velocity and acceleration at + // each subsequent iterate + solver.solve_highest_derivative_and_advance_time_step(assembly); + + // loop over time steps + while (t_step < n_steps) { + + libMesh::out + << "Time step: " << t_step + << " : t = " << sys.time + << " : dt = " << solver.dt + << " : xdot-L2 = " << solver.velocity().l2_norm() + << std::endl; + + // write the time-step + transient_output.write_timestep(transient_output_name, + eq_sys, + t_step+1, + sys.time); + std::ostringstream oss_sol; + oss_sol << output_name << "_sol_t_" << t_step; + sys_init.system().write_out_vector(*sys.solution, "data", oss_sol.str(), true); + + // solve for the time-step + solver.solve(assembly); + solver.advance_time_step(); + + // update time step + t_step++; + } +} + + +// \subsection conduction_transient_sensitivity_analysis Transient sensitivity analysis +void compute_transient_sensitivity(MAST::PhysicsDisciplineBase& discipline, + MAST::SystemInitialization& sys_init, + MAST::Examples::GetPotWrapper& input, + MAST::Parameter& p) { + + std::string + output_name = input("output_file_root", "prefix of output file names", "output"), + transient_output_name = output_name + "_transient_sensitivity_" + p.name() + ".exo", + nonlinear_sol_dir = input("nonlinear_sol_dir", "directory containing the location of nonlinear solutions", "data"); + + MAST::NonlinearSystem& sys = sys_init.system(); + libMesh::EquationSystems& eq_sys = sys.get_equation_systems(); + + // create the nonlinear assembly object + MAST::TransientAssembly assembly; + MAST::HeatConductionTransientAssemblyElemOperations elem_ops; + MAST::FirstOrderNewmarkTransientSolver solver; + + assembly.set_discipline_and_system(discipline, sys_init); + elem_ops.set_discipline_and_system(discipline, sys_init); + solver.set_discipline_and_system(discipline, sys_init); + solver.set_elem_operation_object(elem_ops); + + // initial condition for solution is zero temperature and initial condition + // for sensitivity of temperature is also assumed to be zero. + sys.solution->zero(); + sys.solution->close(); + sys.update(); + + // file to write the solution for visualization + libMesh::ExodusII_IO transient_output(sys.get_mesh()); + + unsigned int + t_step = 0, + n_steps = input("n_transient_steps", "number of transient time-steps", 100); + solver.dt = input("dt", "time-step size", 1.e+3); + sys.time = 0.; + solver.beta = input("beta", "Newmark solver beta parameter ", 0.5); + + // ask the solver to update the initial condition for d2(X)/dt2 + // This is recommended only for the initial time step, since the time + // integration scheme updates the velocity and acceleration at + // each subsequent iterate + solver.solve_highest_derivative_and_advance_time_step(assembly, false); + solver.solve_highest_derivative_and_advance_time_step_with_sensitivity(assembly, p); + + // loop over time steps + while (t_step < n_steps-1) { + + libMesh::out + << "Time step: " << t_step + << " : t = " << sys.time + << " : xdot-L2 = " << solver.velocity_sensitivity().l2_norm() + << std::endl; + + // write the time-step + sys.solution->swap(solver.solution_sensitivity()); + transient_output.write_timestep(transient_output_name, + eq_sys, + t_step+1, + sys.time); + sys.solution->swap(solver.solution_sensitivity()); + + std::ostringstream oss_sol; + oss_sol << output_name << "_sol_t_" << t_step+1; + sys.read_in_vector(*sys.solution, nonlinear_sol_dir, oss_sol.str(), true); + solver.update_velocity(solver.velocity(), *sys.solution); + sys.update(); + + + // solve for the sensitivity time-step + solver.sensitivity_solve(assembly, p); + solver.advance_time_step(false); + solver.advance_time_step_with_sensitivity(); + + // update time step counter + t_step++; + } + +} + diff --git a/examples/fluid/example_1/example_1.cpp b/examples/fluid/example_1/example_1.cpp index adc0f870..fb4bfc38 100644 --- a/examples/fluid/example_1/example_1.cpp +++ b/examples/fluid/example_1/example_1.cpp @@ -1208,10 +1208,11 @@ class FlowAnalysis { // This is recommended only for the initial time step, since the time // integration scheme updates the velocity and acceleration at // each subsequent iterate + solver.solve_highest_derivative_and_advance_time_step(assembly, false); solver.solve_highest_derivative_and_advance_time_step_with_sensitivity(assembly, p); // loop over time steps - while (t_step < n_steps) { + while (t_step < n_steps-1) { libMesh::out << "Time step: " << t_step @@ -1230,11 +1231,12 @@ class FlowAnalysis { _sys->solution->swap(solver.solution_sensitivity()); } - std::ostringstream oss; - oss << nonlinear_sol_root << t_step; + std::ostringstream oss, oss_sens; + oss << nonlinear_sol_root << t_step+1; _sys->read_in_vector(*_sys->solution, nonlinear_sol_dir, oss.str(), true); + solver.update_velocity(solver.velocity(), *_sys->solution); _sys->update(); - oss << "_sens_t"; + oss_sens << nonlinear_sol_root << t_step << "_sens_t"; _sys->write_out_vector(/*solver.dt, _sys->time,*/ solver.solution_sensitivity(), "data", oss.str(), true); // solve for the sensitivity time-step @@ -1247,6 +1249,7 @@ class FlowAnalysis { << std::setw(30) << force.output_sensitivity_total(p) << std::endl; solver.sensitivity_solve(assembly, p); + solver.advance_time_step(false); solver.advance_time_step_with_sensitivity(); // update time step counter diff --git a/src/heat_conduction/heat_conduction_elem_base.cpp b/src/heat_conduction/heat_conduction_elem_base.cpp index c95bb898..11368158 100644 --- a/src/heat_conduction/heat_conduction_elem_base.cpp +++ b/src/heat_conduction/heat_conduction_elem_base.cpp @@ -173,7 +173,8 @@ MAST::HeatConductionElementBase::velocity_residual (bool request_jacobian, mat_n2n2 = RealMatrixX::Zero(n_phi, n_phi); RealVectorX vec1 = RealVectorX::Zero(1), - vec2_n2 = RealVectorX::Zero(n_phi); + vec2_n2 = RealVectorX::Zero(n_phi), + local_f = RealVectorX::Zero(n_phi); std::unique_ptr > capacitance = _property.thermal_capacitance_matrix(*this); @@ -192,9 +193,9 @@ MAST::HeatConductionElementBase::velocity_residual (bool request_jacobian, Bmat.right_multiply(vec1, _vel); // B * T_dot Bmat.vector_mult_transpose(vec2_n2, vec1); // B^T * B * T_dot - f += JxW[qp] * material_mat(0,0) * vec2_n2; // (rho*cp)*JxW B^T B T_dot + local_f += JxW[qp] * material_mat(0,0) * vec2_n2; // (rho*cp)*JxW B^T B T_dot - if (request_jacobian) { + if (request_jacobian || _property.if_diagonal_mass_matrix()) { Bmat.right_multiply_transpose(mat_n2n2, Bmat); // B^T B jac_xdot += JxW[qp] * material_mat(0,0) * mat_n2n2; // B^T B * JxW (rho*cp) @@ -215,7 +216,27 @@ MAST::HeatConductionElementBase::velocity_residual (bool request_jacobian, } } - + // diagonalize the matrix and compute the residual based on that. + if (_property.if_diagonal_mass_matrix()) { + + for (unsigned int i=0; i (_active_sol_function)->clear_element_quadrature_point_solution(); @@ -701,10 +722,13 @@ velocity_residual_sensitivity (const MAST::FunctionBase& p, RealMatrixX material_mat = RealMatrixX::Zero(dim, dim), - mat_n2n2 = RealMatrixX::Zero(n_phi, n_phi); + mat_n2n2 = RealMatrixX::Zero(n_phi, n_phi), + local_jac_xdot = RealMatrixX::Zero(n_phi, n_phi); + RealVectorX vec1 = RealVectorX::Zero(1), - vec2_n2 = RealVectorX::Zero(n_phi); + vec2_n2 = RealVectorX::Zero(n_phi), + local_f = RealVectorX::Zero(n_phi); std::unique_ptr > capacitance = _property.thermal_capacitance_matrix(*this); @@ -723,9 +747,28 @@ velocity_residual_sensitivity (const MAST::FunctionBase& p, Bmat.right_multiply(vec1, _vel); // B * T_dot Bmat.vector_mult_transpose(vec2_n2, vec1); // B^T * B * T_dot - f += JxW[qp] * material_mat(0,0) * vec2_n2; // (rho*cp)*JxW B^T B T_dot + local_f += JxW[qp] * material_mat(0,0) * vec2_n2; // (rho*cp)*JxW B^T B T_dot + + if (_property.if_diagonal_mass_matrix()) { + + Bmat.right_multiply_transpose(mat_n2n2, Bmat); // B^T B + local_jac_xdot += JxW[qp] * material_mat(0,0) * mat_n2n2; // B^T B * JxW (rho*cp) + } } + if (_property.if_diagonal_mass_matrix()) { + + for (unsigned int i=0; i @@ -757,10 +800,12 @@ velocity_residual_boundary_velocity (const MAST::FunctionBase& p, RealMatrixX material_mat = RealMatrixX::Zero(dim, dim), - mat_n2n2 = RealMatrixX::Zero(n_phi, n_phi); + mat_n2n2 = RealMatrixX::Zero(n_phi, n_phi), + local_jac_xdot = RealMatrixX::Zero(n_phi, n_phi); RealVectorX vec1 = RealVectorX::Zero(1), vec2_n2 = RealVectorX::Zero(n_phi), + local_f = RealVectorX::Zero(n_phi), vel = RealVectorX::Zero(dim); std::unique_ptr > capacitance = @@ -793,9 +838,30 @@ velocity_residual_boundary_velocity (const MAST::FunctionBase& p, Bmat.right_multiply(vec1, _vel); // B * T_dot Bmat.vector_mult_transpose(vec2_n2, vec1); // B^T * B * T_dot - f += JxW_Vn[qp] * material_mat(0,0) * vec2_n2; // (rho*cp)*JxW B^T B T_dot + local_f += JxW_Vn[qp] * material_mat(0,0) * vec2_n2; // (rho*cp)*JxW B^T B T_dot + + if (_property.if_diagonal_mass_matrix()) { + + Bmat.right_multiply_transpose(mat_n2n2, Bmat); // B^T B + local_jac_xdot += JxW_Vn[qp] * material_mat(0,0) * mat_n2n2; // B^T B * JxW (rho*cp) + } } + // diagonalize the matrix and compute the residual based on that. + if (_property.if_diagonal_mass_matrix()) { + + for (unsigned int i=0; i diff --git a/src/heat_conduction/heat_conduction_transient_assembly.cpp b/src/heat_conduction/heat_conduction_transient_assembly.cpp index e3e51d32..be5ba52d 100644 --- a/src/heat_conduction/heat_conduction_transient_assembly.cpp +++ b/src/heat_conduction/heat_conduction_transient_assembly.cpp @@ -88,8 +88,27 @@ elem_sensitivity_calculations(const MAST::FunctionBase& f, RealVectorX& f_m, RealVectorX& f_x) { - libmesh_error(); // to be implemented + libmesh_assert(_physics_elem); + + MAST::HeatConductionElementBase& e = + dynamic_cast(*_physics_elem); + + unsigned int + n = (unsigned int)f_m.size(); + RealMatrixX + dummy = RealMatrixX::Zero(n, n); + + f_m.setZero(); + f_x.setZero(); + + // assembly of the flux terms + e.internal_residual_sensitivity(f, f_x); + e.side_external_residual_sensitivity(f, f_x, _discipline->side_loads()); + e.volume_external_residual_sensitivity(f, f_x, _discipline->volume_loads()); + + //assembly of the capacitance term + e.velocity_residual_sensitivity(f, f_m); } diff --git a/src/solver/first_order_newmark_transient_solver.cpp b/src/solver/first_order_newmark_transient_solver.cpp index 461c256b..b4a50c6b 100644 --- a/src/solver/first_order_newmark_transient_solver.cpp +++ b/src/solver/first_order_newmark_transient_solver.cpp @@ -323,7 +323,10 @@ MAST::FirstOrderNewmarkTransientSolver:: elem_sensitivity_contribution_previous_timestep(const std::vector& prev_sols, RealVectorX& vec) { - // make sure that the assembly object is provided + // nothing to be done for highest derivative term + if (_if_highest_derivative_solution) return; + + // make sure that the assembly object is provided libmesh_assert(_assembly_ops); libmesh_assert_equal_to(prev_sols.size(), 2); diff --git a/src/solver/second_order_newmark_transient_solver.cpp b/src/solver/second_order_newmark_transient_solver.cpp index d260048a..ae381902 100644 --- a/src/solver/second_order_newmark_transient_solver.cpp +++ b/src/solver/second_order_newmark_transient_solver.cpp @@ -418,6 +418,9 @@ MAST::SecondOrderNewmarkTransientSolver:: elem_sensitivity_contribution_previous_timestep(const std::vector& prev_sols, RealVectorX& vec) { + // nothing to be done for highest derivative term + if (_if_highest_derivative_solution) return; + // make sure that the assembly object is provided libmesh_assert(_assembly_ops); libmesh_assert_equal_to(prev_sols.size(), 3); diff --git a/src/solver/transient_solver_base.cpp b/src/solver/transient_solver_base.cpp index 45164e08..7bd0c1a1 100644 --- a/src/solver/transient_solver_base.cpp +++ b/src/solver/transient_solver_base.cpp @@ -334,13 +334,16 @@ MAST::TransientSolverBase::sensitivity_solve(MAST::AssemblyBase& assembly, // ask the Newton solver to solve for the system solution _system->system().sensitivity_solve(*this, assembly, f); + + this->solution_sensitivity() = assembly.system().get_sensitivity_solution(); } void MAST::TransientSolverBase:: -solve_highest_derivative_and_advance_time_step(MAST::AssemblyBase& assembly) { +solve_highest_derivative_and_advance_time_step(MAST::AssemblyBase& assembly, + bool increment_time) { libmesh_assert(_first_step); libmesh_assert(_system); @@ -415,7 +418,7 @@ solve_highest_derivative_and_advance_time_step(MAST::AssemblyBase& assembly) { } // finally, update the system time - sys.time += dt; + if (increment_time) sys.time += dt; _first_step = false; } @@ -424,8 +427,9 @@ solve_highest_derivative_and_advance_time_step(MAST::AssemblyBase& assembly) { void MAST::TransientSolverBase:: -solve_highest_derivative_and_advance_time_step_with_sensitivity(MAST::AssemblyBase& assembly, - const MAST::FunctionBase& f) { +solve_highest_derivative_and_advance_time_step_with_sensitivity + (MAST::AssemblyBase& assembly, + const MAST::FunctionBase& f) { libmesh_assert(_first_sensitivity_step); libmesh_assert(_system); @@ -731,7 +735,7 @@ build_perturbed_local_quantities(const libMesh::NumericVector& current_dso void -MAST::TransientSolverBase::advance_time_step() { +MAST::TransientSolverBase::advance_time_step(bool increment_time) { libmesh_assert(_system); libmesh_assert(_discipline); @@ -757,7 +761,7 @@ MAST::TransientSolverBase::advance_time_step() { } // finally, update the system time - sys.time += dt; + if (increment_time) sys.time += dt; _first_step = false; } diff --git a/src/solver/transient_solver_base.h b/src/solver/transient_solver_base.h index 8ee002e5..b15e691e 100644 --- a/src/solver/transient_solver_base.h +++ b/src/solver/transient_solver_base.h @@ -191,9 +191,13 @@ namespace MAST { * Initializes the highest derivative solution using the solution * and its low order time derivatives at the specified time step. * Then advances the time step so that the solver is ready for - * time integration. + * time integration. If \p increment_time is \p true then + * the value fo time will be incremented in System. This may be useful when calling this + * method inside a sensitivity analysis routine. */ - void solve_highest_derivative_and_advance_time_step(MAST::AssemblyBase& assembly); + void + solve_highest_derivative_and_advance_time_step(MAST::AssemblyBase& assembly, + bool increment_time = true); /*! * solves for the sensitivity of highest derivative and advances @@ -207,9 +211,11 @@ namespace MAST { /*! * advances the time step and copies the current solution to old - * solution, and so on. + * solution, and so on. If \p increment_time is \p true then + * the value fo time will be incremented in System. This may be useful when calling this + * method inside a sensitivity analysis routine. */ - virtual void advance_time_step(); + virtual void advance_time_step(bool increment_time = true); /*! * advances the time step and copies the current sensitivity solution