From 97bd72f49631b2d73e7fe55b327f6e78a7486ed9 Mon Sep 17 00:00:00 2001 From: bangerth Date: Tue, 20 May 2008 20:49:10 +0000 Subject: [PATCH] Remove gravity from input files. Move handling boundary conditions to a function in EulerEquations. git-svn-id: https://svn.dealii.org/trunk@16146 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-33/doc/intro.dox | 5 + deal.II/examples/step-33/doc/results.dox | 3 - deal.II/examples/step-33/input.prm | 3 - deal.II/examples/step-33/step-33.cc | 542 ++++++++++++++--------- 4 files changed, 348 insertions(+), 205 deletions(-) diff --git a/deal.II/examples/step-33/doc/intro.dox b/deal.II/examples/step-33/doc/intro.dox index ef696f4c25..f4f6999dbb 100644 --- a/deal.II/examples/step-33/doc/intro.dox +++ b/deal.II/examples/step-33/doc/intro.dox @@ -8,6 +8,11 @@ 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 on different meshes, with differing boundary conditions. +
+The original code and documentation was later slightly modified by Wolfgang +Bangerth to make it more modular and allow replacing the parts that are +specific to the Euler equations by other hyperbolic conservation laws without +too much trouble. Note:The program uses the struct EulerEquations { + // @sect4{Component description} + // First a few variables that // describe the various components of our // solution vector in a generic way. This @@ -188,6 +190,8 @@ struct EulerEquations } + // @sect4{Transformations between variables} + // Next, we define the gas // constant. We will set it to 1.4 // in its definition immediately @@ -206,7 +210,7 @@ struct EulerEquations // and $O_2$. static const double gas_gamma; - + // In the following, we will need to // compute the kinetic energy and the // pressure from a vector of conserved @@ -266,6 +270,8 @@ struct EulerEquations } + // @sect4{EulerEquations::compute_flux_matrix} + // We define the flux function // $F(W)$ as one large matrix. // Each row of this matrix @@ -335,6 +341,8 @@ struct EulerEquations } + // @sect4{EulerEquations::compute_normal_flux} + // On the boundaries of the // domain and across hanging // nodes we use a numerical flux @@ -345,12 +353,12 @@ struct EulerEquations // $\alpha$. It's form has also // been given already in the // introduction: - template + template static void numerical_normal_flux(const Point &normal, - const std::vector &Wplus, - const std::vector &Wminus, - const double alpha, + const InputVector &Wplus, + const InputVector &Wminus, + const double alpha, Sacado::Fad::DFad (&normal_flux)[n_components]) { Sacado::Fad::DFad iflux[n_components][dim]; @@ -363,12 +371,13 @@ struct EulerEquations { normal_flux[di] = 0; for (unsigned int d=0; dWminus will of course be + // modified, so it shouldn't be a + // const argument. Yet it is + // in the implementation below, and needs + // to be in order to allow the code to + // compile. The reason is that we call + // this function at a place where + // Wminus is of type + // Table@<2,Sacado::Fad::DFad@ + // @>, this being 2d table with + // indices representing the quadrature + // point and the vector component, + // respectively. We call this function + // with Wminus[q] as last + // argument; subscripting a 2d table + // yields a temporary accessor object + // representing a 1d vector, just what we + // want here. The problem is that a + // temporary accessor object can't be + // bound to a non-const reference + // argument of a function, as we would + // like here, according to the C++ 1998 + // and 2003 standards (something that + // will be fixed with the next standard + // in the form of rvalue references). We + // get away with making the output + // argument here a constant because it is + // the accessor object that's + // constant, not the table it points to: + // that one can still be written to. The + // hack is unpleasant nevertheless + // because it restricts the kind of data + // types that may be used as template + // argument to this function: a regular + // vector isn't going to do because that + // one can not be written to when marked + // const. With no good + // solution around at the moment, we'll + // go with the pragmatic, even if not + // pretty, solution shown here: + template + static + void + compute_Wminus (const BoundaryKind (&boundary_kind)[n_components], + const Point &normal_vector, + const DataVector &Wplus, + const Vector &boundary_values, + const DataVector &Wminus) + { + for (unsigned int c = 0; c < n_components; c++) + switch (boundary_kind[c]) + { + case inflow_boundary: + { + Wminus[c] = boundary_values(c); + break; + } + + case outflow_boundary: + { + Wminus[c] = Wplus[c]; + break; + } + + // Prescribed pressure boundary + // conditions are a bit more + // complicated by the fact that + // even though the pressure is + // prescribed, we really are + // setting the energy component + // here, which will depend on + // velocity and pressure. So + // even though this seems like + // a Dirichlet type boundary + // condition, we get + // sensitivities of energy to + // velocity and density (unless + // these are also prescribed): + case pressure_boundary: + { + const typename DataVector::value_type + density = (boundary_kind[density_component] == + inflow_boundary + ? + boundary_values(density_component) + : + Wplus[density_component]); + + typename DataVector::value_type kinetic_energy = 0; + for (unsigned int d=0; d vdotn = 0; + for (unsigned int d = 0; d < dim; d++) { + vdotn += Wplus[d]*normal_vector[d]; + } + + Wminus[c] = Wplus[c] - 2.0*vdotn*normal_vector[c]; + break; + } + + default: + Assert (false, ExcNotImplemented()); + } + } + + + // @sect4{EulerEquations::Postprocessor} // Finally, we declare a class that // implements a postprocessing of data @@ -1149,17 +1331,11 @@ namespace Parameters { static const unsigned int max_n_boundaries = 10; - enum BoundaryKind - { - inflow_boundary, - outflow_boundary, - no_penetration_boundary, - pressure_boundary - }; - struct BoundaryConditions { - BoundaryKind kind[EulerEquations::n_components]; + typename EulerEquations::BoundaryKind + kind[EulerEquations::n_components]; + FunctionParser values; BoundaryConditions (); @@ -1169,7 +1345,6 @@ namespace Parameters AllParameters (); double diffusion_power; - double gravity; double time_step, final_time; double theta; @@ -1212,10 +1387,6 @@ namespace Parameters Patterns::Double(), "power of mesh size for diffusion"); - prm.declare_entry("gravity", "0.0", - Patterns::Double(), - "gravity forcing"); - prm.enter_subsection("time stepping"); { prm.declare_entry("time step", "0.1", @@ -1282,7 +1453,6 @@ namespace Parameters { mesh_filename = prm.get("mesh"); diffusion_power = prm.get_double("diffusion power"); - gravity = prm.get_double("gravity"); prm.enter_subsection("time stepping"); { @@ -1318,13 +1488,17 @@ namespace Parameters = prm.get("w_" + Utilities::int_to_string(di)); if ((di < dim) && (no_penetration == true)) - boundary_conditions[boundary_id].kind[di] = no_penetration_boundary; + boundary_conditions[boundary_id].kind[di] + = EulerEquations::no_penetration_boundary; else if (boundary_type == "inflow") - boundary_conditions[boundary_id].kind[di] = inflow_boundary; + boundary_conditions[boundary_id].kind[di] + = EulerEquations::inflow_boundary; else if (boundary_type == "pressure") - boundary_conditions[boundary_id].kind[di] = pressure_boundary; + boundary_conditions[boundary_id].kind[di] + = EulerEquations::pressure_boundary; else if (boundary_type == "outflow") - boundary_conditions[boundary_id].kind[di] = outflow_boundary; + boundary_conditions[boundary_id].kind[di] + = EulerEquations::outflow_boundary; else AssertThrow (false, ExcNotImplemented()); @@ -1397,14 +1571,14 @@ class ConservationLaw void assemble_system (); void assemble_cell_term (const FEValues &fe_v, const std::vector &dofs); - void assemble_face_term(const unsigned int face_no, - const FEFaceValuesBase &fe_v, - const FEFaceValuesBase &fe_v_neighbor, - const std::vector &dofs, - const std::vector &dofs_neighbor, - const bool external_face, - const unsigned int boundary_id, - const double face_diameter); + void assemble_face_term (const unsigned int face_no, + const FEFaceValuesBase &fe_v, + const FEFaceValuesBase &fe_v_neighbor, + const std::vector &dofs, + const std::vector &dofs_neighbor, + const bool external_face, + const unsigned int boundary_id, + const double face_diameter); std::pair solve (Vector &solution); @@ -1836,13 +2010,13 @@ void ConservationLaw::assemble_system () neighbor_child->get_dof_indices (dof_indices_neighbor); - assemble_face_term(face_no, fe_v_subface, - fe_v_face_neighbor, - dof_indices, - dof_indices_neighbor, - false, - numbers::invalid_unsigned_int, - neighbor_child->diameter()); + assemble_face_term (face_no, fe_v_subface, + fe_v_face_neighbor, + dof_indices, + dof_indices_neighbor, + false, + numbers::invalid_unsigned_int, + neighbor_child->diameter()); } } @@ -1882,13 +2056,13 @@ void ConservationLaw::assemble_system () neighbor_face_no, neighbor_subface_no); - assemble_face_term(face_no, fe_v_face, - fe_v_subface_neighbor, - dof_indices, - dof_indices_neighbor, - false, - numbers::invalid_unsigned_int, - cell->face(face_no)->diameter()); + assemble_face_term (face_no, fe_v_face, + fe_v_subface_neighbor, + dof_indices, + dof_indices_neighbor, + false, + numbers::invalid_unsigned_int, + cell->face(face_no)->diameter()); } } } @@ -2264,12 +2438,15 @@ assemble_cell_term (const FEValues &fe_v, // @sect4{ConservationLaw::assemble_face_term} // - // These are either - // boundary terms or terms across differing - // levels of refinement. In the first case, - // fe_v==fe_v_neighbor and dofs==dofs_neighbor. - // The int boundary < 0 if not at a boundary, - // otherwise it is the boundary indicator. + // Here, we do essentially the same as in the + // previous function. t the top, we introduce + // the independent variables. Because the + // current function is also used if we are + // working on an internal face between two + // cells, the independent variables are not + // only the degrees of freedom on the current + // cell but in the case of an interior face + // also the ones on the neighbor. template void ConservationLaw::assemble_face_term(const unsigned int face_no, @@ -2282,160 +2459,126 @@ ConservationLaw::assemble_face_term(const unsigned int face_no, const double face_diameter) { const unsigned int n_q_points = fe_v.n_quadrature_points; - const unsigned int dofs_per_cell = fe_v.get_fe().dofs_per_cell; - const unsigned int ndofs_per_cell = fe_v_neighbor.get_fe().dofs_per_cell; - Assert(dofs_per_cell == ndofs_per_cell, - ExcDimensionMismatch(dofs_per_cell, ndofs_per_cell)); - - // As above, the fad degrees of freedom - std::vector > independent_local_dof_values(dofs_per_cell+ndofs_per_cell); - - // The conservative variables for this cell, - // and for - std::vector > > Wplus (n_q_points, - std::vector >(EulerEquations::n_components)); - std::vector > > Wminus (n_q_points, - std::vector >(EulerEquations::n_components)); - - - const std::vector > &normals = fe_v.get_normal_vectors (); - - - // If we are at a boundary, then - // dofs_neighbor are the same as dofs, so - // we do not want to duplicate them. If - // there is a neighbor cell, then we want - // to include them. - int ndofs = (external_face == false ? dofs_per_cell + ndofs_per_cell : dofs_per_cell); - // Set the local independent_local_dof_valuesS. - for (unsigned int in = 0; in < dofs_per_cell; in++) { - independent_local_dof_values[in] = current_solution(dof_indices[in]); - independent_local_dof_values[in].diff(in, ndofs); - } - // If present, set the neighbor dofs. - if (external_face == false) - for (unsigned int in = 0; in < ndofs_per_cell; in++) { - independent_local_dof_values[in+dofs_per_cell] = current_solution(dof_indices_neighbor[in]); - independent_local_dof_values[in+dofs_per_cell].diff(in+dofs_per_cell, ndofs); - } + const unsigned int dofs_per_cell = fe_v.dofs_per_cell; - // Set the values of the local conservative variables. - // Initialize all variables to zero. - for (unsigned int q = 0; q < n_q_points; q++) { - for (unsigned int di = 0; di < EulerEquations::n_components; di++) { - Wplus[q][di] = 0; - Wminus[q][di] = 0; - } - for (unsigned int sf = 0; sf < dofs_per_cell; sf++) { - int di = fe_v.get_fe().system_to_component_index(sf).first; - Wplus[q][di] += - (parameters.theta*independent_local_dof_values[sf]+(1.0-parameters.theta)*old_solution(dof_indices[sf]))*fe_v.shape_value_component(sf, q, di); + std::vector > + independent_local_dof_values (dofs_per_cell), + independent_neighbor_dof_values (external_face == false ? + dofs_per_cell : + 0); + + 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]); + independent_local_dof_values[i].diff(i, n_independent_variables); } + if (external_face == false) + for (unsigned int i = 0; i < dofs_per_cell; i++) + { + independent_neighbor_dof_values[i] + = current_solution(dof_indices_neighbor[i]); + independent_neighbor_dof_values[i] + .diff(i+dofs_per_cell, n_independent_variables); + } + + + // Next, we need to define the values of + // the conservative variables $\tilde + // {\mathbf W}$ on this side of the face + // ($\tilde {\mathbf W}^+$) and on the + // opposite side ($\tilde {\mathbf + // W}^-$). The former can be computed in + // exactly the same way as in the previous + // function, but note that the + // fe_v variable now is of + // type FEFaceValues or FESubfaceValues: + Table<2,Sacado::Fad::DFad > + Wplus (n_q_points, EulerEquations::n_components), + Wminus (n_q_points, EulerEquations::n_components); - // If there is a cell across, then initialize - // the exterior trace as a function of the other - // cell degrees of freedom. - if (external_face == false) { - for (unsigned int sf = 0; sf < ndofs_per_cell; sf++) { - int di = fe_v_neighbor.get_fe().system_to_component_index(sf).first; - Wminus[q][di] += - (parameters.theta*independent_local_dof_values[sf+dofs_per_cell]+(1.0-parameters.theta)*old_solution(dof_indices_neighbor[sf]))* - fe_v_neighbor.shape_value_component(sf, q, di); + for (unsigned int q=0; q::max_n_boundaries, ExcIndexRange (boundary_id, 0, Parameters::AllParameters::max_n_boundaries)); - // Evaluate the function object. This is - // a bit tricky; a given boundary might - // have both prescribed and implicit - // values. If a particular component is - // not prescribed, the values evaluate to - // zero and are ignored, below. - std::vector > bvals(n_q_points, Vector(EulerEquations::n_components)); - parameters.boundary_conditions[boundary_id].values.vector_value_list(fe_v.get_quadrature_points(), bvals); - - // We loop the quadrature points, and we treat each - // component individualy. - for (unsigned int q = 0; q < n_q_points; q++) { - for (unsigned int di = 0; di < EulerEquations::n_components; di++) { - - // An inflow/dirichlet type of boundary condition - if (parameters.boundary_conditions[boundary_id].kind[di] == Parameters::AllParameters::inflow_boundary) { - Wminus[q][di] = bvals[q](di); - } else if (parameters.boundary_conditions[boundary_id].kind[di] == Parameters::AllParameters::pressure_boundary) { - // A prescribed pressure boundary - // condition. This boundary - // condition is complicated by the - // fact that even though the - // pressure is prescribed, we - // really are setting the energy - // index here, which will depend on - // velocity and pressure. So even - // though this seems like a - // dirichlet type boundary - // condition, we get sensitivities - // of energy to velocity and - // density (unless these are also - // prescribed. - Sacado::Fad::DFad rho_vel_sqr = 0; - Sacado::Fad::DFad dens; - - dens = parameters.boundary_conditions[boundary_id].kind[EulerEquations::density_component] == Parameters::AllParameters::inflow_boundary ? bvals[q](EulerEquations::density_component) : - Wplus[q][EulerEquations::density_component]; - - for (unsigned int d=0; d < dim; d++) { - if (parameters.boundary_conditions[boundary_id].kind[d] == Parameters::AllParameters::inflow_boundary) - rho_vel_sqr += bvals[q](d)*bvals[q](d); - else - rho_vel_sqr += Wplus[q][d]*Wplus[q][d]; - } - rho_vel_sqr /= dens; - // Finally set the energy value as determined by the - // prescribed pressure and the other variables. - Wminus[q][di] = bvals[q](di)/(EulerEquations::gas_gamma-1.0) + - 0.5*rho_vel_sqr; - - } else if (parameters.boundary_conditions[boundary_id].kind[di] == Parameters::AllParameters::outflow_boundary) { - // A free/outflow boundary, very simple. - Wminus[q][di] = Wplus[q][di]; - - } else { - // We must be at a no-penetration - // boundary. We prescribe the - // velocity (we are dealing with a - // particular component here so - // that the average of the - // velocities is orthogonal to the - // surface normal. This creates - // sensitivies of across the - // velocity components. - Sacado::Fad::DFad vdotn = 0; - for (unsigned int d = 0; d < dim; d++) { - vdotn += Wplus[q][d]*normals[q](d); - } - - Wminus[q][di] = Wplus[q][di] - 2.0*vdotn*normals[q](di); - } - } - } // for q - } // b>= 0 - + std::vector > + boundary_values(n_q_points, Vector(EulerEquations::n_components)); + parameters.boundary_conditions[boundary_id] + .values.vector_value_list(fe_v.get_quadrature_points(), + boundary_values); + + for (unsigned int q = 0; q < n_q_points; q++) + EulerEquations::compute_Wminus (parameters.boundary_conditions[boundary_id].kind, + fe_v.normal_vector(q), + Wplus[q], + boundary_values[q], + Wminus[q]); + } + // Determine the Lax-Friedrich's stability parameter, // and evaluate the numerical flux function at the quadrature points typedef Sacado::Fad::DFad NormalFlux[EulerEquations::n_components]; @@ -2457,7 +2600,8 @@ ConservationLaw::assemble_face_term(const unsigned int face_no, } for (unsigned int q=0; q::numerical_normal_flux(normals[q], Wplus[q], Wminus[q], alpha, + EulerEquations::numerical_normal_flux(fe_v.normal_vector(q), + Wplus[q], Wminus[q], alpha, normal_fluxes[q]); // Now assemble the face term -- 2.39.5