/* This plugin is written for perturb the initial temperature field when restarting a model. This implementation is based on the prescribed_velocity.cc file in aspect/cookbooks/prescribed_velocity/ @written : Chameera Silva {2023-03-19} */ #include #include #include #include #include #include namespace aspect { using namespace dealii; // Global variables (to be set by parameters) bool prescribe_internal_temperatures; // Because we do not initially know what dimensions we're in, we need // function parser objects for both 2d and 3d Functions::ParsedFunction<2> prescribed_temperature_indicator_function_2d (2); Functions::ParsedFunction<3> prescribed_temperature_indicator_function_3d (3); Functions::ParsedFunction<2> prescribed_temperature_function_2d (2); Functions::ParsedFunction<3> prescribed_temperature_function_3d (3); /** * Declare additinal parameters. */ void declare_parameters(const unsigned int dim, ParameterHandler &prm) { prm.declare_entry ("Prescribe internal temperatures", "false", Patterns::Bool (), "Whether or not to use any prescribed internal temperatures. " "Locations in which to prescribe temperatures are defined " "in section ``Prescribed temperatures/Indicator function'' " "and the temperatures are defined in section ``Prescribed " "temperatures/Temperature function''. Indicators are evaluated " "at the center of each cell, and all DOFs associated with " "the specified temperature component at the indicated cells " "are constrained." ); prm.enter_subsection ("Prescribed temperatures"); { prm.enter_subsection ("Indicator fucntion"); { if (dim == 2) Functions::ParsedFunction<2>::declare_parameters (prm, 2); else Functions::ParsedFunction<3>::declare_parameters (prm, 3); } prm.leave_subsection (); prm.enter_subsection ("Temperature function"); { if (dim == 2) Functions::ParsedFunction<2>::declare_parameters (prm, 2); else Functions::ParsedFunction<3>::declare_parameters (prm, 3); } prm.leave_subsection (); } prm.leave_subsection (); } template void parse_parameters(const Parameters &, ParameterHandler &prm) { prescribe_internal_temperatures = prm.get_bool ("Prescribe internal temperatures"); prm.enter_subsection ("Prescribed temperatures"); { prm.enter_subsection ("Indicator function"); { try { if (dim == 2) prescribed_temperature_indicator_function_2d.parse_parameters (prm); else prescribed_temperature_indicator_function_3d.parse_parameters (prm); } catch (...) { std::cerr << "ERROR: FunctionParser failed to parse\n" << "\t'Prescribed temperatures.indicator function'\n" << "with expression\n" << "\t'" < as_2d(const Point<3> &/*p*/) { return Point<2>(); } const Point<2> &as_2d(const Point<2> &p) { return p; } const Point<3> as_3d(const Point<2> &/*p*/) { return Point<3>(); } const Point<3> &as_3d(const Point<3> &p) { return p; } } /** * This function is called by a signal which is triggered after the other constraints * have been calculated. This enables us to define additional constraints in the mass * matrix on any arbitrary degree of freedom in the model space. */ template void constrain_internal_temperatures (const SimulatorAccess &simulator_access, AffineConstraints ¤t_constraints) { if (prescribe_internal_temperatures) { const std::vector> points = aspect::Utilities::get_unit_support_points(simulator_access); const Quadrature quadrature (points); FEValues fe_values (simulator_access.get_fe(), quadrature, update_quadrature_points); typename DoFHandler::active_cell_iterator cell; // Loop over all cells for (cell = simulator_access.get_dof_handler().begin_active(); cell != simulator_access.get_dof_handler().end(); ++cell) if (! cell->is_artificial()) { fe_values.reinit (cell); std::vector local_dof_indices(simulator_access.get_fe().dofs_per_cell); cell->get_dof_indices (local_dof_indices); for (unsigned int q=0; q= simulator_access.introspection().component_indices.temperature) && (c_idx <= simulator_access.introspection().component_indices.temperature)) { // Which temperature component is this DOF associated with? const unsigned int component_direction = (c_idx - simulator_access.introspection().component_indices.temperature); // we get time passed as seconds (always) but may want // to reinterpret it in years double time = simulator_access.get_time(); if (simulator_access.convert_output_to_years()) time /= year_in_seconds; prescribed_temperature_indicator_function_2d.set_time (time); prescribed_temperature_indicator_function_3d.set_time (time); prescribed_temperature_function_2d.set_time (time); prescribed_temperature_function_3d.set_time (time); const Point p =fe_values.quadrature_point(q); // Because we defined and parsed our parameter // file differently for 2d and 3d we need to // be sure to query the correct object for // function values. The function parser // objects expect points of a certain // dimension, but Point p will be compiled for // both 2d and 3d, so we need to do some trickery // to make this compile. double indicator, u_i; if (dim == 2) { indicator = prescribed_temperature_indicator_function_2d.value (as_2d(p), component_direction); u_i = prescribed_temperature_function_2d.value (as_2d(p), component_direction); } else { indicator = prescribed_temperature_indicator_function_3d.value (as_3d(p), component_direction); u_i = prescribed_temperature_function_3d.value (as_3d(p), component_direction); } if (indicator > 0.5) { // Add a constraint of the form dof[q] = u_i // to the list of constraints. current_constraints.add_line (local_dof_indices[q]); current_constraints.set_inhomogeneity (local_dof_indices[q], u_i); } } } } } } // Connect declare_parameters and parse_parameters to appropriate signals. void parameter_connector () { SimulatorSignals<2>::declare_additional_parameters.connect (&declare_parameters); SimulatorSignals<3>::declare_additional_parameters.connect (&declare_parameters); SimulatorSignals<2>::parse_additional_parameters.connect (&parse_parameters<2>); SimulatorSignals<3>::parse_additional_parameters.connect (&parse_parameters<3>); } // Connect constraints function to correct signal. template void signal_connector (SimulatorSignals &signals) { signals.post_constraints_creation.connect (&constrain_internal_temperatures); } // Tell ASPECT to send signals to the connector functions ASPECT_REGISTER_SIGNALS_PARAMETER_CONNECTOR(parameter_connector) ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, signal_connector<3>) }