From d5e9078a68eb027a95870414754e99a091c7f588 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 23 May 2008 13:04:06 +0000 Subject: [PATCH] Corrected a few typos in intro and in-code comments. git-svn-id: https://svn.dealii.org/trunk@16177 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-33/doc/intro.dox | 38 ++-- deal.II/examples/step-33/step-33.cc | 260 ++++++++++++------------- 2 files changed, 149 insertions(+), 149 deletions(-) diff --git a/deal.II/examples/step-33/doc/intro.dox b/deal.II/examples/step-33/doc/intro.dox index 256a8a18fe..980f632b4d 100644 --- a/deal.II/examples/step-33/doc/intro.dox +++ b/deal.II/examples/step-33/doc/intro.dox @@ -3,7 +3,7 @@ This program was written for fun by David Neckels (NCAR) while working at Sandia (on the Wyoming Express bus to and from Corrales each day). -The main purpose was to better understand Euler flow. +The main purpose was to better understand Euler flow. The code solves the basic Euler equations of gas dynamics, by using a fully implicit Newton iteration (inspired by Sandia's Aria code). The code may be configured by an input deck to run different simulations @@ -39,7 +39,7 @@ a basic system of conservation laws. In spatial dimension $d$ they read with the solution $\mathbf{w}=(\rho v_1,\ldots,\rho v_d,\rho, E)^{\top}$ consisting of $\rho$ the fluid density, ${\mathbf v}=(v_1,\ldots v_d)^T$ the flow velocity (and thus $\rho\mathbf v$ being the linear momentum -density), and +density), and $E$ the energy density of the gas. We interpret the equations above as $\partial_t \mathbf{w}_i + \nabla \cdot \mathbf{F}_i(\mathbf{w}) = \mathbf G_i(\mathbf w)$, $i=1,\ldots,dim+2$. @@ -80,7 +80,7 @@ With this, the entire system of equations reads: \partial_t (\rho v_i) + \sum_{s=1}^d \frac{\partial(\rho v_i v_s + \delta_{is} p)}{\partial x_s} &=& g_i \rho, \qquad i=1,\dots,d, \\ \partial_t \rho + \sum_{s=1}^d \frac{\partial(\rho v_s)}{\partial x_s} &=& 0, \\ - \partial_t E + \sum_{s=1}^d \frac{\partial((E+p)v_s)}{\partial x_s} &=& + \partial_t E + \sum_{s=1}^d \frac{\partial((E+p)v_s)}{\partial x_s} &=& \rho \mathbf g \cdot \mathbf v. @f} These equations describe, respectively, the conservation of momentum, @@ -114,10 +114,10 @@ The diffusion term $h^{\eta}(\nabla \mathbf{w} , \nabla \mathbf{z})$ is introduc On the boundary, we have to say what the outer trace $\mathbf{w}^-$ is. Depending on the boundary condition, we prescribe either of the following: The reason for this approach is that it separates the various concerns in a @@ -331,7 +331,7 @@ would be relatively straightforward to adapt it to a different set of equations: One would simply re-implement the members of the EulerEquations class for some other hyperbolic equation, or augment the existing equations by additional ones (for example by advecting -additional variables, or my adding chemistry, etc). Such modifications, +additional variables, or by adding chemistry, etc). Such modifications, however, would not affect the time stepping, or the nonlinear solvers if correctly done, and consequently nothing in the ConservationLaw would have to be modified. diff --git a/deal.II/examples/step-33/step-33.cc b/deal.II/examples/step-33/step-33.cc index 4db6ef786f..e2e2c0ab2d 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -110,7 +110,7 @@ using namespace dealii; // introduction. We group all this into a // structure that defines everything that has // to do with the flux. All members of this - // structures are static, i.e. the structure + // structure are static, i.e. the structure // has no actual state specified by instance // member variables. The better way to do // this, rather than a structure with all @@ -124,7 +124,7 @@ template struct EulerEquations { // @sect4{Component description} - + // First a few variables that // describe the various components of our // solution vector in a generic way. This @@ -188,13 +188,13 @@ struct EulerEquations .push_back (DataComponentInterpretation::component_is_scalar); data_component_interpretation .push_back (DataComponentInterpretation::component_is_scalar); - + return data_component_interpretation; } - - + + // @sect4{Transformations between variables} - + // Next, we define the gas // constant. We will set it to 1.4 // in its definition immediately @@ -270,11 +270,11 @@ struct EulerEquations return ((gas_gamma-1.0) * (*(W.begin() + energy_component) - compute_kinetic_energy(W))); - } + } + + // @sect4{EulerEquations::compute_flux_matrix} - // @sect4{EulerEquations::compute_flux_matrix} - // We define the flux function // $F(W)$ as one large matrix. // Each row of this matrix @@ -317,7 +317,7 @@ struct EulerEquations // matrix that correspond to the // momentum terms: const number pressure = compute_pressure (W); - + for (unsigned int d=0; d iflux[n_components][dim]; Sacado::Fad::DFad oflux[n_components][dim]; - + compute_flux_matrix (Wplus, iflux); compute_flux_matrix (Wminus, oflux); - + for (unsigned int di=0; di &dof_handler, const Mapping &mapping, const Vector &solution, - Vector &refinement_indicators) - { + Vector &refinement_indicators) + { const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; std::vector dofs (dofs_per_cell); @@ -634,7 +634,7 @@ struct EulerEquations std::vector > > dU (1, std::vector >(n_components)); - + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -647,13 +647,13 @@ struct EulerEquations = std::log(1+ std::sqrt(dU[0][density_component] * dU[0][density_component])); - } + } } - - - + + + // @sect4{EulerEquations::Postprocessor} - + // Finally, we declare a class that // implements a postprocessing of data // components. The problem this class @@ -673,14 +673,14 @@ struct EulerEquations // possibility to generate schlieren // plots. Schlieren plots are a way to // visualize shocks and other sharp - // interfaces. The word "schlieren" a + // interfaces. The word "schlieren" is a // German word that may be translated as // "striae" -- it may be simpler to // explain it by an example, however: // schlieren is what you see when you, // for example, pour highly concentrated // alcohol, or a transparent saline - // solution into water; the two have the + // solution, into water; the two have the // same color, but they have different // refractive indices and so before they // are fully mixed light goes through the @@ -688,7 +688,7 @@ struct EulerEquations // brightness variations if you look at // it. That's "schlieren". A similar // effect happens in compressible flow - // due because the refractive index + // because the refractive index // depends on the pressure (and therefore // the density) of the gas. // @@ -722,7 +722,7 @@ struct EulerEquations { public: Postprocessor (const bool do_schlieren_plot); - + virtual void compute_derived_quantities_vector (const std::vector > &uh, @@ -732,11 +732,11 @@ struct EulerEquations std::vector > &computed_quantities) const; virtual std::vector get_names () const; - + virtual std::vector get_data_component_interpretation () const; - + virtual UpdateFlags get_needed_update_flags () const; virtual unsigned int n_output_variables() const; @@ -745,7 +745,7 @@ struct EulerEquations const bool do_schlieren_plot; }; }; - + template const double EulerEquations::gas_gamma = 1.4; @@ -895,7 +895,7 @@ UpdateFlags EulerEquations::Postprocessor:: get_needed_update_flags () const { - if (do_schlieren_plot == true) + if (do_schlieren_plot == true) return update_values | update_gradients; else return update_values; @@ -1020,7 +1020,7 @@ namespace Parameters // perturbation that is added to the // diagonal before forming the prec, // and RTOL is a scaling factor $rtol - // >= 1$. + // \geq 1$. // - ilut_drop: The ILUT will // drop any values that // have magnitude less than this value. @@ -1034,11 +1034,11 @@ namespace Parameters // ParameterHandler::declare_entry // call in // declare_parameters(). - struct Solver + struct Solver { enum SolverType { gmres, direct }; SolverType solver; - + enum OutputType { quiet, verbose }; OutputType output; @@ -1089,10 +1089,10 @@ namespace Parameters } prm.leave_subsection(); } - - - + + + void Solver::parse_parameters (ParameterHandler &prm) { prm.enter_subsection("linear solver"); @@ -1102,9 +1102,9 @@ namespace Parameters output = verbose; if (op == "quiet") output = quiet; - + const std::string sv = prm.get("method"); - if (sv == "direct") + if (sv == "direct") solver = direct; else if (sv == "gmres") solver = gmres; @@ -1116,11 +1116,11 @@ namespace Parameters ilut_rtol = prm.get_double("ilut relative tolerance"); ilut_drop = prm.get_double("ilut drop tolerance"); } - prm.leave_subsection(); + prm.leave_subsection(); } - - + + // @sect4{Parameters::Refinement} // // Similarly, here are a few parameters @@ -1167,7 +1167,7 @@ namespace Parameters } prm.leave_subsection(); } - + void Refinement::parse_parameters (ParameterHandler &prm) { @@ -1179,7 +1179,7 @@ namespace Parameters } prm.leave_subsection(); } - + // @sect4{Parameters::Flux} @@ -1199,13 +1199,13 @@ namespace Parameters // dependent value. In the latter case, it // is chosen as $\frac{h}{2\delta T}$ with // $h$ the diameter of the face to which - // the flux is applied, and $\delta T$ + // the flux is applied, and $\delta T$ // the current time step. struct Flux { enum StabilizationKind { constant, mesh_dependent }; StabilizationKind stabilization_kind; - + double stabilization_value; static void declare_parameters (ParameterHandler &prm); @@ -1225,22 +1225,22 @@ namespace Parameters Patterns::Double(), "alpha stabilization"); } - prm.leave_subsection(); + prm.leave_subsection(); } - - + + void Flux::parse_parameters (ParameterHandler &prm) { prm.enter_subsection("flux"); { const std::string stab = prm.get("stab"); - if (stab == "constant") + if (stab == "constant") stabilization_kind = constant; - else if (stab == "mesh") + else if (stab == "mesh") stabilization_kind = mesh_dependent; else AssertThrow (false, ExcNotImplemented()); - + stabilization_value = prm.get_double("stab value"); } prm.leave_subsection(); @@ -1271,7 +1271,7 @@ namespace Parameters void Output::declare_parameters (ParameterHandler &prm) { prm.enter_subsection("output"); - { + { prm.declare_entry("schlieren plot", "true", Patterns::Bool (), "Whether or not to produce schlieren plots"); @@ -1281,7 +1281,7 @@ namespace Parameters } prm.leave_subsection(); } - + void Output::parse_parameters (ParameterHandler &prm) @@ -1396,15 +1396,15 @@ namespace Parameters { typename EulerEquations::BoundaryKind kind[EulerEquations::n_components]; - + FunctionParser values; BoundaryConditions (); }; - - + + AllParameters (); - + double diffusion_power; double time_step, final_time; @@ -1415,7 +1415,7 @@ namespace Parameters FunctionParser initial_conditions; BoundaryConditions boundary_conditions[max_n_boundaries]; - + static void declare_parameters (ParameterHandler &prm); void parse_parameters (ParameterHandler &prm); }; @@ -1434,7 +1434,7 @@ namespace Parameters : initial_conditions (EulerEquations::n_components) {} - + template void @@ -1481,7 +1481,7 @@ namespace Parameters "outflow", Patterns::Selection("inflow|outflow|pressure"), ""); - + prm.declare_entry("w_" + Utilities::int_to_string(di) + " value", "0.0", Patterns::Anything(), @@ -1526,7 +1526,7 @@ namespace Parameters } else is_stationary = false; - + final_time = prm.get_double("final time"); theta = prm.get_double("theta scheme value"); } @@ -1540,7 +1540,7 @@ namespace Parameters { std::vector expressions(EulerEquations::n_components, "0.0"); - + const bool no_penetration = prm.get_bool("no penetration"); for (unsigned int di=0; di::n_components; ++di) @@ -1564,7 +1564,7 @@ namespace Parameters AssertThrow (false, ExcNotImplemented()); expressions[di] = prm.get("w_" + Utilities::int_to_string(di) + - " value"); + " value"); } boundary_conditions[boundary_id].values @@ -1592,11 +1592,11 @@ namespace Parameters Parameters::Refinement::parse_parameters (prm); Parameters::Flux::parse_parameters (prm); Parameters::Output::parse_parameters (prm); - } + } } - - + + // @sect3{Conservation law class} @@ -1625,7 +1625,7 @@ class ConservationLaw public: ConservationLaw (const char *input_filename); void run (); - + private: void setup_system (); @@ -1676,14 +1676,14 @@ class ConservationLaw // not of sufficiently high // order. Triangulation triangulation; - const MappingQ1 mapping; - + const MappingQ1 mapping; + const FESystem fe; DoFHandler dof_handler; const QGauss quadrature; const QGauss face_quadrature; - + // Next come a number of data // vectors that correspond to the // solution of the previous time @@ -1756,7 +1756,7 @@ class ConservationLaw Epetra_SerialComm communicator; std::auto_ptr Map; std::auto_ptr Matrix; - + Parameters::AllParameters parameters; ConditionalOStream verbose_cout; }; @@ -1821,7 +1821,7 @@ void ConservationLaw::setup_system () dof_handler.n_dofs()); DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); sparsity_pattern.compress(); - + std::vector row_lengths (dof_handler.n_dofs()); for (unsigned int i=0; i::setup_system () const unsigned int max_nonzero_entries = *std::max_element (row_lengths.begin(), row_lengths.end()); - + std::vector values(max_nonzero_entries, 0); std::vector row_indices(max_nonzero_entries); - + for (unsigned int row=0; rowInsertGlobalValues(row, row_lengths[row], &values[0], &row_indices[0]); } @@ -1918,7 +1918,7 @@ void ConservationLaw::assemble_system () | update_JxW_values | update_normal_vectors, neighbor_face_update_flags = update_values; - + FEValues fe_v (mapping, fe, quadrature, update_flags); FEFaceValues fe_v_face (mapping, fe, face_quadrature, @@ -1937,7 +1937,7 @@ void ConservationLaw::assemble_system () typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); - for (; cell!=endc; ++cell) + for (; cell!=endc; ++cell) { fe_v.reinit (cell); cell->get_dof_indices (dof_indices); @@ -2040,7 +2040,7 @@ void ConservationLaw::assemble_system () { const unsigned int neighbor2= cell->neighbor_of_neighbor(face_no); - + for (unsigned int subface_no=0; subface_no::subfaces_per_face; ++subface_no) @@ -2066,7 +2066,7 @@ void ConservationLaw::assemble_system () dof_indices_neighbor, false, numbers::invalid_unsigned_int, - neighbor_child->diameter()); + neighbor_child->diameter()); } } @@ -2105,7 +2105,7 @@ void ConservationLaw::assemble_system () fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no, neighbor_subface_no); - + assemble_face_term (face_no, fe_v_face, fe_v_subface_neighbor, dof_indices, @@ -2115,7 +2115,7 @@ void ConservationLaw::assemble_system () cell->face(face_no)->diameter()); } } - } + } // After all this assembling, notify the // Trilinos matrix object that the matrix @@ -2216,7 +2216,7 @@ template void ConservationLaw:: assemble_cell_term (const FEValues &fe_v, - const std::vector &dof_indices) + const std::vector &dof_indices) { const unsigned int dofs_per_cell = fe_v.dofs_per_cell; const unsigned int n_q_points = fe_v.n_quadrature_points; @@ -2229,12 +2229,12 @@ assemble_cell_term (const FEValues &fe_v, Table<2,Sacado::Fad::DFad > W_theta (n_q_points, EulerEquations::n_components); - + Table<3,Sacado::Fad::DFad > grad_W (n_q_points, EulerEquations::n_components, dim); std::vector residual_derivatives (dofs_per_cell); - + // Next, we have to define the independent // variables that we will try to determine // by solving a Newton step. These @@ -2244,7 +2244,7 @@ assemble_cell_term (const FEValues &fe_v, std::vector > independent_local_dof_values(dofs_per_cell); for (unsigned int i=0; i &fe_v, for (unsigned int i=0; i &fe_v, (1-parameters.theta) * old_solution(dof_indices[i])) * fe_v.shape_value_component(i, q, c); - + for (unsigned int d = 0; d < dim; d++) grad_W[q][c][d] += independent_local_dof_values[i] * fe_v.shape_grad_component(i, q, c)[d]; @@ -2343,13 +2343,13 @@ assemble_cell_term (const FEValues &fe_v, typedef Sacado::Fad::DFad ForcingVector[EulerEquations::n_components]; ForcingVector *forcing = new ForcingVector[n_q_points]; - + for (unsigned int q=0; q::compute_flux_matrix (W_theta[q], flux[q]); EulerEquations::compute_forcing_vector (W_theta[q], forcing[q]); } - + // We now have all of the pieces in place, // so perform the assembly. We have an @@ -2397,14 +2397,14 @@ assemble_cell_term (const FEValues &fe_v, // entries. Then, when we sum into the // right_hand_side vector, // we negate this residual. - for (unsigned int i=0; i F_i = 0; const unsigned int component_i = fe_v.get_fe().system_to_component_index(i).first; - // The residual for each row (i) will be accumulating + // The residual for each row (i) will be accumulating // into this fad variable. At the end of the assembly // for this row, we will query for the sensitivities // to this variable and add them into the Jacobian. @@ -2417,7 +2417,7 @@ assemble_cell_term (const FEValues &fe_v, fe_v.shape_value_component(i, point, component_i) * fe_v.JxW(point); - for (unsigned int d=0; d::assemble_face_term(const unsigned int face_no, const std::vector &dof_indices_neighbor, const bool external_face, const unsigned int boundary_id, - const double face_diameter) + const double face_diameter) { const unsigned int n_q_points = fe_v.n_quadrature_points; const unsigned int dofs_per_cell = fe_v.dofs_per_cell; @@ -2520,7 +2520,7 @@ ConservationLaw::assemble_face_term(const unsigned int face_no, const unsigned int n_independent_variables = (external_face == false ? 2 * dofs_per_cell : dofs_per_cell); - + for (unsigned int i = 0; i < dofs_per_cell; i++) { independent_local_dof_values[i] = current_solution(dof_indices[i]); @@ -2588,7 +2588,7 @@ ConservationLaw::assemble_face_term(const unsigned int face_no, // of $W^-$ will be either functions of // $W^+$, or they will be prescribed, // depending on the kind of boundary - // condition imposed here. + // condition imposed here. // // To start the evaluation, let us ensure // that the boundary id specified for this @@ -2629,7 +2629,7 @@ ConservationLaw::assemble_face_term(const unsigned int face_no, Wminus[q]); } - + // Now that we have $\mathbf w^+$ and // $\mathbf w^-$, we can go about computing // the numerical flux function $\mathbf @@ -2680,11 +2680,11 @@ ConservationLaw::assemble_face_term(const unsigned int face_no, { const unsigned int component_i = fe_v.get_fe().system_to_component_index(i).first; - + F_i += normal_fluxes[point][component_i] * fe_v.shape_value_component(i, point, component_i) * fe_v.JxW(point); - } + } for (unsigned int k=0; k::assemble_face_term(const unsigned int face_no, reinterpret_cast( const_cast( &dof_indices[0]))); - + if (external_face == false) { for (unsigned int k=0; k::assemble_face_term(const unsigned int face_no, const_cast( &dof_indices_neighbor[0]))); } - + right_hand_side(dof_indices[i]) -= F_i.val(); } @@ -2737,12 +2737,12 @@ ConservationLaw::assemble_face_term(const unsigned int face_no, template std::pair -ConservationLaw::solve (Vector &newton_update) +ConservationLaw::solve (Vector &newton_update) { Epetra_Vector x(View, *Map, newton_update.begin()); Epetra_Vector b(View, *Map, right_hand_side.begin()); - + switch (parameters.solver) { // If the parameter file specified that @@ -2823,7 +2823,7 @@ ConservationLaw::solve (Vector &newton_update) solver.TrueResidual()); } } - + Assert (false, ExcNotImplemented()); return std::make_pair (0,0); } @@ -2874,7 +2874,7 @@ ConservationLaw::refine_grid (const Vector &refinement_indicators) if ((cell->level() < parameters.shock_levels) && (std::fabs(refinement_indicators(cell_no)) > parameters.shock_val)) cell->set_refine_flag(); - else + else if ((cell->level() > 0) && (std::fabs(refinement_indicators(cell_no)) < 0.75*parameters.shock_val)) cell->set_coarsen_flag(); @@ -2897,7 +2897,7 @@ ConservationLaw::refine_grid (const Vector &refinement_indicators) transfer_in.push_back(predictor); triangulation.prepare_coarsening_and_refinement(); - + SolutionTransfer soltrans(dof_handler); soltrans.prepare_for_coarsening_and_refinement(transfer_in); @@ -2950,12 +2950,12 @@ ConservationLaw::refine_grid (const Vector &refinement_indicators) template void ConservationLaw::output_results () const { - typename EulerEquations::Postprocessor + typename EulerEquations::Postprocessor postprocessor (parameters.schlieren_plot); DataOut data_out; data_out.attach_dof_handler (dof_handler); - + data_out.add_data_vector (current_solution, EulerEquations::component_names (), DataOut::type_dof_data, @@ -2995,7 +2995,7 @@ void ConservationLaw::output_results () const // solution. At the end of this process, we // output the initial solution. template -void ConservationLaw::run () +void ConservationLaw::run () { { GridIn grid_in; @@ -3004,12 +3004,12 @@ void ConservationLaw::run () std::ifstream input_file(parameters.mesh_filename.c_str()); Assert (input_file, ExcFileNotOpen(parameters.mesh_filename.c_str())); - grid_in.read_ucd(input_file); + grid_in.read_ucd(input_file); } - + dof_handler.clear(); dof_handler.distribute_dofs (fe); - + // Size all of the fields. old_solution.reinit (dof_handler.n_dofs()); current_solution.reinit (dof_handler.n_dofs()); @@ -3030,7 +3030,7 @@ void ConservationLaw::run () compute_refinement_indicators(refinement_indicators); refine_grid(refinement_indicators); - + setup_system(); VectorTools::interpolate(dof_handler, @@ -3064,7 +3064,7 @@ void ConservationLaw::run () << dof_handler.n_dofs() << std::endl << std::endl; - + std::cout << " NonLin Res Lin Iter Lin Res" << std::endl << " _____________________________________" << std::endl; @@ -3118,10 +3118,10 @@ void ConservationLaw::run () { Matrix->PutScalar(0); Matrix->FillComplete(); - + right_hand_side = 0; assemble_system (); - + const double res_norm = right_hand_side.l2_norm(); if (std::fabs(res_norm) < 1e-10) { @@ -3134,9 +3134,9 @@ void ConservationLaw::run () std::pair convergence = solve (newton_update); - + current_solution += newton_update; - + std::printf(" %-16.3e %04d %-5.2e\n", res_norm, convergence.first, convergence.second); } @@ -3144,7 +3144,7 @@ void ConservationLaw::run () ++nonlin_iter; AssertThrow (nonlin_iter <= 10, ExcMessage ("No convergence in nonlinear solver")); - } + } // We only get to this point if the // Newton iteration has converged, so @@ -3214,14 +3214,14 @@ void ConservationLaw::run () // that the program aborts if no // input file name is given on the // command line. -int main (int argc, char *argv[]) +int main (int argc, char *argv[]) { if (argc != 2) { std::cout << "Usage:" << argv[0] << " infile" << std::endl; std::exit(1); } - + try { ConservationLaw<2> cons (argv[1]); @@ -3239,7 +3239,7 @@ int main (int argc, char *argv[]) << std::endl; return 1; } - catch (...) + catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------" @@ -3250,7 +3250,7 @@ int main (int argc, char *argv[]) << std::endl; return 1; }; - + return 0; } -- 2.39.5