From ca3dd85986957856d173dbab783ada308893648e Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 16 May 2008 03:41:39 +0000 Subject: [PATCH] Implement the postprocessor. Will hook it in tomorrow. git-svn-id: https://svn.dealii.org/trunk@16108 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-33/step-33.cc | 375 +++++++++++++++++++++++++--- 1 file changed, 343 insertions(+), 32 deletions(-) diff --git a/deal.II/examples/step-33/step-33.cc b/deal.II/examples/step-33/step-33.cc index 4a2edfc2e8..99be279ecd 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -158,6 +158,67 @@ 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 + // variables. This we can do based on the + // energy density and the kinetic energy + // $\frac 12 \rho |\mathbf v|^2 = + // \frac{|\rho \mathbf v|^2}{2\rho}$ + // (note that the independent variables + // contain the momentum components $\rho + // v_i$, not the velocities $v_i$). + // + // There is one slight problem: We will + // need to call the following functions + // with input arguments of type + // std::vector@ and + // Vector@. The + // problem is that the former has an + // access operator + // operator[] whereas the + // latter, for historical reasons, has + // operator(). We wouldn't + // be able to write the function in a + // generic way if we were to use one or + // the other of these. Fortunately, we + // can use the following trick: instead + // of writing v[i] or + // v(i), we can use + // *(v.begin() + i), i.e. we + // generate an iterator that points to + // the ith element, and then + // dereference it. This works for both + // kinds of vectors -- not the prettiest + // solution, but one that works. + template + static + number + compute_kinetic_energy (const InputVector &W) + { + number kinetic_energy = 0; + for (unsigned int d=0; d + static + number + compute_pressure (const InputVector &W) + { + return ((gas_gamma-1.0) * + (*(W.begin() + energy_component) - + compute_kinetic_energy(W))); + } + + + // We define the flux function // $F(W)$ as one large matrix. // Each row of this matrix @@ -197,37 +258,21 @@ struct EulerEquations void flux_matrix (const std::vector &W, number (&flux)[n_components][dim]) { - // First compute the pressure - // that appears in the flux - // matrix, based on the - // energy density and the - // kinetic energy $\frac 12 - // \rho |\mathbf v|^2 = - // \frac{|\rho \mathbf - // v|^2}{2\rho}$ (note that - // the independent variables - // contain the momentum - // components $\rho v_i$, not - // the velocities $v_i$): - number kinetic_energy = 0; - for (unsigned int d=0; ddim columns of the + // matrix that correspond to the + // momentum terms: + const number pressure = compute_pressure (W); - const number pressure = (gas_gamma-1.0)*(W[energy_component] - kinetic_energy); - - // Then compute the first - // dim columns - // of the matrix that - // correspond to the momentum - // terms: for (unsigned int d=0; d + { + public: + Postprocessor (const bool do_schlieren_plot); + + virtual + void + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &duh, + const std::vector > > &dduh, + const std::vector > &normals, + 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; + + private: + const bool do_schlieren_plot; + }; }; template const double EulerEquations::gas_gamma = 1.4; + + +template +EulerEquations::Postprocessor:: +Postprocessor (const bool do_schlieren_plot) + : + do_schlieren_plot (do_schlieren_plot) +{} + + + // This is the only function worth commenting + // on. When generating graphical output, the + // DataOut and related classes will call this + // function on each cell, with values, + // gradients, hessians, and normal vectors + // (in case we're working on faces) at each + // quadrature point. Note that the data at + // each quadrature point is itself + // vector-valued, namely the conserved + // variables. What we're going to do here is + // to compute the quantities we're interested + // in at each quadrature point. Note that for + // this we can ignore the hessians ("dduh") + // and normal vectors; to avoid compiler + // warnings about unused variables, we + // comment out their names. +template +void +EulerEquations::Postprocessor:: +compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &duh, + const std::vector > > &/*dduh*/, + const std::vector > &/*normals*/, + std::vector > &computed_quantities) const +{ + // At the beginning of the function, let us + // make sure that all variables have the + // correct sizes, so that we can access + // individual vector elements without + // having to wonder whether we might read + // or write invalid elements; we also check + // that the duh vector only + // contains data if we really need it (the + // system knows about this because we say + // so in the + // get_needed_update_flags() + // function below). For the inner vectors, + // we check that at least the first element + // of the outer vector has the correct + // inner size: + const unsigned int n_quadrature_points = uh.size(); + + if (do_schlieren_plot == true) + Assert (duh.size() == n_quadrature_points, + ExcInternalError()) + else + Assert (duh.size() == 0, + ExcInternalError()); + + Assert (computed_quantities.size() == n_quadrature_points, + ExcInternalError()); + + Assert (uh[0].size() == n_components, + ExcInternalError()); + + if (do_schlieren_plot == true) + Assert (computed_quantities[0].size() == dim+2, + ExcInternalError()) + else + Assert (computed_quantities[0].size() == dim+1, + ExcInternalError()); + + // Then loop over all quadrature points and + // do our work there. The code should be + // pretty self-explanatory. The order of + // output variables is first + // dim velocities, then the + // pressure, and if so desired the + // schlieren plot. Note that we try to be + // generic about the order of variables in + // the input vector, using the + // first_momentum_component + // and density_component + // information: + for (unsigned int q=0; q (uh[q]); + + if (do_schlieren_plot == true) + computed_quantities[q](dim+1) = duh[q][density_component] * + duh[q][density_component]; + } +} + + +template +std::vector +EulerEquations::Postprocessor:: +get_names () const +{ + std::vector names; + for (unsigned int d=0; d +std::vector +EulerEquations::Postprocessor:: +get_data_component_interpretation () const +{ + std::vector + interpretation; + + for (unsigned int d=0; d +UpdateFlags +EulerEquations::Postprocessor:: +get_needed_update_flags () const +{ + if (do_schlieren_plot == true) + return update_values | update_gradients; + else + return update_values; +} + + + +template +unsigned int +EulerEquations::Postprocessor:: +n_output_variables () const +{ + if (do_schlieren_plot == true) + return dim+2; + else + return dim+1; +} + + +template class EulerEquations<2>::Postprocessor; + + + // @sect3{Run time parameter handling} // Our next job is to define a few @@ -588,8 +897,10 @@ namespace Parameters const std::string stab = prm.get("stab"); 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"); } @@ -1633,10 +1944,10 @@ void ConsLaw::solve (Vector &dsolution, int &niter, double &lin_res } } - // @sect3{Postprocessing and Output} - // Recover the physical variables from the conservative - // variables so that output will be (perhaps) more - // meaningfull. + // @sect3{Postprocessing and Output} Recover + // the physical variables from the + // conservative variables so that output will + // be (perhaps) more meaningfull. template void ConsLaw::postprocess() { const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; -- 2.39.5