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-23.h
Go to the documentation of this file.
1,
661 *   const unsigned int component = 0) const override
662 *   {
663 *   (void)component;
664 *   Assert(component == 0, ExcIndexRange(component, 0, 1));
665 *   return 0;
666 *   }
667 *   };
668 *  
669 *  
670 *  
671 *   template <int dim>
672 *   class InitialValuesV : public Function<dim>
673 *   {
674 *   public:
675 *   virtual double value(const Point<dim> & /*p*/,
676 *   const unsigned int component = 0) const override
677 *   {
678 *   (void)component;
679 *   Assert(component == 0, ExcIndexRange(component, 0, 1));
680 *   return 0;
681 *   }
682 *   };
683 *  
684 *  
685 *  
686 * @endcode
687 *
688 * Secondly, we have the right hand side forcing term. Boring as we are, we
689 * choose zero here as well:
690 *
691 * @code
692 *   template <int dim>
693 *   class RightHandSide : public Function<dim>
694 *   {
695 *   public:
696 *   virtual double value(const Point<dim> & /*p*/,
697 *   const unsigned int component = 0) const override
698 *   {
699 *   (void)component;
700 *   Assert(component == 0, ExcIndexRange(component, 0, 1));
701 *   return 0;
702 *   }
703 *   };
704 *  
705 *  
706 *  
707 * @endcode
708 *
709 * Finally, we have boundary values for @f$u@f$ and @f$v@f$. They are as described
710 * in the introduction, one being the time derivative of the other:
711 *
712 * @code
713 *   template <int dim>
714 *   class BoundaryValuesU : public Function<dim>
715 *   {
716 *   public:
717 *   virtual double value(const Point<dim> &p,
718 *   const unsigned int component = 0) const override
719 *   {
720 *   (void)component;
721 *   Assert(component == 0, ExcIndexRange(component, 0, 1));
722 *  
723 *   if ((this->get_time() <= 0.5) && (p[0] < 0) && (p[1] < 1. / 3) &&
724 *   (p[1] > -1. / 3))
725 *   return std::sin(this->get_time() * 4 * numbers::PI);
726 *   else
727 *   return 0;
728 *   }
729 *   };
730 *  
731 *  
732 *  
733 *   template <int dim>
734 *   class BoundaryValuesV : public Function<dim>
735 *   {
736 *   public:
737 *   virtual double value(const Point<dim> &p,
738 *   const unsigned int component = 0) const override
739 *   {
740 *   (void)component;
741 *   Assert(component == 0, ExcIndexRange(component, 0, 1));
742 *  
743 *   if ((this->get_time() <= 0.5) && (p[0] < 0) && (p[1] < 1. / 3) &&
744 *   (p[1] > -1. / 3))
745 *   return (std::cos(this->get_time() * 4 * numbers::PI) * 4 * numbers::PI);
746 *   else
747 *   return 0;
748 *   }
749 *   };
750 *  
751 *  
752 *  
753 * @endcode
754 *
755 *
756 * <a name="step_23-ImplementationofthecodeWaveEquationcodeclass"></a>
757 * <h3>Implementation of the <code>WaveEquation</code> class</h3>
758 *
759
760 *
761 * The implementation of the actual logic is actually fairly short, since we
762 * relegate things like assembling the matrices and right hand side vectors
763 * to the library. The rest boils down to not much more than 130 lines of
764 * actual code, a significant fraction of which is boilerplate code that can
765 * be taken from previous example programs (e.g. the functions that solve
766 * linear systems, or that generate output).
767 *
768
769 *
770 * Let's start with the constructor (for an explanation of the choice of
771 * time step, see the section on Courant, Friedrichs, and Lewy in the
772 * introduction):
773 *
774 * @code
775 *   template <int dim>
776 *   WaveEquation<dim>::WaveEquation()
777 *   : fe(1)
778 *   , dof_handler(triangulation)
779 *   , time_step(1. / 64)
780 *   , time(time_step)
781 *   , timestep_number(1)
782 *   , theta(0.5)
783 *   {}
784 *  
785 *  
786 * @endcode
787 *
788 *
789 * <a name="step_23-WaveEquationsetup_system"></a>
790 * <h4>WaveEquation::setup_system</h4>
791 *
792
793 *
794 * The next function is the one that sets up the mesh, DoFHandler, and
795 * matrices and vectors at the beginning of the program, i.e. before the
796 * first time step. The first few lines are pretty much standard if you've
797 * read through the tutorial programs at least up to @ref step_6 "step-6":
798 *
799 * @code
800 *   template <int dim>
801 *   void WaveEquation<dim>::setup_system()
802 *   {
804 *   triangulation.refine_global(7);
805 *  
806 *   std::cout << "Number of active cells: " << triangulation.n_active_cells()
807 *   << std::endl;
808 *  
809 *   dof_handler.distribute_dofs(fe);
810 *  
811 *   std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
812 *   << std::endl
813 *   << std::endl;
814 *  
815 *   DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
816 *   DoFTools::make_sparsity_pattern(dof_handler, dsp);
817 *   sparsity_pattern.copy_from(dsp);
818 *  
819 * @endcode
820 *
821 * Then comes a block where we have to initialize the 3 matrices we need
822 * in the course of the program: the mass matrix, the Laplace matrix, and
823 * the matrix @f$M+k^2\theta^2A@f$ used when solving for @f$U^n@f$ in each time
824 * step.
825 *
826
827 *
828 * When setting up these matrices, note that they all make use of the same
829 * sparsity pattern object. Finally, the reason why matrices and sparsity
830 * patterns are separate objects in deal.II (unlike in many other finite
831 * element or linear algebra classes) becomes clear: in a significant
832 * fraction of applications, one has to hold several matrices that happen
833 * to have the same sparsity pattern, and there is no reason for them not
834 * to share this information, rather than re-building and wasting memory
835 * on it several times.
836 *
837
838 *
839 * After initializing all of these matrices, we call library functions
840 * that build the Laplace and mass matrices. All they need is a DoFHandler
841 * object and a quadrature formula object that is to be used for numerical
842 * integration. Note that in many respects these functions are better than
843 * what we would usually do in application programs, for example because
844 * they automatically parallelize building the matrices if multiple
845 * processors are available in a machine: for more information see the
846 * documentation of WorkStream or the
847 * @ref threads "Parallel computing with multiple processors"
848 * topic. The matrices for solving linear systems will be filled in the
849 * run() method because we need to re-apply boundary conditions every time
850 * step.
851 *
852 * @code
853 *   mass_matrix.reinit(sparsity_pattern);
854 *   laplace_matrix.reinit(sparsity_pattern);
855 *   matrix_u.reinit(sparsity_pattern);
856 *   matrix_v.reinit(sparsity_pattern);
857 *  
858 *   MatrixCreator::create_mass_matrix(dof_handler,
859 *   QGauss<dim>(fe.degree + 1),
860 *   mass_matrix);
861 *   MatrixCreator::create_laplace_matrix(dof_handler,
862 *   QGauss<dim>(fe.degree + 1),
863 *   laplace_matrix);
864 *  
865 * @endcode
866 *
867 * The rest of the function is spent on setting vector sizes to the
868 * correct value. The final line closes the hanging node constraints
869 * object. Since we work on a uniformly refined mesh, no constraints exist
870 * or have been computed (i.e. there was no need to call
871 * DoFTools::make_hanging_node_constraints as in other programs), but we
872 * need a constraints object in one place further down below anyway.
873 *
874 * @code
875 *   solution_u.reinit(dof_handler.n_dofs());
876 *   solution_v.reinit(dof_handler.n_dofs());
877 *   old_solution_u.reinit(dof_handler.n_dofs());
878 *   old_solution_v.reinit(dof_handler.n_dofs());
879 *   system_rhs.reinit(dof_handler.n_dofs());
880 *  
881 *   constraints.close();
882 *   }
883 *  
884 *  
885 *  
886 * @endcode
887 *
888 *
889 * <a name="step_23-WaveEquationsolve_uandWaveEquationsolve_v"></a>
890 * <h4>WaveEquation::solve_u and WaveEquation::solve_v</h4>
891 *
892
893 *
894 * The next two functions deal with solving the linear systems associated
895 * with the equations for @f$U^n@f$ and @f$V^n@f$. Both are not particularly
896 * interesting as they pretty much follow the scheme used in all the
897 * previous tutorial programs.
898 *
899
900 *
901 * One can make little experiments with preconditioners for the two matrices
902 * we have to invert. As it turns out, however, for the matrices at hand
903 * here, using Jacobi or SSOR preconditioners reduces the number of
904 * iterations necessary to solve the linear system slightly, but due to the
905 * cost of applying the preconditioner it is no win in terms of run-time. It
906 * is not much of a loss either, but let's keep it simple and just do
907 * without:
908 *
909 * @code
910 *   template <int dim>
911 *   void WaveEquation<dim>::solve_u()
912 *   {
913 *   SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
914 *   SolverCG<Vector<double>> cg(solver_control);
915 *  
916 *   cg.solve(matrix_u, solution_u, system_rhs, PreconditionIdentity());
917 *  
918 *   std::cout << " u-equation: " << solver_control.last_step()
919 *   << " CG iterations." << std::endl;
920 *   }
921 *  
922 *  
923 *  
924 *   template <int dim>
925 *   void WaveEquation<dim>::solve_v()
926 *   {
927 *   SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
928 *   SolverCG<Vector<double>> cg(solver_control);
929 *  
930 *   cg.solve(matrix_v, solution_v, system_rhs, PreconditionIdentity());
931 *  
932 *   std::cout << " v-equation: " << solver_control.last_step()
933 *   << " CG iterations." << std::endl;
934 *   }
935 *  
936 *  
937 *  
938 * @endcode
939 *
940 *
941 * <a name="step_23-WaveEquationoutput_results"></a>
942 * <h4>WaveEquation::output_results</h4>
943 *
944
945 *
946 * Likewise, the following function is pretty much what we've done
947 * before. The only thing worth mentioning is how here we generate a string
948 * representation of the time step number padded with leading zeros to 3
949 * character length using the Utilities::int_to_string function's second
950 * argument.
951 *
952 * @code
953 *   template <int dim>
954 *   void WaveEquation<dim>::output_results() const
955 *   {
956 *   DataOut<dim> data_out;
957 *  
958 *   data_out.attach_dof_handler(dof_handler);
959 *   data_out.add_data_vector(solution_u, "U");
960 *   data_out.add_data_vector(solution_v, "V");
961 *  
962 *   data_out.build_patches();
963 *  
964 *   const std::string filename =
965 *   "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtu";
966 * @endcode
967 *
968 * Like @ref step_15 "step-15", since we write output at every time step (and the system
969 * we have to solve is relatively easy), we instruct DataOut to use the
970 * zlib compression algorithm that is optimized for speed instead of disk
971 * usage since otherwise plotting the output becomes a bottleneck:
972 *
973 * @code
974 *   DataOutBase::VtkFlags vtk_flags;
976 *   data_out.set_flags(vtk_flags);
977 *   std::ofstream output(filename);
978 *   data_out.write_vtu(output);
979 *   }
980 *  
981 *  
982 *  
983 * @endcode
984 *
985 *
986 * <a name="step_23-WaveEquationrun"></a>
987 * <h4>WaveEquation::run</h4>
988 *
989
990 *
991 * The following is really the only interesting function of the program. It
992 * contains the loop over all time steps, but before we get to that we have
993 * to set up the grid, DoFHandler, and matrices. In addition, we have to
994 * somehow get started with initial values. To this end, we use the
995 * VectorTools::project function that takes an object that describes a
996 * continuous function and computes the @f$L^2@f$ projection of this function
997 * onto the finite element space described by the DoFHandler object. Can't
998 * be any simpler than that:
999 *
1000 * @code
1001 *   template <int dim>
1002 *   void WaveEquation<dim>::run()
1003 *   {
1004 *   setup_system();
1005 *  
1006 *   VectorTools::project(dof_handler,
1007 *   constraints,
1008 *   QGauss<dim>(fe.degree + 1),
1009 *   InitialValuesU<dim>(),
1010 *   old_solution_u);
1011 *   VectorTools::project(dof_handler,
1012 *   constraints,
1013 *   QGauss<dim>(fe.degree + 1),
1014 *   InitialValuesV<dim>(),
1015 *   old_solution_v);
1016 *  
1017 * @endcode
1018 *
1019 * The next thing is to loop over all the time steps until we reach the
1020 * end time (@f$T=5@f$ in this case). In each time step, we first have to
1021 * solve for @f$U^n@f$, using the equation @f$(M^n + k^2\theta^2 A^n)U^n =@f$
1022 * @f$(M^{n,n-1} - k^2\theta(1-\theta) A^{n,n-1})U^{n-1} + kM^{n,n-1}V^{n-1}
1023 * +@f$ @f$k\theta \left[k \theta F^n + k(1-\theta) F^{n-1} \right]@f$. Note
1024 * that we use the same mesh for all time steps, so that @f$M^n=M^{n,n-1}=M@f$
1025 * and @f$A^n=A^{n,n-1}=A@f$. What we therefore have to do first is to add up
1026 * @f$MU^{n-1} - k^2\theta(1-\theta) AU^{n-1} + kMV^{n-1}@f$ and the forcing
1027 * terms, and put the result into the <code>system_rhs</code> vector. (For
1028 * these additions, we need a temporary vector that we declare before the
1029 * loop to avoid repeated memory allocations in each time step.)
1030 *
1031
1032 *
1033 * The one thing to realize here is how we communicate the time variable
1034 * to the object describing the right hand side: each object derived from
1035 * the Function class has a time field that can be set using the
1036 * Function::set_time and read by Function::get_time. In essence, using
1037 * this mechanism, all functions of space and time are therefore
1038 * considered functions of space evaluated at a particular time. This
1039 * matches well what we typically need in finite element programs, where
1040 * we almost always work on a single time step at a time, and where it
1041 * never happens that, for example, one would like to evaluate a
1042 * space-time function for all times at any given spatial location.
1043 *
1044 * @code
1045 *   Vector<double> tmp(solution_u.size());
1046 *   Vector<double> forcing_terms(solution_u.size());
1047 *  
1048 *   for (; time <= 5; time += time_step, ++timestep_number)
1049 *   {
1050 *   std::cout << "Time step " << timestep_number << " at t=" << time
1051 *   << std::endl;
1052 *  
1053 *   mass_matrix.vmult(system_rhs, old_solution_u);
1054 *  
1055 *   mass_matrix.vmult(tmp, old_solution_v);
1056 *   system_rhs.add(time_step, tmp);
1057 *  
1058 *   laplace_matrix.vmult(tmp, old_solution_u);
1059 *   system_rhs.add(-theta * (1 - theta) * time_step * time_step, tmp);
1060 *  
1061 *   RightHandSide<dim> rhs_function;
1062 *   rhs_function.set_time(time);
1063 *   VectorTools::create_right_hand_side(dof_handler,
1064 *   QGauss<dim>(fe.degree + 1),
1065 *   rhs_function,
1066 *   tmp);
1067 *   forcing_terms = tmp;
1068 *   forcing_terms *= theta * time_step;
1069 *  
1070 *   rhs_function.set_time(time - time_step);
1071 *   VectorTools::create_right_hand_side(dof_handler,
1072 *   QGauss<dim>(fe.degree + 1),
1073 *   rhs_function,
1074 *   tmp);
1075 *  
1076 *   forcing_terms.add((1 - theta) * time_step, tmp);
1077 *  
1078 *   system_rhs.add(theta * time_step, forcing_terms);
1079 *  
1080 * @endcode
1081 *
1082 * After so constructing the right hand side vector of the first
1083 * equation, all we have to do is apply the correct boundary
1084 * values. As for the right hand side, this is a space-time function
1085 * evaluated at a particular time, which we interpolate at boundary
1086 * nodes and then use the result to apply boundary values as we
1087 * usually do. The result is then handed off to the solve_u()
1088 * function:
1089 *
1090 * @code
1091 *   {
1092 *   BoundaryValuesU<dim> boundary_values_u_function;
1093 *   boundary_values_u_function.set_time(time);
1094 *  
1095 *   std::map<types::global_dof_index, double> boundary_values;
1096 *   VectorTools::interpolate_boundary_values(dof_handler,
1097 *   0,
1098 *   boundary_values_u_function,
1099 *   boundary_values);
1100 *  
1101 * @endcode
1102 *
1103 * The matrix for solve_u() is the same in every time steps, so one
1104 * could think that it is enough to do this only once at the
1105 * beginning of the simulation. However, since we need to apply
1106 * boundary values to the linear system (which eliminate some matrix
1107 * rows and columns and give contributions to the right hand side),
1108 * we have to refill the matrix in every time steps before we
1109 * actually apply boundary data. The actual content is very simple:
1110 * it is the sum of the mass matrix and a weighted Laplace matrix:
1111 *
1112 * @code
1113 *   matrix_u.copy_from(mass_matrix);
1114 *   matrix_u.add(theta * theta * time_step * time_step, laplace_matrix);
1115 *   MatrixTools::apply_boundary_values(boundary_values,
1116 *   matrix_u,
1117 *   solution_u,
1118 *   system_rhs);
1119 *   }
1120 *   solve_u();
1121 *  
1122 *  
1123 * @endcode
1124 *
1125 * The second step, i.e. solving for @f$V^n@f$, works similarly, except
1126 * that this time the matrix on the left is the mass matrix (which we
1127 * copy again in order to be able to apply boundary conditions, and
1128 * the right hand side is @f$MV^{n-1} - k\left[ \theta A U^n +
1129 * (1-\theta) AU^{n-1}\right]@f$ plus forcing terms. Boundary values
1130 * are applied in the same way as before, except that now we have to
1131 * use the BoundaryValuesV class:
1132 *
1133 * @code
1134 *   laplace_matrix.vmult(system_rhs, solution_u);
1135 *   system_rhs *= -theta * time_step;
1136 *  
1137 *   mass_matrix.vmult(tmp, old_solution_v);
1138 *   system_rhs += tmp;
1139 *  
1140 *   laplace_matrix.vmult(tmp, old_solution_u);
1141 *   system_rhs.add(-time_step * (1 - theta), tmp);
1142 *  
1143 *   system_rhs += forcing_terms;
1144 *  
1145 *   {
1146 *   BoundaryValuesV<dim> boundary_values_v_function;
1147 *   boundary_values_v_function.set_time(time);
1148 *  
1149 *   std::map<types::global_dof_index, double> boundary_values;
1150 *   VectorTools::interpolate_boundary_values(dof_handler,
1151 *   0,
1152 *   boundary_values_v_function,
1153 *   boundary_values);
1154 *   matrix_v.copy_from(mass_matrix);
1155 *   MatrixTools::apply_boundary_values(boundary_values,
1156 *   matrix_v,
1157 *   solution_v,
1158 *   system_rhs);
1159 *   }
1160 *   solve_v();
1161 *  
1162 * @endcode
1163 *
1164 * Finally, after both solution components have been computed, we
1165 * output the result, compute the energy in the solution, and go on to
1166 * the next time step after shifting the present solution into the
1167 * vectors that hold the solution at the previous time step. Note the
1168 * function SparseMatrix::matrix_norm_square that can compute
1169 * @f$\left<V^n,MV^n\right>@f$ and @f$\left<U^n,AU^n\right>@f$ in one step,
1170 * saving us the expense of a temporary vector and several lines of
1171 * code:
1172 *
1173 * @code
1174 *   output_results();
1175 *  
1176 *   std::cout << " Total energy: "
1177 *   << (mass_matrix.matrix_norm_square(solution_v) +
1178 *   laplace_matrix.matrix_norm_square(solution_u)) /
1179 *   2
1180 *   << std::endl;
1181 *  
1182 *   old_solution_u = solution_u;
1183 *   old_solution_v = solution_v;
1184 *   }
1185 *   }
1186 *   } // namespace Step23
1187 *  
1188 *  
1189 * @endcode
1190 *
1191 *
1192 * <a name="step_23-Thecodemaincodefunction"></a>
1193 * <h3>The <code>main</code> function</h3>
1194 *
1195
1196 *
1197 * What remains is the main function of the program. There is nothing here
1198 * that hasn't been shown in several of the previous programs:
1199 *
1200 * @code
1201 *   int main()
1202 *   {
1203 *   try
1204 *   {
1205 *   using namespace Step23;
1206 *  
1207 *   WaveEquation<2> wave_equation_solver;
1208 *   wave_equation_solver.run();
1209 *   }
1210 *   catch (std::exception &exc)
1211 *   {
1212 *   std::cerr << std::endl
1213 *   << std::endl
1214 *   << "----------------------------------------------------"
1215 *   << std::endl;
1216 *   std::cerr << "Exception on processing: " << std::endl
1217 *   << exc.what() << std::endl
1218 *   << "Aborting!" << std::endl
1219 *   << "----------------------------------------------------"
1220 *   << std::endl;
1221 *  
1222 *   return 1;
1223 *   }
1224 *   catch (...)
1225 *   {
1226 *   std::cerr << std::endl
1227 *   << std::endl
1228 *   << "----------------------------------------------------"
1229 *   << std::endl;
1230 *   std::cerr << "Unknown exception!" << std::endl
1231 *   << "Aborting!" << std::endl
1232 *   << "----------------------------------------------------"
1233 *   << std::endl;
1234 *   return 1;
1235 *   }
1236 *  
1237 *   return 0;
1238 *   }
1239 * @endcode
1240<a name="step_23-Results"></a><h1>Results</h1>
1241
1242
1243When the program is run, it produces the following output:
1244@code
1245Number of active cells: 16384
1246Number of degrees of freedom: 16641
1247
1248Time step 1 at t=0.015625
1249 u-equation: 8 CG iterations.
1250 v-equation: 22 CG iterations.
1251 Total energy: 1.17887
1252Time step 2 at t=0.03125
1253 u-equation: 8 CG iterations.
1254 v-equation: 20 CG iterations.
1255 Total energy: 2.9655
1256Time step 3 at t=0.046875
1257 u-equation: 8 CG iterations.
1258 v-equation: 21 CG iterations.
1259 Total energy: 4.33761
1260Time step 4 at t=0.0625
1261 u-equation: 7 CG iterations.
1262 v-equation: 21 CG iterations.
1263 Total energy: 5.35499
1264Time step 5 at t=0.078125
1265 u-equation: 7 CG iterations.
1266 v-equation: 21 CG iterations.
1267 Total energy: 6.18652
1268Time step 6 at t=0.09375
1269 u-equation: 7 CG iterations.
1270 v-equation: 20 CG iterations.
1271 Total energy: 6.6799
1272
1273...
1274
1275Time step 31 at t=0.484375
1276 u-equation: 7 CG iterations.
1277 v-equation: 20 CG iterations.
1278 Total energy: 21.9068
1279Time step 32 at t=0.5
1280 u-equation: 7 CG iterations.
1281 v-equation: 20 CG iterations.
1282 Total energy: 23.3394
1283Time step 33 at t=0.515625
1284 u-equation: 7 CG iterations.
1285 v-equation: 20 CG iterations.
1286 Total energy: 23.1019
1287
1288...
1289
1290Time step 319 at t=4.98438
1291 u-equation: 7 CG iterations.
1292 v-equation: 20 CG iterations.
1293 Total energy: 23.1019
1294Time step 320 at t=5
1295 u-equation: 7 CG iterations.
1296 v-equation: 20 CG iterations.
1297 Total energy: 23.1019
1298@endcode
1299
1300What we see immediately is that the energy is a constant at least after
1301@f$t=\frac 12@f$ (until which the boundary source term @f$g@f$ is nonzero, injecting
1302energy into the system).
1303
1304In addition to the screen output, the program writes the solution of each time
1305step to an output file. If we process them adequately and paste them into a
1306movie, we get the following:
1307
1308<img src="https://www.dealii.org/images/steps/developer/step-23.movie.gif" alt="Animation of the solution of step 23.">
1309
1310The movie shows the generated wave nice traveling through the domain and back,
1311being reflected at the clamped boundary. Some numerical noise is trailing the
1312wave, an artifact of a too-large mesh size that can be reduced by reducing the
1313mesh width and the time step.
1314
1315
1316<a name="step-23-extensions"></a>
1317<a name="step_23-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1318
1319
1320If you want to explore a bit, try out some of the following things:
1321<ul>
1322 <li>Varying @f$\theta@f$. This gives different time stepping schemes, some of
1323 which are stable while others are not. Take a look at how the energy
1324 evolves.
1325
1326 <li>Different initial and boundary conditions, right hand sides.
1327
1328 <li>More complicated domains or more refined meshes. Remember that the time
1329 step needs to be bounded by the mesh width, so changing the mesh should
1330 always involve also changing the time step. We will come back to this issue
1331 in @ref step_24 "step-24".
1332
1333 <li>Variable coefficients: In real media, the wave speed is often
1334 variable. In particular, the "real" wave equation in realistic media would
1335 read
1336 @f[
1337 \rho(x) \frac{\partial^2 u}{\partial t^2}
1338 -
1339 \nabla \cdot
1340 a(x) \nabla u = f,
1341 @f]
1342 where @f$\rho(x)@f$ is the density of the material, and @f$a(x)@f$ is related to the
1343 stiffness coefficient. The wave speed is then @f$c=\sqrt{a/\rho}@f$.
1344
1345 To make such a change, we would have to compute the mass and Laplace
1346 matrices with a variable coefficient. Fortunately, this isn't too hard: the
1347 functions MatrixCreator::create_laplace_matrix and
1348 MatrixCreator::create_mass_matrix have additional default parameters that can
1349 be used to pass non-constant coefficient functions to them. The required
1350 changes are therefore relatively small. On the other hand, care must be
1351 taken again to make sure the time step is within the allowed range.
1352
1353 <li>In the in-code comments, we discussed the fact that the matrices for
1354 solving for @f$U^n@f$ and @f$V^n@f$ need to be reset in every time because of
1355 boundary conditions, even though the actual content does not change. It is
1356 possible to avoid copying by not eliminating columns in the linear systems,
1357 which is implemented by appending a @p false argument to the call:
1358 @code
1359 MatrixTools::apply_boundary_values(boundary_values,
1360 matrix_u,
1361 solution_u,
1362 system_rhs,
1363 false);
1364 @endcode
1365
1366 <li>deal.II being a library that supports adaptive meshes it would of course be
1367 nice if this program supported change the mesh every few time steps. Given the
1368 structure of the solution &mdash; a wave that travels through the domain &mdash;
1369 it would seem appropriate if we only refined the mesh where the wave currently is,
1370 and not simply everywhere. It is intuitively clear that we should be able to
1371 save a significant amount of cells this way. (Though upon further thought one
1372 realizes that this is really only the case in the initial stages of the simulation.
1373 After some time, for wave phenomena, the domain is filled with reflections of
1374 the initial wave going in every direction and filling every corner of the domain.
1375 At this point, there is in general little one can gain using local mesh
1376 refinement.)
1377
1378 To make adaptively changing meshes possible, there are basically two routes.
1379 The "correct" way would be to go back to the weak form we get using Rothe's
1380 method. For example, the first of the two equations to be solved in each time
1381 step looked like this:
1382 \f{eqnarray*}
1383 (u^n,\varphi) + k^2\theta^2(\nabla u^n,\nabla \varphi) &=&
1384 (u^{n-1},\varphi) - k^2\theta(1-\theta)(\nabla u^{n-1},\nabla \varphi)
1385 +
1386 k(v^{n-1},\varphi)
1387 + k^2\theta
1388 \left[
1389 \theta (f^n,\varphi) + (1-\theta) (f^{n-1},\varphi)
1390 \right].
1391 \f}
1392 Now, note that we solve for @f$u^n@f$ on mesh @f${\mathbb T}^n@f$, and
1393 consequently the test functions @f$\varphi@f$ have to be from the space
1394 @f$V_h^n@f$ as well. As discussed in the introduction, terms like
1395 @f$(u^{n-1},\varphi)@f$ then require us to integrate the solution of the
1396 previous step (which may have been computed on a different mesh
1397 @f${\mathbb T}^{n-1}@f$) against the test functions of the current mesh,
1398 leading to a matrix @f$M^{n,n-1}@f$. This process of integrating shape
1399 functions from different meshes is, at best, awkward. It can be done
1400 but because it is difficult to ensure that @f${\mathbb T}^{n-1}@f$ and
1401 @f${\mathbb T}^{n}@f$ differ by at most one level of refinement, one
1402 has to recursively match cells from both meshes. It is feasible to
1403 do this, but it leads to lengthy and not entirely obvious code.
1404
1405 The second approach is the following: whenever we change the mesh,
1406 we simply interpolate the solution from the last time step on the old
1407 mesh to the new mesh, using the SolutionTransfer class. In other words,
1408 instead of the equation above, we would solve
1409 \f{eqnarray*}
1410 (u^n,\varphi) + k^2\theta^2(\nabla u^n,\nabla \varphi) &=&
1411 (I^n u^{n-1},\varphi) - k^2\theta(1-\theta)(\nabla I^n u^{n-1},\nabla \varphi)
1412 +
1413 k(I^n v^{n-1},\varphi)
1414 + k^2\theta
1415 \left[
1416 \theta (f^n,\varphi) + (1-\theta) (f^{n-1},\varphi)
1417 \right],
1418 \f}
1419 where @f$I^n@f$ interpolates a given function onto mesh @f${\mathbb T}^n@f$.
1420 This is a much simpler approach because, in each time step, we no
1421 longer have to worry whether @f$u^{n-1},v^{n-1}@f$ were computed on the
1422 same mesh as we are using now or on a different mesh. Consequently,
1423 the only changes to the code necessary are the addition of a function
1424 that computes the error, marks cells for refinement, sets up a
1425 SolutionTransfer object, transfers the solution to the new mesh, and
1426 rebuilds matrices and right hand side vectors on the new mesh. Neither
1427 the functions building the matrices and right hand sides, nor the
1428 solvers need to be changed.
1429
1430 While this second approach is, strictly speaking,
1431 not quite correct in the Rothe framework (it introduces an addition source
1432 of error, namely the interpolation), it is nevertheless what
1433 almost everyone solving time dependent equations does. We will use this
1434 method in @ref step_31 "step-31", for example.
1435</ul>
1436 *
1437 *
1438<a name="step_23-PlainProg"></a>
1439<h1> The plain program</h1>
1440@include "step-23.cc"
1441*/
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
Number get_time() const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Definition point.h:111
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
unsigned int level
Definition grid_out.cc:4626
#define Assert(cond, exc)
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
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_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)
const Event initial
Definition event.cc:64
void interpolate(const DoFHandler< dim, spacedim > &dof1, const InVector &u1, const DoFHandler< dim, spacedim > &dof2, OutVector &u2)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
void project(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const AffineConstraints< typename VectorType::value_type > &constraints, const Quadrature< dim > &quadrature, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const bool enforce_zero_boundary=false, const Quadrature< dim - 1 > &q_boundary=(dim > 1 ? QGauss< dim - 1 >(2) :Quadrature< dim - 1 >()), const bool project_to_boundary_first=false)
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)
static constexpr double PI
Definition numbers.h:254
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DataOutBase::CompressionLevel compression_level
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)