Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-78.h
Go to the documentation of this file.
1
595) const
596 *   {
597 *   #ifdef MMS
598 *   return -Utilities::fixed_power<2, double>(this->get_time()) + 6;
599 *   #else
600 *   return 0.;
601 *   #endif
602 *   }
603 *  
604 *  
605 *  
606 * @endcode
607 *
608 * Then, we handle the right boundary condition.
609 *
610 * @code
611 *   template <int dim>
612 *   class RightBoundaryValues : public Function<dim>
613 *   {
614 *   public:
615 *   RightBoundaryValues(const double strike_price, const double interest_rate);
616 *  
617 *   virtual double value(const Point<dim> &p,
618 *   const unsigned int component = 0) const override;
619 *  
620 *   private:
621 *   const double strike_price;
622 *   const double interest_rate;
623 *   };
624 *  
625 *  
626 *   template <int dim>
627 *   RightBoundaryValues<dim>::RightBoundaryValues(const double strike_price,
628 *   const double interest_rate)
629 *   : strike_price(strike_price)
630 *   , interest_rate(interest_rate)
631 *   {}
632 *  
633 *  
634 *   template <int dim>
635 *   double RightBoundaryValues<dim>::value(const Point<dim> &p,
636 *   const unsigned int component) const
637 *   {
638 *   #ifdef MMS
639 *   return -Utilities::fixed_power<2, double>(p[component]) -
641 *   #else
642 *   return (p[component] - strike_price) *
643 *   exp((-interest_rate) * (this->get_time()));
644 *   #endif
645 *   }
646 *  
647 *  
648 *  
649 * @endcode
650 *
651 * Finally, we handle the right hand side.
652 *
653 * @code
654 *   template <int dim>
655 *   class RightHandSide : public Function<dim>
656 *   {
657 *   public:
658 *   RightHandSide(const double asset_volatility, const double interest_rate);
659 *  
660 *   virtual double value(const Point<dim> &p,
661 *   const unsigned int component = 0) const override;
662 *  
663 *   private:
664 *   const double asset_volatility;
665 *   const double interest_rate;
666 *   };
667 *  
668 *  
669 *   template <int dim>
670 *   RightHandSide<dim>::RightHandSide(const double asset_volatility,
671 *   const double interest_rate)
672 *   : asset_volatility(asset_volatility)
673 *   , interest_rate(interest_rate)
674 *   {}
675 *  
676 *  
677 *   template <int dim>
678 *   double RightHandSide<dim>::value(const Point<dim> &p,
679 *   const unsigned int component) const
680 *   {
681 *   #ifdef MMS
682 *   return 2 * (this->get_time()) -
683 *   Utilities::fixed_power<2, double>(asset_volatility * p[component]) -
684 *   2 * interest_rate * Utilities::fixed_power<2, double>(p[component]) -
685 *   interest_rate *
686 *   (-Utilities::fixed_power<2, double>(p[component]) -
688 *   #else
689 *   (void)p;
690 *   (void)component;
691 *   return 0.0;
692 *   #endif
693 *   }
694 *  
695 *  
696 *  
697 * @endcode
698 *
699 *
700 * <a name="step_78-ThecodeBlackScholescodeClass"></a>
701 * <h3>The <code>BlackScholes</code> Class</h3>
702 *
703
704 *
705 * The next piece is the declaration of the main class of this program. This
706 * is very similar to the @ref step_26 "step-26" tutorial, with some modifications. New
707 * matrices had to be added to calculate the A and B matrices, as well as the
708 * @f$V_{diff}@f$ vector mentioned in the introduction. We also define the
709 * parameters used in the problem.
710 *
711
712 *
713 * - <code>maximum_stock_price</code>: The imposed upper bound on the spatial
714 * domain. This is the maximum allowed stock price.
715 * - <code>maturity_time</code>: The upper bound on the time domain. This is
716 * when the option expires.\n
717 * - <code>asset_volatility</code>: The volatility of the stock price.\n
718 * - <code>interest_rate</code>: The risk free interest rate.\n
719 * - <code>strike_price</code>: The agreed upon price that the buyer will
720 * have the option of purchasing the stocks at the expiration time.
721 *
722
723 *
724 * Some slight differences between this program and @ref step_26 "step-26" are the creation
725 * of the <code>a_matrix</code> and the <code>b_matrix</code>, which is
726 * described in the introduction. We then also need to store the current time,
727 * the size of the time step, and the number of the current time step.
728 * Next, we will store the output into a <code>DataOutStack</code>
729 * variable because we will be layering the solution at each time on top of
730 * one another to create the solution manifold. Then, we have a variable that
731 * stores the current cycle and number of cycles that we will run when
732 * calculating the solution. The cycle is one full solution calculation given
733 * a mesh. We refine the mesh once in between each cycle to exhibit the
734 * convergence properties of our program. Finally, we store the convergence
735 * data into a convergence table.
736 *
737
738 *
739 * As far as member functions are concerned, we have a function that
740 * calculates the convergence information for each cycle, called
741 * <code>process_solution</code>. This is just like what is done in @ref step_7 "step-7".
742 *
743 * @code
744 *   template <int dim>
745 *   class BlackScholes
746 *   {
747 *   public:
748 *   BlackScholes();
749 *  
750 *   void run();
751 *  
752 *   private:
753 *   void setup_system();
754 *   void solve_time_step();
755 *   void refine_grid();
756 *   void process_solution();
757 *   void add_results_for_output();
758 *   void write_convergence_table();
759 *  
760 *   const double maximum_stock_price;
761 *   const double maturity_time;
762 *   const double asset_volatility;
763 *   const double interest_rate;
764 *   const double strike_price;
765 *  
767 *   const FE_Q<dim> fe;
768 *   DoFHandler<dim> dof_handler;
769 *  
770 *   AffineConstraints<double> constraints;
771 *  
772 *   SparsityPattern sparsity_pattern;
774 *   SparseMatrix<double> laplace_matrix;
775 *   SparseMatrix<double> a_matrix;
776 *   SparseMatrix<double> b_matrix;
777 *   SparseMatrix<double> system_matrix;
778 *  
779 *   Vector<double> solution;
780 *   Vector<double> system_rhs;
781 *  
782 *   double time;
783 *   double time_step;
784 *  
785 *   const double theta;
786 *   const unsigned int n_cycles;
787 *   const unsigned int n_time_steps;
788 *  
789 *   DataOutStack<dim> data_out_stack;
790 *   std::vector<std::string> solution_names;
791 *  
792 *   ConvergenceTable convergence_table;
793 *   };
794 *  
795 * @endcode
796 *
797 *
798 * <a name="step_78-ThecodeBlackScholescodeImplementation"></a>
799 * <h3>The <code>BlackScholes</code> Implementation</h3>
800 *
801
802 *
803 * Now, we get to the implementation of the main class. We will set the values
804 * for the various parameters used in the problem. These were chosen because
805 * they are fairly normal values for these parameters. Although the stock
806 * price has no upper bound in reality (it is in fact infinite), we impose
807 * an upper bound that is twice the strike price. This is a somewhat arbitrary
808 * choice to be twice the strike price, but it is large enough to see the
809 * interesting parts of the solution.
810 *
811 * @code
812 *   template <int dim>
813 *   BlackScholes<dim>::BlackScholes()
814 *   : maximum_stock_price(1.)
815 *   , maturity_time(1.)
816 *   , asset_volatility(.2)
817 *   , interest_rate(0.05)
818 *   , strike_price(0.5)
819 *   , fe(1)
820 *   , dof_handler(triangulation)
821 *   , time(0.0)
822 *   , theta(0.5)
823 *   , n_cycles(4)
824 *   , n_time_steps(5000)
825 *   {
826 *   Assert(dim == 1, ExcNotImplemented());
827 *   }
828 *  
829 * @endcode
830 *
831 *
832 * <a name="step_78-codeBlackScholessetup_systemcode"></a>
833 * <h4><code>BlackScholes::setup_system</code></h4>
834 *
835
836 *
837 * The next function sets up the DoFHandler object, computes
838 * the constraints, and sets the linear algebra objects to their correct
839 * sizes. We also compute the @ref GlossMassMatrix "mass matrix" here by calling a function from the
840 * library. We will compute the other 3 matrices next, because these need to
841 * be computed 'by hand'.
842 *
843
844 *
845 * Note, that the time step is initialized here because the maturity time was
846 * needed to compute the time step.
847 *
848 * @code
849 *   template <int dim>
850 *   void BlackScholes<dim>::setup_system()
851 *   {
852 *   dof_handler.distribute_dofs(fe);
853 *  
854 *   time_step = maturity_time / n_time_steps;
855 *  
856 *   constraints.clear();
857 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
858 *   constraints.close();
859 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
860 *   DoFTools::make_sparsity_pattern(dof_handler,
861 *   dsp,
862 *   constraints,
863 *   /*keep_constrained_dofs = */ true);
864 *   sparsity_pattern.copy_from(dsp);
865 *  
866 *   mass_matrix.reinit(sparsity_pattern);
867 *   laplace_matrix.reinit(sparsity_pattern);
868 *   a_matrix.reinit(sparsity_pattern);
869 *   b_matrix.reinit(sparsity_pattern);
870 *   system_matrix.reinit(sparsity_pattern);
871 *  
872 *   MatrixCreator::create_mass_matrix(dof_handler,
873 *   QGauss<dim>(fe.degree + 1),
874 *   mass_matrix);
875 *  
876 * @endcode
877 *
878 * Below is the code to create the Laplace matrix with non-constant
879 * coefficients. This corresponds to the matrix D in the introduction. This
880 * non-constant coefficient is represented in the
881 * <code>current_coefficient</code> variable.
882 *
883 * @code
884 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
885 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
886 *   const QGauss<dim> quadrature_formula(fe.degree + 1);
887 *   FEValues<dim> fe_values(fe,
888 *   quadrature_formula,
891 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
892 *   for (const auto &cell : dof_handler.active_cell_iterators())
893 *   {
894 *   cell_matrix = 0.;
895 *   fe_values.reinit(cell);
896 *   for (const unsigned int q_index : fe_values.quadrature_point_indices())
897 *   {
898 *   const double current_coefficient =
899 *   fe_values.quadrature_point(q_index).square();
900 *   for (const unsigned int i : fe_values.dof_indices())
901 *   {
902 *   for (const unsigned int j : fe_values.dof_indices())
903 *   cell_matrix(i, j) +=
904 *   (current_coefficient * // (x_q)^2
905 *   fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
906 *   fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
907 *   fe_values.JxW(q_index)); // dx
908 *   }
909 *   }
910 *   cell->get_dof_indices(local_dof_indices);
911 *   for (const unsigned int i : fe_values.dof_indices())
912 *   {
913 *   for (const unsigned int j : fe_values.dof_indices())
914 *   laplace_matrix.add(local_dof_indices[i],
915 *   local_dof_indices[j],
916 *   cell_matrix(i, j));
917 *   }
918 *   }
919 *  
920 * @endcode
921 *
922 * Now we will create the A matrix. Below is the code to create the matrix A
923 * as discussed in the introduction. The non constant coefficient is again
924 * represented in the <code>current_coefficient</code> variable.
925 *
926 * @code
927 *   for (const auto &cell : dof_handler.active_cell_iterators())
928 *   {
929 *   cell_matrix = 0.;
930 *   fe_values.reinit(cell);
931 *   for (const unsigned int q_index : fe_values.quadrature_point_indices())
932 *   {
933 *   const Tensor<1, dim> current_coefficient =
934 *   fe_values.quadrature_point(q_index);
935 *   for (const unsigned int i : fe_values.dof_indices())
936 *   {
937 *   for (const unsigned int j : fe_values.dof_indices())
938 *   {
939 *   cell_matrix(i, j) +=
940 *   (current_coefficient * // x_q
941 *   fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
942 *   fe_values.shape_value(j, q_index) * // phi_j(x_q)
943 *   fe_values.JxW(q_index)); // dx
944 *   }
945 *   }
946 *   }
947 *   cell->get_dof_indices(local_dof_indices);
948 *   for (const unsigned int i : fe_values.dof_indices())
949 *   {
950 *   for (const unsigned int j : fe_values.dof_indices())
951 *   a_matrix.add(local_dof_indices[i],
952 *   local_dof_indices[j],
953 *   cell_matrix(i, j));
954 *   }
955 *   }
956 *  
957 * @endcode
958 *
959 * Finally we will create the matrix B. Below is the code to create the
960 * matrix B as discussed in the introduction. The non constant coefficient
961 * is again represented in the <code>current_coefficient</code> variable.
962 *
963 * @code
964 *   for (const auto &cell : dof_handler.active_cell_iterators())
965 *   {
966 *   cell_matrix = 0.;
967 *   fe_values.reinit(cell);
968 *   for (const unsigned int q_index : fe_values.quadrature_point_indices())
969 *   {
970 *   const Tensor<1, dim> current_coefficient =
971 *   fe_values.quadrature_point(q_index);
972 *   for (const unsigned int i : fe_values.dof_indices())
973 *   {
974 *   for (const unsigned int j : fe_values.dof_indices())
975 *   cell_matrix(i, j) +=
976 *   (current_coefficient * // x_q
977 *   fe_values.shape_value(i, q_index) * // phi_i(x_q)
978 *   fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
979 *   fe_values.JxW(q_index)); // dx
980 *   }
981 *   }
982 *   cell->get_dof_indices(local_dof_indices);
983 *   for (const unsigned int i : fe_values.dof_indices())
984 *   {
985 *   for (const unsigned int j : fe_values.dof_indices())
986 *   b_matrix.add(local_dof_indices[i],
987 *   local_dof_indices[j],
988 *   cell_matrix(i, j));
989 *   }
990 *   }
991 *  
992 *   solution.reinit(dof_handler.n_dofs());
993 *   system_rhs.reinit(dof_handler.n_dofs());
994 *   }
995 *  
996 * @endcode
997 *
998 *
999 * <a name="step_78-codeBlackScholessolve_time_stepcode"></a>
1000 * <h4><code>BlackScholes::solve_time_step</code></h4>
1001 *
1002
1003 *
1004 * The next function is the one that solves the actual linear system for a
1005 * single time step. The only interesting thing here is that the matrices
1006 * we have built are symmetric positive definite, so we can use the
1007 * conjugate gradient method.
1008 *
1009 * @code
1010 *   template <int dim>
1011 *   void BlackScholes<dim>::solve_time_step()
1012 *   {
1013 *   SolverControl solver_control(1000, 1e-12);
1014 *   SolverCG<Vector<double>> cg(solver_control);
1015 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
1016 *   preconditioner.initialize(system_matrix, 1.0);
1017 *   cg.solve(system_matrix, solution, system_rhs, preconditioner);
1018 *   constraints.distribute(solution);
1019 *   }
1020 *  
1021 * @endcode
1022 *
1023 *
1024 * <a name="step_78-codeBlackScholesadd_results_for_outputcode"></a>
1025 * <h4><code>BlackScholes::add_results_for_output</code></h4>
1026 *
1027
1028 *
1029 * This is simply the function to stitch the solution pieces together. For
1030 * this, we create a new layer at each time, and then add the solution vector
1031 * for that timestep. The function then stitches this together with the old
1032 * solutions using 'build_patches'.
1033 *
1034 * @code
1035 *   template <int dim>
1036 *   void BlackScholes<dim>::add_results_for_output()
1037 *   {
1038 *   data_out_stack.new_parameter_value(time, time_step);
1039 *   data_out_stack.attach_dof_handler(dof_handler);
1040 *   data_out_stack.add_data_vector(solution, solution_names);
1041 *   data_out_stack.build_patches(2);
1042 *   data_out_stack.finish_parameter_value();
1043 *   }
1044 *  
1045 * @endcode
1046 *
1047 *
1048 * <a name="step_78-codeBlackScholesrefine_gridcode"></a>
1049 * <h4><code>BlackScholes::refine_grid</code></h4>
1050 *
1051
1052 *
1053 * It is somewhat unnecessary to have a function for the global refinement
1054 * that we do. The reason for the function is to allow for the possibility of
1055 * an adaptive refinement later.
1056 *
1057 * @code
1058 *   template <int dim>
1059 *   void BlackScholes<dim>::refine_grid()
1060 *   {
1061 *   triangulation.refine_global(1);
1062 *   }
1063 *  
1064 * @endcode
1065 *
1066 *
1067 * <a name="step_78-codeBlackScholesprocess_solutioncode"></a>
1068 * <h4><code>BlackScholes::process_solution</code></h4>
1069 *
1070
1071 *
1072 * This is where we calculate the convergence and error data to evaluate the
1073 * effectiveness of the program. Here, we calculate the @f$L^2@f$, @f$H^1@f$ and
1074 * @f$L^{\infty}@f$ norms.
1075 *
1076 * @code
1077 *   template <int dim>
1078 *   void BlackScholes<dim>::process_solution()
1079 *   {
1080 *   Solution<dim> sol(maturity_time);
1081 *   sol.set_time(time);
1082 *   Vector<float> difference_per_cell(triangulation.n_active_cells());
1083 *   VectorTools::integrate_difference(dof_handler,
1084 *   solution,
1085 *   sol,
1086 *   difference_per_cell,
1087 *   QGauss<dim>(fe.degree + 1),
1088 *   VectorTools::L2_norm);
1089 *   const double L2_error =
1090 *   VectorTools::compute_global_error(triangulation,
1091 *   difference_per_cell,
1092 *   VectorTools::L2_norm);
1093 *   VectorTools::integrate_difference(dof_handler,
1094 *   solution,
1095 *   sol,
1096 *   difference_per_cell,
1097 *   QGauss<dim>(fe.degree + 1),
1098 *   VectorTools::H1_seminorm);
1099 *   const double H1_error =
1100 *   VectorTools::compute_global_error(triangulation,
1101 *   difference_per_cell,
1102 *   VectorTools::H1_seminorm);
1103 *   const QTrapezoid<1> q_trapezoid;
1104 *   const QIterated<dim> q_iterated(q_trapezoid, fe.degree * 2 + 1);
1105 *   VectorTools::integrate_difference(dof_handler,
1106 *   solution,
1107 *   sol,
1108 *   difference_per_cell,
1109 *   q_iterated,
1110 *   VectorTools::Linfty_norm);
1111 *   const double Linfty_error =
1112 *   VectorTools::compute_global_error(triangulation,
1113 *   difference_per_cell,
1114 *   VectorTools::Linfty_norm);
1115 *   const unsigned int n_active_cells = triangulation.n_active_cells();
1116 *   const unsigned int n_dofs = dof_handler.n_dofs();
1117 *   convergence_table.add_value("cells", n_active_cells);
1118 *   convergence_table.add_value("dofs", n_dofs);
1119 *   convergence_table.add_value("L2", L2_error);
1120 *   convergence_table.add_value("H1", H1_error);
1121 *   convergence_table.add_value("Linfty", Linfty_error);
1122 *   }
1123 *  
1124 * @endcode
1125 *
1126 *
1127 * <a name="step_78-codeBlackScholeswrite_convergence_tablecode"></a>
1128 * <h4><code>BlackScholes::write_convergence_table</code> </h4>
1129 *
1130
1131 *
1132 * This next part is building the convergence and error tables. By this, we
1133 * need to set the settings for how to output the data that was calculated
1134 * during <code>BlackScholes::process_solution</code>. First, we will create
1135 * the headings and set up the cells properly. During this, we will also
1136 * prescribe the precision of our results. Then we will write the calculated
1137 * errors based on the @f$L^2@f$, @f$H^1@f$, and @f$L^{\infty}@f$ norms to the console and
1138 * to the error LaTeX file.
1139 *
1140 * @code
1141 *   template <int dim>
1142 *   void BlackScholes<dim>::write_convergence_table()
1143 *   {
1144 *   convergence_table.set_precision("L2", 3);
1145 *   convergence_table.set_precision("H1", 3);
1146 *   convergence_table.set_precision("Linfty", 3);
1147 *   convergence_table.set_scientific("L2", true);
1148 *   convergence_table.set_scientific("H1", true);
1149 *   convergence_table.set_scientific("Linfty", true);
1150 *   convergence_table.set_tex_caption("cells", "\\# cells");
1151 *   convergence_table.set_tex_caption("dofs", "\\# dofs");
1152 *   convergence_table.set_tex_caption("L2", "@fL^2@f-error");
1153 *   convergence_table.set_tex_caption("H1", "@fH^1@f-error");
1154 *   convergence_table.set_tex_caption("Linfty", "@fL^\\infty@f-error");
1155 *   convergence_table.set_tex_format("cells", "r");
1156 *   convergence_table.set_tex_format("dofs", "r");
1157 *   std::cout << std::endl;
1158 *   convergence_table.write_text(std::cout);
1159 *   std::string error_filename = "error";
1160 *   error_filename += "-global";
1161 *   error_filename += ".tex";
1162 *   std::ofstream error_table_file(error_filename);
1163 *   convergence_table.write_tex(error_table_file);
1164 *  
1165 * @endcode
1166 *
1167 * Next, we will make the convergence table. We will again write this to
1168 * the console and to the convergence LaTeX file.
1169 *
1170 * @code
1171 *   convergence_table.add_column_to_supercolumn("cells", "n cells");
1172 *   std::vector<std::string> new_order;
1173 *   new_order.emplace_back("n cells");
1174 *   new_order.emplace_back("H1");
1175 *   new_order.emplace_back("L2");
1176 *   convergence_table.set_column_order(new_order);
1177 *   convergence_table.evaluate_convergence_rates(
1178 *   "L2", ConvergenceTable::reduction_rate);
1179 *   convergence_table.evaluate_convergence_rates(
1180 *   "L2", ConvergenceTable::reduction_rate_log2);
1181 *   convergence_table.evaluate_convergence_rates(
1182 *   "H1", ConvergenceTable::reduction_rate);
1183 *   convergence_table.evaluate_convergence_rates(
1184 *   "H1", ConvergenceTable::reduction_rate_log2);
1185 *   std::cout << std::endl;
1186 *   convergence_table.write_text(std::cout);
1187 *   std::string conv_filename = "convergence";
1188 *   conv_filename += "-global";
1189 *   switch (fe.degree)
1190 *   {
1191 *   case 1:
1192 *   conv_filename += "-q1";
1193 *   break;
1194 *   case 2:
1195 *   conv_filename += "-q2";
1196 *   break;
1197 *   default:
1198 *   DEAL_II_NOT_IMPLEMENTED();
1199 *   }
1200 *   conv_filename += ".tex";
1201 *   std::ofstream table_file(conv_filename);
1202 *   convergence_table.write_tex(table_file);
1203 *   }
1204 *  
1205 * @endcode
1206 *
1207 *
1208 * <a name="step_78-codeBlackScholesruncode"></a>
1209 * <h4><code>BlackScholes::run</code></h4>
1210 *
1211
1212 *
1213 * Now we get to the main driver of the program. This is where we do all the
1214 * work of looping through the time steps and calculating the solution vector
1215 * each time. Here at the top, we set the initial refinement value and then
1216 * create a mesh. Then we refine this mesh once. Next, we set up the
1217 * data_out_stack object to store our solution. Finally, we start a for loop
1218 * to loop through the cycles. This lets us recalculate a solution for each
1219 * successive mesh refinement. At the beginning of each iteration, we need to
1220 * reset the time and time step number. We introduce an if statement to
1221 * accomplish this because we don't want to do this on the first iteration.
1222 *
1223 * @code
1224 *   template <int dim>
1225 *   void BlackScholes<dim>::run()
1226 *   {
1227 *   GridGenerator::hyper_cube(triangulation, 0.0, maximum_stock_price, true);
1229 *  
1230 *   solution_names.emplace_back("u");
1231 *   data_out_stack.declare_data_vector(solution_names,
1233 *  
1234 *   Vector<double> vmult_result;
1235 *   Vector<double> forcing_terms;
1236 *  
1237 *   for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
1238 *   {
1239 *   if (cycle != 0)
1240 *   {
1241 *   refine_grid();
1242 *   time = 0.0;
1243 *   }
1244 *  
1245 *   setup_system();
1246 *  
1247 *   std::cout << std::endl
1248 *   << "===========================================" << std::endl
1249 *   << "Cycle " << cycle << ':' << std::endl
1250 *   << "Number of active cells: "
1251 *   << triangulation.n_active_cells() << std::endl
1252 *   << "Number of degrees of freedom: " << dof_handler.n_dofs()
1253 *   << std::endl
1254 *   << std::endl;
1255 *  
1256 *   VectorTools::interpolate(dof_handler,
1257 *   InitialConditions<dim>(strike_price),
1258 *   solution);
1259 *  
1260 *   if (cycle == (n_cycles - 1))
1261 *   {
1262 *   add_results_for_output();
1263 *   }
1264 *  
1265 * @endcode
1266 *
1267 * Next, we run the main loop which runs until we exceed the maturity
1268 * time. We first compute the right hand side of the equation, which is
1269 * described in the introduction. Recall that it contains the term
1270 * @f$\left[-\frac{1}{4}k_n\sigma^2\mathbf{D}-k_nr\mathbf{M}+k_n\sigma^2
1271 * \mathbf{B}-k_nr\mathbf{A}+\mathbf{M}\right]V^{n-1}@f$. We put these
1272 * terms into the variable system_rhs, with the help of a temporary
1273 * vector:
1274 *
1275 * @code
1276 *   vmult_result.reinit(dof_handler.n_dofs());
1277 *   forcing_terms.reinit(dof_handler.n_dofs());
1278 *   for (unsigned int timestep_number = 0; timestep_number < n_time_steps;
1279 *   ++timestep_number)
1280 *   {
1281 *   time += time_step;
1282 *  
1283 *   if (timestep_number % 1000 == 0)
1284 *   std::cout << "Time step " << timestep_number << " at t=" << time
1285 *   << std::endl;
1286 *  
1287 *   mass_matrix.vmult(system_rhs, solution);
1288 *  
1289 *   laplace_matrix.vmult(vmult_result, solution);
1290 *   system_rhs.add(
1291 *   (-1) * (1 - theta) * time_step *
1292 *   Utilities::fixed_power<2, double>(asset_volatility) * 0.5,
1293 *   vmult_result);
1294 *   mass_matrix.vmult(vmult_result, solution);
1295 *  
1296 *   system_rhs.add((-1) * (1 - theta) * time_step * interest_rate * 2,
1297 *   vmult_result);
1298 *  
1299 *   a_matrix.vmult(vmult_result, solution);
1300 *   system_rhs.add((-1) * time_step * interest_rate, vmult_result);
1301 *  
1302 *   b_matrix.vmult(vmult_result, solution);
1303 *   system_rhs.add(
1304 *   (-1) * Utilities::fixed_power<2, double>(asset_volatility) *
1305 *   time_step * 1,
1306 *   vmult_result);
1307 *  
1308 * @endcode
1309 *
1310 * The second piece is to compute the contributions of the source
1311 * terms. This corresponds to the term @f$-k_n\left[\frac{1}{2}F^{n-1}
1312 * +\frac{1}{2}F^n\right]@f$. The following code calls
1313 * VectorTools::create_right_hand_side to compute the vectors @f$F@f$,
1314 * where we set the time of the right hand side (source) function
1315 * before we evaluate it. The result of this all ends up in the
1316 * forcing_terms variable:
1317 *
1318 * @code
1319 *   RightHandSide<dim> rhs_function(asset_volatility, interest_rate);
1320 *   rhs_function.set_time(time);
1321 *   VectorTools::create_right_hand_side(dof_handler,
1322 *   QGauss<dim>(fe.degree + 1),
1323 *   rhs_function,
1324 *   forcing_terms);
1325 *   forcing_terms *= time_step * theta;
1326 *   system_rhs -= forcing_terms;
1327 *  
1328 *   rhs_function.set_time(time - time_step);
1329 *   VectorTools::create_right_hand_side(dof_handler,
1330 *   QGauss<dim>(fe.degree + 1),
1331 *   rhs_function,
1332 *   forcing_terms);
1333 *   forcing_terms *= time_step * (1 - theta);
1334 *   system_rhs -= forcing_terms;
1335 *  
1336 * @endcode
1337 *
1338 * Next, we add the forcing terms to the ones that come from the
1339 * time stepping, and also build the matrix @f$\left[\mathbf{M}+
1340 * \frac{1}{4}k_n\sigma^2\mathbf{D}+k_nr\mathbf{M}\right]@f$ that we
1341 * have to invert in each time step. The final piece of these
1342 * operations is to eliminate hanging node constrained degrees of
1343 * freedom from the linear system:
1344 *
1345 * @code
1346 *   system_matrix.copy_from(mass_matrix);
1347 *   system_matrix.add(
1348 *   (theta)*time_step *
1349 *   Utilities::fixed_power<2, double>(asset_volatility) * 0.5,
1350 *   laplace_matrix);
1351 *   system_matrix.add((time_step)*interest_rate * theta * (1 + 1),
1352 *   mass_matrix);
1353 *  
1354 *   constraints.condense(system_matrix, system_rhs);
1355 *  
1356 * @endcode
1357 *
1358 * There is one more operation we need to do before we can solve it:
1359 * boundary values. To this end, we create a boundary value object,
1360 * set the proper time to the one of the current time step, and
1361 * evaluate it as we have done many times before. The result is used
1362 * to also set the correct boundary values in the linear system:
1363 *
1364 * @code
1365 *   {
1366 *   RightBoundaryValues<dim> right_boundary_function(strike_price,
1367 *   interest_rate);
1368 *   LeftBoundaryValues<dim> left_boundary_function;
1369 *   right_boundary_function.set_time(time);
1370 *   left_boundary_function.set_time(time);
1371 *   std::map<types::global_dof_index, double> boundary_values;
1373 *   0,
1374 *   left_boundary_function,
1375 *   boundary_values);
1377 *   1,
1378 *   right_boundary_function,
1379 *   boundary_values);
1380 *   MatrixTools::apply_boundary_values(boundary_values,
1381 *   system_matrix,
1382 *   solution,
1383 *   system_rhs);
1384 *   }
1385 *  
1386 * @endcode
1387 *
1388 * With this out of the way, all we have to do is solve the system,
1389 * generate graphical data on the last cycle, and create the
1390 * convergence table data.
1391 *
1392 * @code
1393 *   solve_time_step();
1394 *  
1395 *   if (cycle == (n_cycles - 1))
1396 *   {
1397 *   add_results_for_output();
1398 *   }
1399 *   }
1400 *   #ifdef MMS
1401 *   process_solution();
1402 *   #endif
1403 *   }
1404 *  
1405 *   const std::string filename = "solution.vtk";
1406 *   std::ofstream output(filename);
1407 *   data_out_stack.write_vtk(output);
1408 *  
1409 *   #ifdef MMS
1410 *   write_convergence_table();
1411 *   #endif
1412 *   }
1413 *  
1414 *   } // namespace BlackScholesSolver
1415 *  
1416 * @endcode
1417 *
1418 *
1419 * <a name="step_78-ThecodemaincodeFunction"></a>
1420 * <h3>The <code>main</code> Function</h3>
1421 *
1422
1423 *
1424 * Having made it this far, there is, again, nothing much to discuss for the
1425 * main function of this program: it looks like all such functions since @ref step_6 "step-6".
1426 *
1427 * @code
1428 *   int main()
1429 *   {
1430 *   try
1431 *   {
1432 *   using namespace BlackScholesSolver;
1433 *  
1434 *   BlackScholes<1> black_scholes_solver;
1435 *   black_scholes_solver.run();
1436 *   }
1437 *   catch (std::exception &exc)
1438 *   {
1439 *   std::cerr << std::endl
1440 *   << std::endl
1441 *   << "----------------------------------------------------"
1442 *   << std::endl;
1443 *   std::cerr << "Exception on processing: " << std::endl
1444 *   << exc.what() << std::endl
1445 *   << "Aborting!" << std::endl
1446 *   << "----------------------------------------------------"
1447 *   << std::endl;
1448 *   return 1;
1449 *   }
1450 *   catch (...)
1451 *   {
1452 *   std::cerr << std::endl
1453 *   << std::endl
1454 *   << "----------------------------------------------------"
1455 *   << std::endl;
1456 *   std::cerr << "Unknown exception!" << std::endl
1457 *   << "Aborting!" << std::endl
1458 *   << "----------------------------------------------------"
1459 *   << std::endl;
1460 *   return 1;
1461 *   }
1462 *   return 0;
1463 *   }
1464 * @endcode
1465<a name="step_78-Results"></a><h1>Results</h1>
1466
1467
1468
1469Below is the output of the program:
1470@code
1471===========================================
1472Number of active cells: 1
1473Number of degrees of freedom: 2
1474
1475Time step 0 at t=0.0002
1476[...]
1477
1478Cycle 7:
1479Number of active cells: 128
1480Number of degrees of freedom: 129
1481
1482Time step 0 at t=0.0002
1483Time step 1000 at t=0.2002
1484Time step 2000 at t=0.4002
1485Time step 3000 at t=0.6002
1486Time step 4000 at t=0.8002
1487
1488cells dofs L2 H1 Linfty
1489 1 2 1.667e-01 5.774e-01 2.222e-01
1490 2 3 3.906e-02 2.889e-01 5.380e-02
1491 4 5 9.679e-03 1.444e-01 1.357e-02
1492 8 9 2.405e-03 7.218e-02 3.419e-03
1493 16 17 5.967e-04 3.609e-02 8.597e-04
1494 32 33 1.457e-04 1.804e-02 2.155e-04
1495 64 65 3.307e-05 9.022e-03 5.388e-05
1496 128 129 5.016e-06 4.511e-03 1.342e-05
1497
1498n cells H1 L2
1499 1 5.774e-01 - - 1.667e-01 - -
1500 2 2.889e-01 2.00 1.00 3.906e-02 4.27 2.09
1501 4 1.444e-01 2.00 1.00 9.679e-03 4.04 2.01
1502 8 7.218e-02 2.00 1.00 2.405e-03 4.02 2.01
1503 16 3.609e-02 2.00 1.00 5.967e-04 4.03 2.01
1504 32 1.804e-02 2.00 1.00 1.457e-04 4.10 2.03
1505 64 9.022e-03 2.00 1.00 3.307e-05 4.41 2.14
1506 128 4.511e-03 2.00 1.00 5.016e-06 6.59 2.72
1507@endcode
1508
1509What is more interesting is the output of the convergence tables. They are
1510outputted into the console, as well into a LaTeX file. The convergence tables
1511are shown above. Here, you can see that the solution has a convergence rate
1512of @f$\mathcal{O}(h)@f$ with respect to the @f$H^1@f$-norm, and the solution has a convergence rate
1513of @f$\mathcal{O}(h^2)@f$ with respect to the @f$L^2@f$-norm.
1514
1515
1516Below is the visualization of the solution.
1517
1518<div style="text-align:center;">
1519 <img src="https://www.dealii.org/images/steps/developer/step-78.mms-solution.png"
1520 alt="Solution of the MMS problem.">
1521</div>
1522 *
1523 *
1524<a name="step_78-PlainProg"></a>
1525<h1> The plain program</h1>
1526@include "step-78.cc"
1527*/
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
Definition fe_q.h:554
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Definition point.h:111
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
__global__ void set(Number *val, const Number s, const size_type N)
#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())
Definition loop.h:564
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.)
Definition advection.h:74
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition l2.h:159
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition l2.h:57
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 >())
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * end(VectorType &V)
void free(T *&pointer)
Definition cuda.h:96
std::string get_time()
constexpr T fixed_power(const T t)
Definition utilities.h:942
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={})
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})
void create_right_hand_side(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, const Function< spacedim, typename VectorType::value_type > &rhs, VectorType &rhs_vector, const AffineConstraints< typename VectorType::value_type > &constraints=AffineConstraints< typename VectorType::value_type >())
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 > &)