deal.II version GIT relicensing-1721-g8100761196 2024-08-31 12:30:00+00:00
\(\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) const
597 *   {
598 *   #ifdef MMS
599 *   return -Utilities::fixed_power<2, double>(this->get_time()) + 6;
600 *   #else
601 *   return 0.;
602 *   #endif
603 *   }
604 *  
605 *  
606 *  
607 * @endcode
608 *
609 * Then, we handle the right boundary condition.
610 *
611 * @code
612 *   template <int dim>
613 *   class RightBoundaryValues : public Function<dim>
614 *   {
615 *   public:
616 *   RightBoundaryValues(const double strike_price, const double interest_rate);
617 *  
618 *   virtual double value(const Point<dim> &p,
619 *   const unsigned int component = 0) const override;
620 *  
621 *   private:
622 *   const double strike_price;
623 *   const double interest_rate;
624 *   };
625 *  
626 *  
627 *   template <int dim>
628 *   RightBoundaryValues<dim>::RightBoundaryValues(const double strike_price,
629 *   const double interest_rate)
630 *   : strike_price(strike_price)
631 *   , interest_rate(interest_rate)
632 *   {}
633 *  
634 *  
635 *   template <int dim>
636 *   double RightBoundaryValues<dim>::value(const Point<dim> &p,
637 *   const unsigned int component) const
638 *   {
639 *   #ifdef MMS
640 *   return -Utilities::fixed_power<2, double>(p[component]) -
641 *   Utilities::fixed_power<2, double>(this->get_time()) + 6;
642 *   #else
643 *   return (p[component] - strike_price) *
644 *   exp((-interest_rate) * (this->get_time()));
645 *   #endif
646 *   }
647 *  
648 *  
649 *  
650 * @endcode
651 *
652 * Finally, we handle the right hand side.
653 *
654 * @code
655 *   template <int dim>
656 *   class RightHandSide : public Function<dim>
657 *   {
658 *   public:
659 *   RightHandSide(const double asset_volatility, const double interest_rate);
660 *  
661 *   virtual double value(const Point<dim> &p,
662 *   const unsigned int component = 0) const override;
663 *  
664 *   private:
665 *   const double asset_volatility;
666 *   const double interest_rate;
667 *   };
668 *  
669 *  
670 *   template <int dim>
671 *   RightHandSide<dim>::RightHandSide(const double asset_volatility,
672 *   const double interest_rate)
673 *   : asset_volatility(asset_volatility)
674 *   , interest_rate(interest_rate)
675 *   {}
676 *  
677 *  
678 *   template <int dim>
679 *   double RightHandSide<dim>::value(const Point<dim> &p,
680 *   const unsigned int component) const
681 *   {
682 *   #ifdef MMS
683 *   return 2 * (this->get_time()) -
684 *   Utilities::fixed_power<2, double>(asset_volatility * p[component]) -
685 *   2 * interest_rate * Utilities::fixed_power<2, double>(p[component]) -
686 *   interest_rate *
687 *   (-Utilities::fixed_power<2, double>(p[component]) -
688 *   Utilities::fixed_power<2, double>(this->get_time()) + 6);
689 *   #else
690 *   (void)p;
691 *   (void)component;
692 *   return 0.0;
693 *   #endif
694 *   }
695 *  
696 *  
697 *  
698 * @endcode
699 *
700 *
701 * <a name="step_78-ThecodeBlackScholescodeClass"></a>
702 * <h3>The <code>BlackScholes</code> Class</h3>
703 *
704
705 *
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.
711 *
712
713 *
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.
722 *
723
724 *
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.
737 *
738
739 *
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".
743 *
744 * @code
745 *   template <int dim>
746 *   class BlackScholes
747 *   {
748 *   public:
749 *   BlackScholes();
750 *  
751 *   void run();
752 *  
753 *   private:
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();
760 *  
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;
766 *  
768 *   const FE_Q<dim> fe;
769 *   DoFHandler<dim> dof_handler;
770 *  
771 *   AffineConstraints<double> constraints;
772 *  
773 *   SparsityPattern sparsity_pattern;
775 *   SparseMatrix<double> laplace_matrix;
776 *   SparseMatrix<double> a_matrix;
777 *   SparseMatrix<double> b_matrix;
778 *   SparseMatrix<double> system_matrix;
779 *  
780 *   Vector<double> solution;
781 *   Vector<double> system_rhs;
782 *  
783 *   double time;
784 *   double time_step;
785 *  
786 *   const double theta;
787 *   const unsigned int n_cycles;
788 *   const unsigned int n_time_steps;
789 *  
790 *   DataOutStack<dim> data_out_stack;
791 *   std::vector<std::string> solution_names;
792 *  
793 *   ConvergenceTable convergence_table;
794 *   };
795 *  
796 * @endcode
797 *
798 *
799 * <a name="step_78-ThecodeBlackScholescodeImplementation"></a>
800 * <h3>The <code>BlackScholes</code> Implementation</h3>
801 *
802
803 *
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.
811 *
812 * @code
813 *   template <int dim>
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)
820 *   , fe(1)
821 *   , dof_handler(triangulation)
822 *   , time(0.0)
823 *   , theta(0.5)
824 *   , n_cycles(4)
825 *   , n_time_steps(5000)
826 *   {
827 *   Assert(dim == 1, ExcNotImplemented());
828 *   }
829 *  
830 * @endcode
831 *
832 *
833 * <a name="step_78-codeBlackScholessetup_systemcode"></a>
834 * <h4><code>BlackScholes::setup_system</code></h4>
835 *
836
837 *
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'.
843 *
844
845 *
846 * Note, that the time step is initialized here because the maturity time was
847 * needed to compute the time step.
848 *
849 * @code
850 *   template <int dim>
851 *   void BlackScholes<dim>::setup_system()
852 *   {
853 *   dof_handler.distribute_dofs(fe);
854 *  
855 *   time_step = maturity_time / n_time_steps;
856 *  
857 *   constraints.clear();
858 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
859 *   constraints.close();
860 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
861 *   DoFTools::make_sparsity_pattern(dof_handler,
862 *   dsp,
863 *   constraints,
864 *   /*keep_constrained_dofs = */ true);
865 *   sparsity_pattern.copy_from(dsp);
866 *  
867 *   mass_matrix.reinit(sparsity_pattern);
868 *   laplace_matrix.reinit(sparsity_pattern);
869 *   a_matrix.reinit(sparsity_pattern);
870 *   b_matrix.reinit(sparsity_pattern);
871 *   system_matrix.reinit(sparsity_pattern);
872 *  
873 *   MatrixCreator::create_mass_matrix(dof_handler,
874 *   QGauss<dim>(fe.degree + 1),
875 *   mass_matrix);
876 *  
877 * @endcode
878 *
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.
883 *
884 * @code
885 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
886 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
887 *   const QGauss<dim> quadrature_formula(fe.degree + 1);
888 *   FEValues<dim> fe_values(fe,
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())
894 *   {
895 *   cell_matrix = 0.;
896 *   fe_values.reinit(cell);
897 *   for (const unsigned int q_index : fe_values.quadrature_point_indices())
898 *   {
899 *   const double current_coefficient =
900 *   fe_values.quadrature_point(q_index).square();
901 *   for (const unsigned int i : fe_values.dof_indices())
902 *   {
903 *   for (const unsigned int j : fe_values.dof_indices())
904 *   cell_matrix(i, j) +=
905 *   (current_coefficient * // (x_q)^2
906 *   fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
907 *   fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
908 *   fe_values.JxW(q_index)); // dx
909 *   }
910 *   }
911 *   cell->get_dof_indices(local_dof_indices);
912 *   for (const unsigned int i : fe_values.dof_indices())
913 *   {
914 *   for (const unsigned int j : fe_values.dof_indices())
915 *   laplace_matrix.add(local_dof_indices[i],
916 *   local_dof_indices[j],
917 *   cell_matrix(i, j));
918 *   }
919 *   }
920 *  
921 * @endcode
922 *
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.
926 *
927 * @code
928 *   for (const auto &cell : dof_handler.active_cell_iterators())
929 *   {
930 *   cell_matrix = 0.;
931 *   fe_values.reinit(cell);
932 *   for (const unsigned int q_index : fe_values.quadrature_point_indices())
933 *   {
934 *   const Tensor<1, dim> current_coefficient =
935 *   fe_values.quadrature_point(q_index);
936 *   for (const unsigned int i : fe_values.dof_indices())
937 *   {
938 *   for (const unsigned int j : fe_values.dof_indices())
939 *   {
940 *   cell_matrix(i, j) +=
941 *   (current_coefficient * // x_q
942 *   fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
943 *   fe_values.shape_value(j, q_index) * // phi_j(x_q)
944 *   fe_values.JxW(q_index)); // dx
945 *   }
946 *   }
947 *   }
948 *   cell->get_dof_indices(local_dof_indices);
949 *   for (const unsigned int i : fe_values.dof_indices())
950 *   {
951 *   for (const unsigned int j : fe_values.dof_indices())
952 *   a_matrix.add(local_dof_indices[i],
953 *   local_dof_indices[j],
954 *   cell_matrix(i, j));
955 *   }
956 *   }
957 *  
958 * @endcode
959 *
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.
963 *
964 * @code
965 *   for (const auto &cell : dof_handler.active_cell_iterators())
966 *   {
967 *   cell_matrix = 0.;
968 *   fe_values.reinit(cell);
969 *   for (const unsigned int q_index : fe_values.quadrature_point_indices())
970 *   {
971 *   const Tensor<1, dim> current_coefficient =
972 *   fe_values.quadrature_point(q_index);
973 *   for (const unsigned int i : fe_values.dof_indices())
974 *   {
975 *   for (const unsigned int j : fe_values.dof_indices())
976 *   cell_matrix(i, j) +=
977 *   (current_coefficient * // x_q
978 *   fe_values.shape_value(i, q_index) * // phi_i(x_q)
979 *   fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
980 *   fe_values.JxW(q_index)); // dx
981 *   }
982 *   }
983 *   cell->get_dof_indices(local_dof_indices);
984 *   for (const unsigned int i : fe_values.dof_indices())
985 *   {
986 *   for (const unsigned int j : fe_values.dof_indices())
987 *   b_matrix.add(local_dof_indices[i],
988 *   local_dof_indices[j],
989 *   cell_matrix(i, j));
990 *   }
991 *   }
992 *  
993 *   solution.reinit(dof_handler.n_dofs());
994 *   system_rhs.reinit(dof_handler.n_dofs());
995 *   }
996 *  
997 * @endcode
998 *
999 *
1000 * <a name="step_78-codeBlackScholessolve_time_stepcode"></a>
1001 * <h4><code>BlackScholes::solve_time_step</code></h4>
1002 *
1003
1004 *
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
1007 * we have built are symmetric positive definite, so we can use the
1008 * conjugate gradient method.
1009 *
1010 * @code
1011 *   template <int dim>
1012 *   void BlackScholes<dim>::solve_time_step()
1013 *   {
1014 *   SolverControl solver_control(1000, 1e-12);
1015 *   SolverCG<Vector<double>> cg(solver_control);
1016 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
1017 *   preconditioner.initialize(system_matrix, 1.0);
1018 *   cg.solve(system_matrix, solution, system_rhs, preconditioner);
1019 *   constraints.distribute(solution);
1020 *   }
1021 *  
1022 * @endcode
1023 *
1024 *
1025 * <a name="step_78-codeBlackScholesadd_results_for_outputcode"></a>
1026 * <h4><code>BlackScholes::add_results_for_output</code></h4>
1027 *
1028
1029 *
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'.
1034 *
1035 * @code
1036 *   template <int dim>
1037 *   void BlackScholes<dim>::add_results_for_output()
1038 *   {
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();
1044 *   }
1045 *  
1046 * @endcode
1047 *
1048 *
1049 * <a name="step_78-codeBlackScholesrefine_gridcode"></a>
1050 * <h4><code>BlackScholes::refine_grid</code></h4>
1051 *
1052
1053 *
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.
1057 *
1058 * @code
1059 *   template <int dim>
1060 *   void BlackScholes<dim>::refine_grid()
1061 *   {
1062 *   triangulation.refine_global(1);
1063 *   }
1064 *  
1065 * @endcode
1066 *
1067 *
1068 * <a name="step_78-codeBlackScholesprocess_solutioncode"></a>
1069 * <h4><code>BlackScholes::process_solution</code></h4>
1070 *
1071
1072 *
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.
1076 *
1077 * @code
1078 *   template <int dim>
1079 *   void BlackScholes<dim>::process_solution()
1080 *   {
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,
1085 *   solution,
1086 *   sol,
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,
1095 *   solution,
1096 *   sol,
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,
1107 *   solution,
1108 *   sol,
1109 *   difference_per_cell,
1110 *   q_iterated,
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);
1123 *   }
1124 *  
1125 * @endcode
1126 *
1127 *
1128 * <a name="step_78-codeBlackScholeswrite_convergence_tablecode"></a>
1129 * <h4><code>BlackScholes::write_convergence_table</code> </h4>
1130 *
1131
1132 *
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.
1140 *
1141 * @code
1142 *   template <int dim>
1143 *   void BlackScholes<dim>::write_convergence_table()
1144 *   {
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);
1165 *  
1166 * @endcode
1167 *
1168 * Next, we will make the convergence table. We will again write this to
1169 * the console and to the convergence LaTeX file.
1170 *
1171 * @code
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)
1191 *   {
1192 *   case 1:
1193 *   conv_filename += "-q1";
1194 *   break;
1195 *   case 2:
1196 *   conv_filename += "-q2";
1197 *   break;
1198 *   default:
1199 *   DEAL_II_NOT_IMPLEMENTED();
1200 *   }
1201 *   conv_filename += ".tex";
1202 *   std::ofstream table_file(conv_filename);
1203 *   convergence_table.write_tex(table_file);
1204 *   }
1205 *  
1206 * @endcode
1207 *
1208 *
1209 * <a name="step_78-codeBlackScholesruncode"></a>
1210 * <h4><code>BlackScholes::run</code></h4>
1211 *
1212
1213 *
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.
1223 *
1224 * @code
1225 *   template <int dim>
1226 *   void BlackScholes<dim>::run()
1227 *   {
1228 *   GridGenerator::hyper_cube(triangulation, 0.0, maximum_stock_price, true);
1229 *   triangulation.refine_global(0);
1230 *  
1231 *   solution_names.emplace_back("u");
1232 *   data_out_stack.declare_data_vector(solution_names,
1234 *  
1235 *   Vector<double> vmult_result;
1236 *   Vector<double> forcing_terms;
1237 *  
1238 *   for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
1239 *   {
1240 *   if (cycle != 0)
1241 *   {
1242 *   refine_grid();
1243 *   time = 0.0;
1244 *   }
1245 *  
1246 *   setup_system();
1247 *  
1248 *   std::cout << std::endl
1249 *   << "===========================================" << std::endl
1250 *   << "Cycle " << cycle << ':' << std::endl
1251 *   << "Number of active cells: "
1252 *   << triangulation.n_active_cells() << std::endl
1253 *   << "Number of degrees of freedom: " << dof_handler.n_dofs()
1254 *   << std::endl
1255 *   << std::endl;
1256 *  
1257 *   VectorTools::interpolate(dof_handler,
1258 *   InitialConditions<dim>(strike_price),
1259 *   solution);
1260 *  
1261 *   if (cycle == (n_cycles - 1))
1262 *   {
1263 *   add_results_for_output();
1264 *   }
1265 *  
1266 * @endcode
1267 *
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
1274 * vector:
1275 *
1276 * @code
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)
1281 *   {
1282 *   time += time_step;
1283 *  
1284 *   if (timestep_number % 1000 == 0)
1285 *   std::cout << "Time step " << timestep_number << " at t=" << time
1286 *   << std::endl;
1287 *  
1288 *   mass_matrix.vmult(system_rhs, solution);
1289 *  
1290 *   laplace_matrix.vmult(vmult_result, solution);
1291 *   system_rhs.add(
1292 *   (-1) * (1 - theta) * time_step *
1293 *   Utilities::fixed_power<2, double>(asset_volatility) * 0.5,
1294 *   vmult_result);
1295 *   mass_matrix.vmult(vmult_result, solution);
1296 *  
1297 *   system_rhs.add((-1) * (1 - theta) * time_step * interest_rate * 2,
1298 *   vmult_result);
1299 *  
1300 *   a_matrix.vmult(vmult_result, solution);
1301 *   system_rhs.add((-1) * time_step * interest_rate, vmult_result);
1302 *  
1303 *   b_matrix.vmult(vmult_result, solution);
1304 *   system_rhs.add(
1305 *   (-1) * Utilities::fixed_power<2, double>(asset_volatility) *
1306 *   time_step * 1,
1307 *   vmult_result);
1308 *  
1309 * @endcode
1310 *
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
1314 * VectorTools::create_right_hand_side to compute the vectors @f$F@f$,
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:
1318 *
1319 * @code
1320 *   RightHandSide<dim> rhs_function(asset_volatility, interest_rate);
1321 *   rhs_function.set_time(time);
1322 *   VectorTools::create_right_hand_side(dof_handler,
1323 *   QGauss<dim>(fe.degree + 1),
1324 *   rhs_function,
1325 *   forcing_terms);
1326 *   forcing_terms *= time_step * theta;
1327 *   system_rhs -= forcing_terms;
1328 *  
1329 *   rhs_function.set_time(time - time_step);
1330 *   VectorTools::create_right_hand_side(dof_handler,
1331 *   QGauss<dim>(fe.degree + 1),
1332 *   rhs_function,
1333 *   forcing_terms);
1334 *   forcing_terms *= time_step * (1 - theta);
1335 *   system_rhs -= forcing_terms;
1336 *  
1337 * @endcode
1338 *
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:
1345 *
1346 * @code
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,
1351 *   laplace_matrix);
1352 *   system_matrix.add((time_step)*interest_rate * theta * (1 + 1),
1353 *   mass_matrix);
1354 *  
1355 *   constraints.condense(system_matrix, system_rhs);
1356 *  
1357 * @endcode
1358 *
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:
1364 *
1365 * @code
1366 *   {
1367 *   RightBoundaryValues<dim> right_boundary_function(strike_price,
1368 *   interest_rate);
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;
1374 *   0,
1375 *   left_boundary_function,
1376 *   boundary_values);
1378 *   1,
1379 *   right_boundary_function,
1380 *   boundary_values);
1381 *   MatrixTools::apply_boundary_values(boundary_values,
1382 *   system_matrix,
1383 *   solution,
1384 *   system_rhs);
1385 *   }
1386 *  
1387 * @endcode
1388 *
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.
1392 *
1393 * @code
1394 *   solve_time_step();
1395 *  
1396 *   if (cycle == (n_cycles - 1))
1397 *   {
1398 *   add_results_for_output();
1399 *   }
1400 *   }
1401 *   #ifdef MMS
1402 *   process_solution();
1403 *   #endif
1404 *   }
1405 *  
1406 *   const std::string filename = "solution.vtk";
1407 *   std::ofstream output(filename);
1408 *   data_out_stack.write_vtk(output);
1409 *  
1410 *   #ifdef MMS
1411 *   write_convergence_table();
1412 *   #endif
1413 *   }
1414 *  
1415 *   } // namespace BlackScholesSolver
1416 *  
1417 * @endcode
1418 *
1419 *
1420 * <a name="step_78-ThecodemaincodeFunction"></a>
1421 * <h3>The <code>main</code> Function</h3>
1422 *
1423
1424 *
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".
1427 *
1428 * @code
1429 *   int main()
1430 *   {
1431 *   try
1432 *   {
1433 *   using namespace BlackScholesSolver;
1434 *  
1435 *   BlackScholes<1> black_scholes_solver;
1436 *   black_scholes_solver.run();
1437 *   }
1438 *   catch (std::exception &exc)
1439 *   {
1440 *   std::cerr << std::endl
1441 *   << std::endl
1442 *   << "----------------------------------------------------"
1443 *   << std::endl;
1444 *   std::cerr << "Exception on processing: " << std::endl
1445 *   << exc.what() << std::endl
1446 *   << "Aborting!" << std::endl
1447 *   << "----------------------------------------------------"
1448 *   << std::endl;
1449 *   return 1;
1450 *   }
1451 *   catch (...)
1452 *   {
1453 *   std::cerr << std::endl
1454 *   << std::endl
1455 *   << "----------------------------------------------------"
1456 *   << std::endl;
1457 *   std::cerr << "Unknown exception!" << std::endl
1458 *   << "Aborting!" << std::endl
1459 *   << "----------------------------------------------------"
1460 *   << std::endl;
1461 *   return 1;
1462 *   }
1463 *   return 0;
1464 *   }
1465 * @endcode
1466<a name="step_78-Results"></a><h1>Results</h1>
1467
1468
1469
1470Below is the output of the program:
1471@code
1472===========================================
1473Number of active cells: 1
1474Number of degrees of freedom: 2
1475
1476Time step 0 at t=0.0002
1477[...]
1478
1479Cycle 7:
1480Number of active cells: 128
1481Number of degrees of freedom: 129
1482
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
1488
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
1498
1499n cells H1 L2
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
1508@endcode
1509
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.
1515
1516
1517Below is the visualization of the solution.
1518
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.">
1522</div>
1523 *
1524 *
1525<a name="step_78-PlainProg"></a>
1526<h1> The plain program</h1>
1527@include "step-78.cc"
1528*/
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())
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
#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)
std::string get_time()
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={}, const unsigned int level=numbers::invalid_unsigned_int)
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 > &)