599 *
return -Utilities::fixed_power<2, double>(this->
get_time()) + 6;
609 * Then, we handle the right boundary condition.
613 *
class RightBoundaryValues :
public Function<dim>
616 * RightBoundaryValues(
const double strike_price,
const double interest_rate);
619 *
const unsigned int component = 0)
const override;
622 *
const double strike_price;
623 *
const double interest_rate;
628 * RightBoundaryValues<dim>::RightBoundaryValues(
const double strike_price,
629 *
const double interest_rate)
630 * : strike_price(strike_price)
631 * , interest_rate(interest_rate)
636 *
double RightBoundaryValues<dim>::value(
const Point<dim> &p,
637 *
const unsigned int component)
const
640 *
return -Utilities::fixed_power<2, double>(p[component]) -
641 * Utilities::fixed_power<2, double>(this->
get_time()) + 6;
643 *
return (p[component] - strike_price) *
652 * Finally, we handle the right hand side.
656 *
class RightHandSide :
public Function<dim>
659 * RightHandSide(
const double asset_volatility,
const double interest_rate);
662 *
const unsigned int component = 0)
const override;
665 *
const double asset_volatility;
666 *
const double interest_rate;
671 * RightHandSide<dim>::RightHandSide(
const double asset_volatility,
672 *
const double interest_rate)
673 * : asset_volatility(asset_volatility)
674 * , interest_rate(interest_rate)
679 *
double RightHandSide<dim>::value(
const Point<dim> &p,
680 *
const unsigned int component)
const
684 * Utilities::fixed_power<2, double>(asset_volatility * p[component]) -
685 * 2 * interest_rate * Utilities::fixed_power<2, double>(p[component]) -
687 * (-Utilities::fixed_power<2, double>(p[component]) -
688 * Utilities::fixed_power<2, double>(this->
get_time()) + 6);
701 * <a name=
"step_78-ThecodeBlackScholescodeClass"></a>
702 * <h3>The <code>BlackScholes</code> Class</h3>
706 * The next piece is the declaration of the main
class of this program. This
707 * is very similar to the @ref step_26
"step-26" tutorial, with some modifications. New
708 * matrices had to be added to calculate the A and B matrices, as well as the
709 * @f$V_{diff}@f$ vector mentioned in the introduction. We also define the
710 * parameters used in the problem.
714 * - <code>maximum_stock_price</code>: The imposed upper bound on the spatial
715 * domain. This is the maximum allowed stock price.
716 * - <code>maturity_time</code>: The upper bound on the time domain. This is
717 * when the option expires.\n
718 * - <code>asset_volatility</code>: The volatility of the stock price.\n
719 * - <code>interest_rate</code>: The risk free interest rate.\n
720 * - <code>strike_price</code>: The agreed upon price that the buyer will
721 * have the option of purchasing the stocks at the expiration time.
725 * Some slight differences between
this program and @ref step_26
"step-26" are the creation
726 * of the <code>a_matrix</code> and the <code>b_matrix</code>, which is
727 * described in the introduction. We then also need to store the current time,
728 * the size of the time step, and the number of the current time step.
729 * Next, we will store the output into a <code>
DataOutStack</code>
730 * variable because we will be layering the solution at each time on top of
731 * one another to create the solution manifold. Then, we have a variable that
732 * stores the current cycle and number of cycles that we will
run when
733 * calculating the solution. The cycle is one full solution calculation given
734 * a mesh. We
refine the mesh once in between each cycle to exhibit the
735 * convergence properties of our program. Finally, we store the convergence
736 * data into a convergence table.
740 * As far as member
functions are concerned, we have a function that
741 * calculates the convergence information
for each cycle, called
742 * <code>process_solution</code>. This is just like what is done in @ref step_7
"step-7".
754 *
void setup_system();
755 *
void solve_time_step();
756 *
void refine_grid();
757 *
void process_solution();
758 *
void add_results_for_output();
759 *
void write_convergence_table();
761 *
const double maximum_stock_price;
762 *
const double maturity_time;
763 *
const double asset_volatility;
764 *
const double interest_rate;
765 *
const double strike_price;
786 *
const double theta;
787 *
const unsigned int n_cycles;
788 *
const unsigned int n_time_steps;
791 * std::vector<std::string> solution_names;
799 * <a name=
"step_78-ThecodeBlackScholescodeImplementation"></a>
800 * <h3>The <code>BlackScholes</code> Implementation</h3>
804 * Now, we get to the implementation of the main
class. We will set the
values
805 *
for the various parameters used in the problem. These were chosen because
806 * they are fairly normal
values for these parameters. Although the stock
807 * price has no upper bound in reality (it is in fact infinite), we impose
808 * an upper bound that is twice the strike price. This is a somewhat arbitrary
809 * choice to be twice the strike price, but it is large enough to see the
810 * interesting parts of the solution.
814 * BlackScholes<dim>::BlackScholes()
815 * : maximum_stock_price(1.)
816 * , maturity_time(1.)
817 * , asset_volatility(.2)
818 * , interest_rate(0.05)
819 * , strike_price(0.5)
825 * , n_time_steps(5000)
827 *
Assert(dim == 1, ExcNotImplemented());
833 * <a name=
"step_78-codeBlackScholessetup_systemcode"></a>
834 * <h4><code>BlackScholes::setup_system</code></h4>
838 * The next function sets up the
DoFHandler object, computes
839 * the constraints, and sets the linear algebra objects to their correct
840 * sizes. We also compute the @ref GlossMassMatrix
"mass matrix" here by calling a function from the
841 * library. We will compute the other 3 matrices next, because these need to
842 * be computed
'by hand'.
846 * Note, that the time step is initialized here because the maturity time was
847 * needed to compute the time step.
851 *
void BlackScholes<dim>::setup_system()
855 * time_step = maturity_time / n_time_steps;
857 * constraints.clear();
859 * constraints.close();
865 * sparsity_pattern.copy_from(dsp);
868 * laplace_matrix.reinit(sparsity_pattern);
869 * a_matrix.reinit(sparsity_pattern);
870 * b_matrix.reinit(sparsity_pattern);
871 * system_matrix.reinit(sparsity_pattern);
879 * Below is the code to create the Laplace
matrix with non-
constant
880 * coefficients. This corresponds to the
matrix D in the introduction. This
881 * non-
constant coefficient is represented in the
882 * <code>current_coefficient</code> variable.
885 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
887 *
const QGauss<dim> quadrature_formula(fe.degree + 1);
889 * quadrature_formula,
892 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
893 *
for (
const auto &cell : dof_handler.active_cell_iterators())
896 * fe_values.reinit(cell);
897 *
for (
const unsigned int q_index : fe_values.quadrature_point_indices())
899 * const double current_coefficient =
900 * fe_values.quadrature_point(q_index).square();
901 *
for (
const unsigned int i : fe_values.dof_indices())
903 * for (const unsigned
int j : fe_values.dof_indices())
905 * (current_coefficient *
906 * fe_values.shape_grad(i, q_index) *
907 * fe_values.shape_grad(j, q_index) *
908 * fe_values.JxW(q_index));
911 * cell->get_dof_indices(local_dof_indices);
912 *
for (
const unsigned int i : fe_values.dof_indices())
914 * for (const unsigned
int j : fe_values.dof_indices())
915 * laplace_matrix.add(local_dof_indices[i],
916 * local_dof_indices[j],
923 * Now we will create the A
matrix. Below is the code to create the
matrix A
924 * as discussed in the introduction. The non
constant coefficient is again
925 * represented in the <code>current_coefficient</code> variable.
928 *
for (
const auto &cell : dof_handler.active_cell_iterators())
931 * fe_values.reinit(cell);
932 *
for (
const unsigned int q_index : fe_values.quadrature_point_indices())
934 * const
Tensor<1, dim> current_coefficient =
935 * fe_values.quadrature_point(q_index);
936 *
for (
const unsigned int i : fe_values.dof_indices())
938 * for (const unsigned
int j : fe_values.dof_indices())
941 * (current_coefficient *
942 * fe_values.shape_grad(i, q_index) *
943 * fe_values.shape_value(j, q_index) *
944 * fe_values.JxW(q_index));
948 * cell->get_dof_indices(local_dof_indices);
949 *
for (
const unsigned int i : fe_values.dof_indices())
951 * for (const unsigned
int j : fe_values.dof_indices())
952 * a_matrix.add(local_dof_indices[i],
953 * local_dof_indices[j],
960 * Finally we will create the
matrix B. Below is the code to create the
961 *
matrix B as discussed in the introduction. The non
constant coefficient
962 * is again represented in the <code>current_coefficient</code> variable.
965 *
for (
const auto &cell : dof_handler.active_cell_iterators())
968 * fe_values.reinit(cell);
969 *
for (
const unsigned int q_index : fe_values.quadrature_point_indices())
971 * const
Tensor<1, dim> current_coefficient =
972 * fe_values.quadrature_point(q_index);
973 *
for (
const unsigned int i : fe_values.dof_indices())
975 * for (const unsigned
int j : fe_values.dof_indices())
977 * (current_coefficient *
978 * fe_values.shape_value(i, q_index) *
979 * fe_values.shape_grad(j, q_index) *
980 * fe_values.JxW(q_index));
983 * cell->get_dof_indices(local_dof_indices);
984 *
for (
const unsigned int i : fe_values.dof_indices())
986 * for (const unsigned
int j : fe_values.dof_indices())
987 * b_matrix.add(local_dof_indices[i],
988 * local_dof_indices[j],
993 * solution.reinit(dof_handler.n_dofs());
994 * system_rhs.reinit(dof_handler.n_dofs());
1000 * <a name=
"step_78-codeBlackScholessolve_time_stepcode"></a>
1001 * <h4><code>BlackScholes::solve_time_step</code></h4>
1005 * The next function is the one that solves the actual linear system
for a
1006 * single time step. The only interesting thing here is that the matrices
1011 *
template <
int dim>
1012 *
void BlackScholes<dim>::solve_time_step()
1017 * preconditioner.
initialize(system_matrix, 1.0);
1018 * cg.solve(system_matrix, solution, system_rhs, preconditioner);
1019 * constraints.distribute(solution);
1025 * <a name=
"step_78-codeBlackScholesadd_results_for_outputcode"></a>
1026 * <h4><code>BlackScholes::add_results_for_output</code></h4>
1030 * This is simply the function to stitch the solution pieces together. For
1031 *
this, we create a
new layer at each time, and then add the solution vector
1032 *
for that timestep. The function then stitches
this together with the old
1033 * solutions
using 'build_patches
'.
1036 * template <int dim>
1037 * void BlackScholes<dim>::add_results_for_output()
1039 * data_out_stack.new_parameter_value(time, time_step);
1040 * data_out_stack.attach_dof_handler(dof_handler);
1041 * data_out_stack.add_data_vector(solution, solution_names);
1042 * data_out_stack.build_patches(2);
1043 * data_out_stack.finish_parameter_value();
1049 * <a name="step_78-codeBlackScholesrefine_gridcode"></a>
1050 * <h4><code>BlackScholes::refine_grid</code></h4>
1054 * It is somewhat unnecessary to have a function for the global refinement
1055 * that we do. The reason for the function is to allow for the possibility of
1056 * an adaptive refinement later.
1059 * template <int dim>
1060 * void BlackScholes<dim>::refine_grid()
1062 * triangulation.refine_global(1);
1068 * <a name="step_78-codeBlackScholesprocess_solutioncode"></a>
1069 * <h4><code>BlackScholes::process_solution</code></h4>
1073 * This is where we calculate the convergence and error data to evaluate the
1074 * effectiveness of the program. Here, we calculate the @f$L^2@f$, @f$H^1@f$ and
1075 * @f$L^{\infty}@f$ norms.
1078 * template <int dim>
1079 * void BlackScholes<dim>::process_solution()
1081 * Solution<dim> sol(maturity_time);
1082 * sol.set_time(time);
1083 * Vector<float> difference_per_cell(triangulation.n_active_cells());
1084 * VectorTools::integrate_difference(dof_handler,
1087 * difference_per_cell,
1088 * QGauss<dim>(fe.degree + 1),
1089 * VectorTools::L2_norm);
1090 * const double L2_error =
1091 * VectorTools::compute_global_error(triangulation,
1092 * difference_per_cell,
1093 * VectorTools::L2_norm);
1094 * VectorTools::integrate_difference(dof_handler,
1097 * difference_per_cell,
1098 * QGauss<dim>(fe.degree + 1),
1099 * VectorTools::H1_seminorm);
1100 * const double H1_error =
1101 * VectorTools::compute_global_error(triangulation,
1102 * difference_per_cell,
1103 * VectorTools::H1_seminorm);
1104 * const QTrapezoid<1> q_trapezoid;
1105 * const QIterated<dim> q_iterated(q_trapezoid, fe.degree * 2 + 1);
1106 * VectorTools::integrate_difference(dof_handler,
1109 * difference_per_cell,
1111 * VectorTools::Linfty_norm);
1112 * const double Linfty_error =
1113 * VectorTools::compute_global_error(triangulation,
1114 * difference_per_cell,
1115 * VectorTools::Linfty_norm);
1116 * const unsigned int n_active_cells = triangulation.n_active_cells();
1117 * const unsigned int n_dofs = dof_handler.n_dofs();
1118 * convergence_table.add_value("cells", n_active_cells);
1119 * convergence_table.add_value("dofs", n_dofs);
1120 * convergence_table.add_value("L2", L2_error);
1121 * convergence_table.add_value("H1", H1_error);
1122 * convergence_table.add_value("Linfty", Linfty_error);
1128 * <a name="step_78-codeBlackScholeswrite_convergence_tablecode"></a>
1129 * <h4><code>BlackScholes::write_convergence_table</code> </h4>
1133 * This next part is building the convergence and error tables. By this, we
1134 * need to set the settings for how to output the data that was calculated
1135 * during <code>BlackScholes::process_solution</code>. First, we will create
1136 * the headings and set up the cells properly. During this, we will also
1137 * prescribe the precision of our results. Then we will write the calculated
1138 * errors based on the @f$L^2@f$, @f$H^1@f$, and @f$L^{\infty}@f$ norms to the console and
1139 * to the error LaTeX file.
1142 * template <int dim>
1143 * void BlackScholes<dim>::write_convergence_table()
1145 * convergence_table.set_precision("L2", 3);
1146 * convergence_table.set_precision("H1", 3);
1147 * convergence_table.set_precision("Linfty", 3);
1148 * convergence_table.set_scientific("L2", true);
1149 * convergence_table.set_scientific("H1", true);
1150 * convergence_table.set_scientific("Linfty", true);
1151 * convergence_table.set_tex_caption("cells", "\\# cells");
1152 * convergence_table.set_tex_caption("dofs", "\\# dofs");
1153 * convergence_table.set_tex_caption("L2", "@fL^2@f-error");
1154 * convergence_table.set_tex_caption("H1", "@fH^1@f-error");
1155 * convergence_table.set_tex_caption("Linfty", "@fL^\\infty@f-error");
1156 * convergence_table.set_tex_format("cells", "r");
1157 * convergence_table.set_tex_format("dofs", "r");
1158 * std::cout << std::endl;
1159 * convergence_table.write_text(std::cout);
1160 * std::string error_filename = "error";
1161 * error_filename += "-global";
1162 * error_filename += ".tex";
1163 * std::ofstream error_table_file(error_filename);
1164 * convergence_table.write_tex(error_table_file);
1168 * Next, we will make the convergence table. We will again write this to
1169 * the console and to the convergence LaTeX file.
1172 * convergence_table.add_column_to_supercolumn("cells", "n cells");
1173 * std::vector<std::string> new_order;
1174 * new_order.emplace_back("n cells");
1175 * new_order.emplace_back("H1");
1176 * new_order.emplace_back("L2");
1177 * convergence_table.set_column_order(new_order);
1178 * convergence_table.evaluate_convergence_rates(
1179 * "L2", ConvergenceTable::reduction_rate);
1180 * convergence_table.evaluate_convergence_rates(
1181 * "L2", ConvergenceTable::reduction_rate_log2);
1182 * convergence_table.evaluate_convergence_rates(
1183 * "H1", ConvergenceTable::reduction_rate);
1184 * convergence_table.evaluate_convergence_rates(
1185 * "H1", ConvergenceTable::reduction_rate_log2);
1186 * std::cout << std::endl;
1187 * convergence_table.write_text(std::cout);
1188 * std::string conv_filename = "convergence";
1189 * conv_filename += "-global";
1190 * switch (fe.degree)
1193 * conv_filename += "-q1";
1196 * conv_filename += "-q2";
1199 * DEAL_II_NOT_IMPLEMENTED();
1201 * conv_filename += ".tex";
1202 * std::ofstream table_file(conv_filename);
1203 * convergence_table.write_tex(table_file);
1209 * <a name="step_78-codeBlackScholesruncode"></a>
1210 * <h4><code>BlackScholes::run</code></h4>
1214 * Now we get to the main driver of the program. This is where we do all the
1215 * work of looping through the time steps and calculating the solution vector
1216 * each time. Here at the top, we set the initial refinement value and then
1217 * create a mesh. Then we refine this mesh once. Next, we set up the
1218 * data_out_stack object to store our solution. Finally, we start a for loop
1219 * to loop through the cycles. This lets us recalculate a solution for each
1220 * successive mesh refinement. At the beginning of each iteration, we need to
1221 * reset the time and time step number. We introduce an if statement to
1222 * accomplish this because we don't want to
do this on the
first iteration.
1225 *
template <
int dim>
1226 *
void BlackScholes<dim>::run()
1231 * solution_names.emplace_back(
"u");
1232 * data_out_stack.declare_data_vector(solution_names,
1238 *
for (
unsigned int cycle = 0; cycle < n_cycles; ++cycle)
1248 * std::cout << std::endl
1249 * <<
"===========================================" << std::endl
1250 * <<
"Cycle " << cycle <<
':' << std::endl
1251 * <<
"Number of active cells: "
1253 * <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
1258 * InitialConditions<dim>(strike_price),
1261 *
if (cycle == (n_cycles - 1))
1263 * add_results_for_output();
1268 * Next, we
run the main
loop which runs until we exceed the maturity
1269 * time. We
first compute the right hand side of the equation, which is
1270 * described in the introduction. Recall that it contains the term
1271 * @f$\left[-\frac{1}{4}k_n\sigma^2\mathbf{D}-k_nr\mathbf{M}+k_n\sigma^2
1272 * \mathbf{B}-k_nr\mathbf{A}+\mathbf{M}\right]V^{n-1}@f$. We put these
1273 * terms into the variable system_rhs, with the help of a temporary
1277 * vmult_result.reinit(dof_handler.n_dofs());
1278 * forcing_terms.reinit(dof_handler.n_dofs());
1279 *
for (
unsigned int timestep_number = 0; timestep_number < n_time_steps;
1280 * ++timestep_number)
1282 * time += time_step;
1284 *
if (timestep_number % 1000 == 0)
1285 * std::cout <<
"Time step " << timestep_number <<
" at t=" << time
1290 * laplace_matrix.vmult(vmult_result, solution);
1292 * (-1) * (1 - theta) * time_step *
1293 * Utilities::fixed_power<2, double>(asset_volatility) * 0.5,
1297 * system_rhs.add((-1) * (1 - theta) * time_step * interest_rate * 2,
1300 * a_matrix.vmult(vmult_result, solution);
1301 * system_rhs.add((-1) * time_step * interest_rate, vmult_result);
1303 * b_matrix.vmult(vmult_result, solution);
1305 * (-1) * Utilities::fixed_power<2, double>(asset_volatility) *
1311 * The
second piece is to compute the contributions of the source
1312 * terms. This corresponds to the term @f$-k_n\left[\frac{1}{2}
F^{n-1}
1313 * +\frac{1}{2}
F^n\right]@f$. The following code calls
1315 * where we set the time of the right hand side (source) function
1316 * before we evaluate it. The result of
this all ends up in the
1317 * forcing_terms variable:
1320 * RightHandSide<dim> rhs_function(asset_volatility, interest_rate);
1321 * rhs_function.set_time(time);
1326 * forcing_terms *= time_step * theta;
1327 * system_rhs -= forcing_terms;
1329 * rhs_function.set_time(time - time_step);
1334 * forcing_terms *= time_step * (1 - theta);
1335 * system_rhs -= forcing_terms;
1339 * Next, we add the forcing terms to the ones that come from the
1340 * time stepping, and also build the
matrix @f$\left[\mathbf{M}+
1341 * \frac{1}{4}k_n\sigma^2\mathbf{D}+k_nr\mathbf{M}\right]@f$ that we
1342 * have to
invert in each time step. The
final piece of these
1343 * operations is to eliminate hanging node constrained degrees of
1344 * freedom from the linear system:
1347 * system_matrix.copy_from(mass_matrix);
1348 * system_matrix.add(
1349 * (theta)*time_step *
1350 * Utilities::fixed_power<2, double>(asset_volatility) * 0.5,
1352 * system_matrix.add((time_step)*interest_rate * theta * (1 + 1),
1355 * constraints.condense(system_matrix, system_rhs);
1359 * There is one more operation we need to
do before we can solve it:
1360 * boundary
values. To
this end, we create a boundary
value object,
1361 * set the proper time to the one of the current time step, and
1362 * evaluate it as we have done many times before. The result is used
1363 * to also set the correct boundary
values in the linear system:
1367 * RightBoundaryValues<dim> right_boundary_function(strike_price,
1369 * LeftBoundaryValues<dim> left_boundary_function;
1370 * right_boundary_function.set_time(time);
1371 * left_boundary_function.set_time(time);
1372 * std::map<types::global_dof_index, double> boundary_values;
1375 * left_boundary_function,
1379 * right_boundary_function,
1389 * With
this out of the way, all we have to
do is solve the system,
1390 * generate graphical data on the last cycle, and create the
1391 * convergence table data.
1394 * solve_time_step();
1396 *
if (cycle == (n_cycles - 1))
1398 * add_results_for_output();
1402 * process_solution();
1406 *
const std::string filename =
"solution.vtk";
1407 * std::ofstream output(filename);
1408 * data_out_stack.write_vtk(output);
1411 * write_convergence_table();
1420 * <a name=
"step_78-ThecodemaincodeFunction"></a>
1421 * <h3>The <code>main</code>
Function</h3>
1425 * Having made it
this far, there is, again,
nothing much to discuss
for the
1426 * main function of
this program: it looks like all such
functions since @ref step_6
"step-6".
1433 *
using namespace BlackScholesSolver;
1435 * BlackScholes<1> black_scholes_solver;
1436 * black_scholes_solver.run();
1438 *
catch (std::exception &exc)
1440 * std::cerr << std::endl
1442 * <<
"----------------------------------------------------"
1444 * std::cerr <<
"Exception on processing: " << std::endl
1445 * << exc.what() << std::endl
1446 * <<
"Aborting!" << std::endl
1447 * <<
"----------------------------------------------------"
1453 * std::cerr << std::endl
1455 * <<
"----------------------------------------------------"
1457 * std::cerr <<
"Unknown exception!" << std::endl
1458 * <<
"Aborting!" << std::endl
1459 * <<
"----------------------------------------------------"
1466<a name=
"step_78-Results"></a><h1>Results</h1>
1470Below is the output of the program:
1472===========================================
1473Number of active cells: 1
1474Number of degrees of freedom: 2
1476Time step 0 at t=0.0002
1480Number of active cells: 128
1481Number of degrees of freedom: 129
1483Time step 0 at t=0.0002
1484Time step 1000 at t=0.2002
1485Time step 2000 at t=0.4002
1486Time step 3000 at t=0.6002
1487Time step 4000 at t=0.8002
1489cells dofs
L2 H1 Linfty
1490 1 2 1.667e-01 5.774e-01 2.222e-01
1491 2 3 3.906e-02 2.889e-01 5.380e-02
1492 4 5 9.679e-03 1.444e-01 1.357e-02
1493 8 9 2.405e-03 7.218e-02 3.419e-03
1494 16 17 5.967e-04 3.609e-02 8.597e-04
1495 32 33 1.457e-04 1.804e-02 2.155e-04
1496 64 65 3.307e-05 9.022e-03 5.388e-05
1497 128 129 5.016e-06 4.511e-03 1.342e-05
1500 1 5.774e-01 - - 1.667e-01 - -
1501 2 2.889e-01 2.00 1.00 3.906e-02 4.27 2.09
1502 4 1.444e-01 2.00 1.00 9.679e-03 4.04 2.01
1503 8 7.218e-02 2.00 1.00 2.405e-03 4.02 2.01
1504 16 3.609e-02 2.00 1.00 5.967e-04 4.03 2.01
1505 32 1.804e-02 2.00 1.00 1.457e-04 4.10 2.03
1506 64 9.022e-03 2.00 1.00 3.307e-05 4.41 2.14
1507 128 4.511e-03 2.00 1.00 5.016e-06 6.59 2.72
1510What is more interesting is the output of the convergence tables. They are
1511outputted into the console, as well into a LaTeX file. The convergence tables
1512are shown above. Here, you can see that the solution has a convergence rate
1513of @f$\mathcal{O}(h)@f$ with respect to the @f$H^1@f$-norm, and the solution has a convergence rate
1514of @f$\mathcal{O}(h^2)@f$ with respect to the @f$L^2@f$-
norm.
1517Below is the visualization of the solution.
1519<div style=
"text-align:center;">
1520 <img src=
"https://www.dealii.org/images/steps/developer/step-78.mms-solution.png"
1521 alt=
"Solution of the MMS problem.">
1525<a name=
"step_78-PlainProg"></a>
1526<h1> The plain program</h1>
1527@include
"step-78.cc"
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
#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.
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.
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.)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
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 > &)