600 *
return -Utilities::fixed_power<2, double>(this->
get_time()) + 6;
620 *
const unsigned int component = 0)
const override;
638 *
const unsigned int component)
const
641 *
return -Utilities::fixed_power<2, double>(p[component]) -
642 *
Utilities::fixed_power<2, double>(this->
get_time()) + 6;
663 *
const unsigned int component = 0)
const override;
681 *
const unsigned int component)
const
686 *
2 *
interest_rate * Utilities::fixed_power<2, double>(p[component]) -
688 *
(-Utilities::fixed_power<2, double>(p[component]) -
689 *
Utilities::fixed_power<2, double>(this->
get_time()) + 6);
702 * <a name=
"step_78-ThecodeBlackScholescodeClass"></a>
757 *
void refine_grid();
787 *
const double theta;
788 *
const unsigned int n_cycles;
800 * <a name=
"step_78-ThecodeBlackScholescodeImplementation"></a>
834 * <a name=
"step_78-codeBlackScholessetup_systemcode"></a>
854 *
dof_handler.distribute_dofs(fe);
858 *
constraints.
clear();
860 *
constraints.close();
866 *
sparsity_pattern.copy_from(
dsp);
870 *
a_matrix.reinit(sparsity_pattern);
871 *
b_matrix.reinit(sparsity_pattern);
872 *
system_matrix.reinit(sparsity_pattern);
886 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
893 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
894 *
for (
const auto &cell : dof_handler.active_cell_iterators())
897 *
fe_values.reinit(cell);
898 *
for (
const unsigned int q_index : fe_values.quadrature_point_indices())
901 *
fe_values.quadrature_point(q_index).square();
902 *
for (
const unsigned int i : fe_values.dof_indices())
907 *
fe_values.shape_grad(i, q_index) *
908 *
fe_values.shape_grad(
j, q_index) *
909 *
fe_values.JxW(q_index));
912 *
cell->get_dof_indices(local_dof_indices);
913 *
for (
const unsigned int i : fe_values.dof_indices())
917 *
local_dof_indices[
j],
929 *
for (
const auto &cell : dof_handler.active_cell_iterators())
932 *
fe_values.reinit(cell);
933 *
for (
const unsigned int q_index : fe_values.quadrature_point_indices())
936 *
fe_values.quadrature_point(q_index);
937 *
for (
const unsigned int i : fe_values.dof_indices())
943 *
fe_values.shape_grad(i, q_index) *
944 *
fe_values.shape_value(
j, q_index) *
945 *
fe_values.JxW(q_index));
949 *
cell->get_dof_indices(local_dof_indices);
950 *
for (
const unsigned int i : fe_values.dof_indices())
953 *
a_matrix.add(local_dof_indices[i],
954 *
local_dof_indices[
j],
966 *
for (
const auto &cell : dof_handler.active_cell_iterators())
969 *
fe_values.reinit(cell);
970 *
for (
const unsigned int q_index : fe_values.quadrature_point_indices())
973 *
fe_values.quadrature_point(q_index);
974 *
for (
const unsigned int i : fe_values.dof_indices())
979 *
fe_values.shape_value(i, q_index) *
980 *
fe_values.shape_grad(
j, q_index) *
981 *
fe_values.JxW(q_index));
984 *
cell->get_dof_indices(local_dof_indices);
985 *
for (
const unsigned int i : fe_values.dof_indices())
988 *
b_matrix.add(local_dof_indices[i],
989 *
local_dof_indices[
j],
994 *
solution.reinit(dof_handler.n_dofs());
1001 * <a name=
"step_78-codeBlackScholessolve_time_stepcode"></a>
1012 *
template <
int dim>
1018 *
preconditioner.initialize(system_matrix, 1.0);
1019 *
cg.solve(system_matrix, solution,
system_rhs, preconditioner);
1020 *
constraints.distribute(solution);
1026 * <a name=
"step_78-codeBlackScholesadd_results_for_outputcode"></a>
1037 * template <int dim>
1038 * void BlackScholes<dim>::add_results_for_output()
1040 * data_out_stack.new_parameter_value(time, time_step);
1041 * data_out_stack.attach_dof_handler(dof_handler);
1042 * data_out_stack.add_data_vector(solution, solution_names);
1043 * data_out_stack.build_patches(2);
1044 * data_out_stack.finish_parameter_value();
1050 * <a name="step_78-codeBlackScholesrefine_gridcode"></a>
1051 * <h4><code>BlackScholes::refine_grid</code></h4>
1055 * It is somewhat unnecessary to have a function for the global refinement
1056 * that we do. The reason for the function is to allow for the possibility of
1057 * an adaptive refinement later.
1060 * template <int dim>
1061 * void BlackScholes<dim>::refine_grid()
1063 * triangulation.refine_global(1);
1069 * <a name="step_78-codeBlackScholesprocess_solutioncode"></a>
1070 * <h4><code>BlackScholes::process_solution</code></h4>
1074 * This is where we calculate the convergence and error data to evaluate the
1075 * effectiveness of the program. Here, we calculate the @f$L^2@f$, @f$H^1@f$ and
1076 * @f$L^{\infty}@f$ norms.
1079 * template <int dim>
1080 * void BlackScholes<dim>::process_solution()
1082 * Solution<dim> sol(maturity_time);
1083 * sol.set_time(time);
1084 * Vector<float> difference_per_cell(triangulation.n_active_cells());
1085 * VectorTools::integrate_difference(dof_handler,
1088 * difference_per_cell,
1089 * QGauss<dim>(fe.degree + 1),
1090 * VectorTools::L2_norm);
1091 * const double L2_error =
1092 * VectorTools::compute_global_error(triangulation,
1093 * difference_per_cell,
1094 * VectorTools::L2_norm);
1095 * VectorTools::integrate_difference(dof_handler,
1098 * difference_per_cell,
1099 * QGauss<dim>(fe.degree + 1),
1100 * VectorTools::H1_seminorm);
1101 * const double H1_error =
1102 * VectorTools::compute_global_error(triangulation,
1103 * difference_per_cell,
1104 * VectorTools::H1_seminorm);
1105 * const QTrapezoid<1> q_trapezoid;
1106 * const QIterated<dim> q_iterated(q_trapezoid, fe.degree * 2 + 1);
1107 * VectorTools::integrate_difference(dof_handler,
1110 * difference_per_cell,
1112 * VectorTools::Linfty_norm);
1113 * const double Linfty_error =
1114 * VectorTools::compute_global_error(triangulation,
1115 * difference_per_cell,
1116 * VectorTools::Linfty_norm);
1117 * const unsigned int n_active_cells = triangulation.n_active_cells();
1118 * const unsigned int n_dofs = dof_handler.n_dofs();
1119 * convergence_table.add_value("cells", n_active_cells);
1120 * convergence_table.add_value("dofs", n_dofs);
1121 * convergence_table.add_value("L2", L2_error);
1122 * convergence_table.add_value("H1", H1_error);
1123 * convergence_table.add_value("Linfty", Linfty_error);
1129 * <a name="step_78-codeBlackScholeswrite_convergence_tablecode"></a>
1130 * <h4><code>BlackScholes::write_convergence_table</code> </h4>
1134 * This next part is building the convergence and error tables. By this, we
1135 * need to set the settings for how to output the data that was calculated
1136 * during <code>BlackScholes::process_solution</code>. First, we will create
1137 * the headings and set up the cells properly. During this, we will also
1138 * prescribe the precision of our results. Then we will write the calculated
1139 * errors based on the @f$L^2@f$, @f$H^1@f$, and @f$L^{\infty}@f$ norms to the console and
1140 * to the error LaTeX file.
1143 * template <int dim>
1144 * void BlackScholes<dim>::write_convergence_table()
1146 * convergence_table.set_precision("L2", 3);
1147 * convergence_table.set_precision("H1", 3);
1148 * convergence_table.set_precision("Linfty", 3);
1149 * convergence_table.set_scientific("L2", true);
1150 * convergence_table.set_scientific("H1", true);
1151 * convergence_table.set_scientific("Linfty", true);
1152 * convergence_table.set_tex_caption("cells", "\\# cells");
1153 * convergence_table.set_tex_caption("dofs", "\\# dofs");
1154 * convergence_table.set_tex_caption("L2", "@fL^2@f-error");
1155 * convergence_table.set_tex_caption("H1", "@fH^1@f-error");
1156 * convergence_table.set_tex_caption("Linfty", "@fL^\\infty@f-error");
1157 * convergence_table.set_tex_format("cells", "r");
1158 * convergence_table.set_tex_format("dofs", "r");
1159 * std::cout << std::endl;
1160 * convergence_table.write_text(std::cout);
1161 * std::string error_filename = "error";
1162 * error_filename += "-global";
1163 * error_filename += ".tex";
1164 * std::ofstream error_table_file(error_filename);
1165 * convergence_table.write_tex(error_table_file);
1169 * Next, we will make the convergence table. We will again write this to
1170 * the console and to the convergence LaTeX file.
1173 * convergence_table.add_column_to_supercolumn("cells", "n cells");
1174 * std::vector<std::string> new_order;
1175 * new_order.emplace_back("n cells");
1176 * new_order.emplace_back("H1");
1177 * new_order.emplace_back("L2");
1178 * convergence_table.set_column_order(new_order);
1179 * convergence_table.evaluate_convergence_rates(
1180 * "L2", ConvergenceTable::reduction_rate);
1181 * convergence_table.evaluate_convergence_rates(
1182 * "L2", ConvergenceTable::reduction_rate_log2);
1183 * convergence_table.evaluate_convergence_rates(
1184 * "H1", ConvergenceTable::reduction_rate);
1185 * convergence_table.evaluate_convergence_rates(
1186 * "H1", ConvergenceTable::reduction_rate_log2);
1187 * std::cout << std::endl;
1188 * convergence_table.write_text(std::cout);
1189 * std::string conv_filename = "convergence";
1190 * conv_filename += "-global";
1191 * switch (fe.degree)
1194 * conv_filename += "-q1";
1197 * conv_filename += "-q2";
1200 * DEAL_II_NOT_IMPLEMENTED();
1202 * conv_filename += ".tex";
1203 * std::ofstream table_file(conv_filename);
1204 * convergence_table.write_tex(table_file);
1210 * <a name="step_78-codeBlackScholesruncode"></a>
1211 * <h4><code>BlackScholes::run</code></h4>
1215 * Now we get to the main driver of the program. This is where we do all the
1216 * work of looping through the time steps and calculating the solution vector
1217 * each time. Here at the top, we set the initial refinement value and then
1218 * create a mesh. Then we refine this mesh once. Next, we set up the
1219 * data_out_stack object to store our solution. Finally, we start a for loop
1220 * to loop through the cycles. This lets us recalculate a solution for each
1221 * successive mesh refinement. At the beginning of each iteration, we need to
1222 * reset the time and time step number. We introduce an if statement to
1226 *
template <
int dim>
1239 *
for (
unsigned int cycle = 0; cycle < n_cycles; ++cycle)
1249 *
std::cout << std::endl
1250 *
<<
"===========================================" << std::endl
1251 *
<<
"Cycle " << cycle <<
':' << std::endl
1252 *
<<
"Number of active cells: "
1254 *
<<
"Number of degrees of freedom: " << dof_handler.n_dofs()
1262 *
if (cycle == (n_cycles - 1))
1348 *
system_matrix.copy_from(mass_matrix);
1349 *
system_matrix.add(
1356 *
constraints.condense(system_matrix,
system_rhs);
1397 *
if (cycle == (n_cycles - 1))
1407 *
const std::string
filename =
"solution.vtk";
1421 * <a name=
"step_78-ThecodemaincodeFunction"></a>
1439 *
catch (std::exception &exc)
1441 *
std::cerr << std::endl
1443 *
<<
"----------------------------------------------------"
1445 *
std::cerr <<
"Exception on processing: " << std::endl
1446 *
<< exc.what() << std::endl
1447 *
<<
"Aborting!" << std::endl
1448 *
<<
"----------------------------------------------------"
1454 *
std::cerr << std::endl
1456 *
<<
"----------------------------------------------------"
1458 *
std::cerr <<
"Unknown exception!" << std::endl
1459 *
<<
"Aborting!" << std::endl
1460 *
<<
"----------------------------------------------------"
1473===========================================
1477Time step 0 at t=0.0002
1484Time step 0 at t=0.0002
1485Time step 1000 at t=0.2002
1486Time step 2000 at t=0.4002
1487Time step 3000 at t=0.6002
1488Time step 4000 at t=0.8002
1491 1 2 1.667e-01 5.774e-01 2.222e-01
1492 2 3 3.906e-02 2.889e-01 5.380e-02
1493 4 5 9.679e-03 1.444e-01 1.357e-02
1494 8 9 2.405e-03 7.218e-02 3.419e-03
1495 16 17 5.967e-04 3.609e-02 8.597e-04
1496 32 33 1.457e-04 1.804e-02 2.155e-04
1497 64 65 3.307e-05 9.022e-03 5.388e-05
1498 128 129 5.016e-06 4.511e-03 1.342e-05
1501 1 5.774e-01 - - 1.667e-01 - -
1502 2 2.889e-01 2.00 1.00 3.906e-02 4.27 2.09
1503 4 1.444e-01 2.00 1.00 9.679e-03 4.04 2.01
1504 8 7.218e-02 2.00 1.00 2.405e-03 4.02 2.01
1505 16 3.609e-02 2.00 1.00 5.967e-04 4.03 2.01
1506 32 1.804e-02 2.00 1.00 1.457e-04 4.10 2.03
1507 64 9.022e-03 2.00 1.00 3.307e-05 4.41 2.14
1508 128 4.511e-03 2.00 1.00 5.016e-06 6.59 2.72
1521 <
img src=
"https://www.dealii.org/images/steps/developer/step-78.mms-solution.png"
1522 alt=
"Solution of the MMS problem.">
1526<a name=
"step_78-PlainProg"></a>
numbers::NumberTraits< Number >::real_type norm() const
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
constexpr types::blas_int one
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
void create_mass_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, MatrixType &matrix, const Function< spacedim, typename MatrixType::value_type > *const a=nullptr, const AffineConstraints< typename MatrixType::value_type > &constraints=AffineConstraints< typename MatrixType::value_type >())
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * end(VectorType &V)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
int(&) functions(const void *v1, const void *v2)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)