-/* step-12.cc,v 1.33 2005/08/08 16:41:40 wolf Exp */\r
-/* Author: David Neckels, Boulder Colorado 2007 */\r
-/* step-12.cc,v 1.33 2005/08/08 16:41:40 wolf Exp */\r
-/* Version: Version-5-2-0 */\r
-/* */\r
-/* Copyright (C) 2001, 2002, 2003, 2004, 2005 by the deal.II authors */\r
-/* */\r
-/* This file is subject to QPL and may not be distributed */\r
-/* without copyright and license information. Please refer */\r
-/* to the file deal.II/doc/license.html for the text and */\r
-/* further information on this license. */\r
-\r
- // This program solves the Euler equations\r
- // of gas dynamics for a given configuration\r
- // file. It uses a standard Galerkin approach\r
- // with weakly applied boundary conditions.\r
- \r
- // @sect3{Include files}\r
-\r
- // Aztecoo require mpi (even though we run on only\r
- // one processor in this example).\r
-#include <mpi.h>\r
-\r
- // Here we have the necessary TRILINOS includes.\r
- //\r
- // Epetra is the basic trilinos vector/matrix library.\r
-#include <Epetra_MpiComm.h>\r
-#include <Epetra_Map.h>\r
-#include <Epetra_CrsGraph.h>\r
-#include <Epetra_CrsMatrix.h>\r
-#include <Epetra_Vector.h>\r
- // Teuchos is a Trilinos utility library that is used\r
- // to set parameters within the Aztec solver library.\r
-#include "Teuchos_ParameterList.hpp"\r
- // Aztec is the iterative solver library.\r
-#include <AztecOO.h>\r
-#include <AztecOO_Operator.h>\r
-#define HAVE_IFPACK_TEUCHOS\r
-#include <Ifpack.h>\r
-\r
- // Amesos is a direct solver package within Trilinos.\r
-#include <Amesos.h>\r
- // Sacado is the automatic differentiation package, which\r
- // is used to find the jacobian for a fully implicit Newton\r
- // iteration.\r
-#include <Sacado.hpp>\r
-\r
- // A standard set of dealii includes. Nothing special to\r
- // comment on here.\r
-#include <base/quadrature_lib.h>\r
-#include <base/function.h>\r
-#include <base/parameter_handler.h>\r
-#include <base/function_parser.h>\r
-\r
-#include <lac/vector.h>\r
-#include <lac/sparse_matrix.h>\r
-#include <lac/vector_memory.h>\r
-\r
-#include <grid/tria.h>\r
-#include <grid/grid_generator.h>\r
-#include <grid/grid_out.h>\r
-#include <grid/grid_refinement.h>\r
-#include <grid/tria_accessor.h>\r
-#include <grid/tria_iterator.h>\r
-#include <grid/grid_in.h>\r
-\r
-#include <fe/fe_values.h>\r
-#include <fe/fe_system.h>\r
-\r
-#include <dofs/dof_handler.h>\r
-#include <dofs/dof_accessor.h>\r
-#include <dofs/dof_tools.h>\r
-\r
-#include <numerics/data_out.h>\r
-#include <numerics/vectors.h>\r
-#include <numerics/solution_transfer.h>\r
-\r
-#include <fe/mapping_q1.h>\r
-#include <fe/fe_q.h>\r
-#include <numerics/derivative_approximation.h>\r
- // And this again is C++:\r
-#include <iostream>\r
-#include <fstream>\r
-#include <vector>\r
-\r
- // Introduce the dealii library into the current namespace.\r
-using namespace dealii;\r
-\r
- // We define a shorter name for the automatic differentiation\r
- // type.\r
-typedef Sacado::Fad::DFad<double> fad_double;\r
-typedef unsigned int UInt;\r
- // The Epetra library requires a 'communicator', which describes\r
- // the layout of a parallel (or serial) set of processors.\r
-Epetra_MpiComm *Comm;\r
-\r
- //@sect3{Flux function definition}\r
- // Here we define the flux function for this system of conservation\r
- // laws. Note: it would be terribly difficult to use this example\r
- // to solve some other system of conservation laws.\r
- //\r
- // We define the number of components in the system. Euler's has\r
- // one entry for momenta in each spatial direction, plus the energy\r
- // and density components.\r
-#define N_COMP (2 + DIMENSION)\r
- // Define a handle to the density and energy indices. We have arrange\r
- // the momenta to be first, then density, and, lastly, energy.\r
-#define DENS_IDX DIMENSION\r
-#define ENERGY_IDX (DIMENSION+1)\r
-\r
- // The gas constant. This value is representative of air.\r
-const double GAMMA = 1.4;\r
- // We define the flux functions as one large matrix. Each row of this\r
- // matrix represents a scalar conservation law for the component in\r
- // that row. We template the numerical type of the flux function\r
- // so that we may use the automatic differentiation type here.\r
- // The flux functions are defined in terms of the\r
- // conserved variables $\rho w_0, \dots, \rho w_{d-1}, \rho, E$,\r
- // so they do not look exactly like the Euler equations one is\r
- // used to seeing. We evaluate the flux at a single quadrature\r
- // point.\r
-template <typename number, int dim>\r
-void Flux(std::vector<std::vector<number> > &flux, \r
- const Point<dim> &point, \r
- const std::vector<number> &W)\r
-{\r
-\r
- // Pressure is a dependent variable: $p = \r
- // (\gamma - 1)(E-\frac{1}{2} \rho |v|^2)$.\r
- number rho_normVsqr;\r
- for (int d0 = 0; d0 < dim; d0++) rho_normVsqr += W[d0]*W[d0];\r
- // Since W are $\rho v$, we get a $\rho^2$ in the\r
- // numerator, so dividing a $\rho$ out gives the desired $ \rho |v|^2$.\r
- rho_normVsqr /= W[DENS_IDX];\r
-\r
- number pressure = (GAMMA-1.0)*(W[ENERGY_IDX] - number(0.5)*(rho_normVsqr));\r
-\r
- // We compute the momentum terms. We divide by the\r
- // density here to get $v_i \rho v_j$\r
- for (int d = 0; d < dim; d++) {\r
- for (int d1 = 0; d1 < dim; d1++) {\r
- flux[d][d1] = W[d]*W[d1]/W[DENS_IDX];\r
- }\r
- // The pressure contribution, along the diagonal:\r
- flux[d][d] += pressure;\r
- // Advection/conservation of density:\r
- flux[DENS_IDX][d] = W[d]; \r
- // And, lastly, conservation of energy.\r
- flux[ENERGY_IDX][d] = W[d]/W[DENS_IDX]*\r
- (W[ENERGY_IDX] + pressure); // energy\r
- }\r
-}\r
-\r
- // On the boundaries of the domain and across `hanging nodes` we use\r
- // a numerical flux function to enforce boundary conditions. This routine\r
- // is the basic Lax-Friedrich's flux with a stabilization parameter\r
- // $\alpha$.\r
-template <typename number, int dim>\r
-void LFNumFlux(\r
- std::vector<std::vector<fad_double> > &nflux,\r
- const std::vector<Point<dim> > &points, \r
- const std::vector<Point<dim> > &normals,\r
- const std::vector<std::vector<number> > &Wplus,\r
- const std::vector<std::vector<number> > &Wminus,\r
- double alpha)\r
-{\r
- int n_q_points = points.size();\r
-\r
- // We evaluate the flux at each of the quadrature points.\r
- for (int q = 0; q < n_q_points; q++) {\r
- std::vector<std::vector<fad_double> > iflux(N_COMP,\r
- std::vector<fad_double>(dim, 0));\r
- std::vector<std::vector<fad_double> > oflux(N_COMP,\r
- std::vector<fad_double>(dim, 0));\r
-\r
- Flux<number, dim>(iflux, points[q], Wplus[q]);\r
- Flux<number, dim>(oflux, points[q], Wminus[q]);\r
-\r
- for (int di = 0; di < N_COMP; di++) {\r
- nflux[q][di] = 0;\r
- for (int d = 0; d < dim; d++) {\r
- nflux[q][di] += 0.5*(iflux[di][d] + oflux[di][d])*normals[q](d);\r
- }\r
- nflux[q][di] += 0.5*alpha*(Wplus[q][di] - Wminus[q][di]);\r
- }\r
- }\r
-\r
-}\r
-\r
- // @sect3{Initial and side condition parsing}\r
- // For the initial condition we use the expression parser function\r
- // object.\r
-template <int dim>\r
-class InitialCondition : public FunctionParser<dim> \r
-{\r
- public:\r
- InitialCondition ();\r
- \r
- // This function should be called after parsing, but before using\r
- // the object. It formalizes the expressions and initializes the\r
- // function parser with the appropriate expressions.\r
- void Init();\r
-\r
- // During parsing we call this function as the initial condition\r
- // for one of the $\mathbf{w}$ variables is encountered.\r
-\r
- void set_ic(int _row, std::string &expr) {\r
- expressions[_row] = expr;\r
- }\r
-\r
- virtual void vector_value_list (const std::vector<Point<dim> > &points,\r
- std::vector<Vector<double> > &value_list) const;\r
- private:\r
- std::vector<std::string> expressions;\r
-};\r
-\r
-template <int dim>\r
-InitialCondition<dim>::InitialCondition () :\r
- FunctionParser<dim> (N_COMP),\r
- expressions(N_COMP, "0.0")\r
-{}\r
-\r
- // Here we set up x,y,z as the variables that one should use in the input\r
- // deck to describe their initial condition.\r
-template<int dim>\r
-void InitialCondition<dim>::Init() {\r
- std::map<std::string, double> constants;\r
- constants["M_PI"] = M_PI;\r
- std::string variables = (dim == 2 ? "x,y" : "x,y,z");\r
-\r
- FunctionParser<dim>::initialize(variables, expressions, constants);\r
-\r
-}\r
-\r
-template <int dim>\r
-void InitialCondition<dim>::vector_value_list (const std::vector<Point<dim> > &points,\r
- std::vector<Vector<double> > &value_list) const \r
-{\r
- const unsigned int n_points = points.size();\r
-\r
- Assert (value_list.size() == n_points, \r
- ExcDimensionMismatch (value_list.size(), n_points));\r
-\r
- for (unsigned int p=0; p<n_points; ++p)\r
- InitialCondition<dim>::vector_value (points[p],\r
- value_list[p]);\r
-}\r
-\r
- // As above, we use the expression function parser for boundary conditions.\r
-template <int dim>\r
-class SideCondition : public FunctionParser<dim> \r
-{\r
- public:\r
- SideCondition (int ncomp);\r
- ~SideCondition ();\r
-\r
- // As above.\r
- void Init();\r
- // As above.\r
- void set_coeff_row(int _row_n, std::string &expr);\r
- \r
- virtual void vector_value_list (const std::vector<Point<dim> > &points,\r
- std::vector<Vector<double> > &value_list) const;\r
- private:\r
- std::vector<std::string> expressions;\r
-};\r
-\r
-template <int dim>\r
-SideCondition<dim>::SideCondition (int ncomp) :\r
- FunctionParser<dim> (ncomp),\r
- expressions(ncomp, "0.0")\r
-{\r
-}\r
-template <int dim>\r
-void SideCondition<dim>::set_coeff_row (int _row_n, std::string &expr) \r
-{\r
- expressions[_row_n] = expr;\r
-}\r
-\r
-template <int dim>\r
-void SideCondition<dim>::Init() {\r
- std::map<std::string, double> constants;\r
- constants["M_PI"] = M_PI;\r
- std::string variables = (dim == 2 ? "x,y" : "x,y,z");\r
-\r
- FunctionParser<dim>::initialize(variables, expressions, constants);\r
-\r
-}\r
-\r
-template <int dim>\r
-SideCondition<dim>::~SideCondition () \r
-{\r
-}\r
-\r
-template <int dim>\r
-void SideCondition<dim>::vector_value_list (const std::vector<Point<dim> > &points,\r
- std::vector<Vector<double> > &value_list) const \r
-{\r
- const unsigned int n_points = points.size();\r
-\r
- Assert (value_list.size() == n_points, \r
- ExcDimensionMismatch (value_list.size(), n_points));\r
-\r
- for (unsigned int p=0; p<n_points; ++p)\r
- SideCondition<dim>::vector_value (points[p],\r
- value_list[p]);\r
-}\r
- //@sect3{Conservation Law class}\r
- // Here we define a Conservation Law class that helps group\r
- // operations and data for our Euler equations into a manageable\r
- // entity. Functions will be described as their definitions appear.\r
-template <int dim>\r
-class ConsLaw\r
-{\r
- public:\r
- ConsLaw ();\r
- ~ConsLaw ();\r
-\r
- void run ();\r
- void declare_parameters();\r
- void load_parameters(const char *);\r
- \r
- private:\r
- void build_fe();\r
- void setup_system ();\r
- void initialize_system ();\r
- void assemble_system (double &res_norm);\r
- void solve (Vector<double> &solution, int &, double &);\r
- void refine_grid ();\r
- void output_results (const unsigned int cycle) const;\r
- void initialize();\r
- void zero_matrix();\r
- void estimate();\r
- void postprocess();\r
- void compute_predictor();\r
- \r
- Triangulation<dim> triangulation;\r
- const MappingQ1<dim> mapping;\r
- \r
- \r
- FESystem<dim> *fe_ptr;\r
-\r
- DoFHandler<dim> dof_handler;\r
-\r
- SparsityPattern sparsity_pattern;\r
- const QGauss<dim> quadrature;\r
- const QGauss<dim-1> face_quadrature;\r
- \r
- // The actual solution to the Euler equation\r
- Vector<double> solution;\r
- // The current value of the solution during the Newton iteration\r
- Vector<double> nlsolution;\r
- // An estimate of the next time value; used for adaptivity and as a\r
- // guess for the next Newton iteration.\r
- Vector<double> predictor;\r
- // Values after post-processing (used to output the physical variables).\r
- Vector<double> ppsolution;\r
- // The solution to the linear problem during the Newton iteration\r
- Vector<double> dsolution;\r
- Vector<double> right_hand_side;\r
- \r
- public:\r
-\r
- void assemble_cell_term(const FEValues<dim>& fe_v,\r
- std::vector<unsigned int> &dofs,\r
- unsigned int cell_no\r
- );\r
- \r
- void assemble_face_term(\r
- int face_no,\r
- const FEFaceValuesBase<dim>& fe_v,\r
- const FEFaceValuesBase<dim>& fe_v_neighbor,\r
- std::vector<unsigned int> &dofs,\r
- std::vector<unsigned int> &dofs_neighbor,\r
- int boundary = -1\r
- );\r
-\r
- unsigned int get_n_components() const { return N_COMP;}\r
-\r
- private:\r
- // T = current time, dT = time step, TF = final time.\r
- double T, dT, TF;\r
- double face_diameter;\r
- double cell_diameter;\r
- // An object to handle parsing the input deck.\r
- ParameterHandler prm;\r
- // Name of the mesh to read in.\r
- string mesh;\r
- InitialCondition<dim> ic;\r
-\r
- // Enums for the various supported boundary conditions.\r
- typedef enum {INFLOW_BC = 1, OUTFLOW_BC=2, NO_PENETRATION_BC=3, PRESSURE_BC=4} bc_type;\r
-\r
- // For each boundary we store a map from boundary # to the type\r
- // of boundary condition. If the boundary condition is prescribed,\r
- // we store a pointer to a function object that will hold the expression\r
- // for that boundary condition.\r
- typedef typename std::map<unsigned int, std::pair<std::vector<bc_type>, Function<dim>*> > bdry_map_type;\r
- bdry_map_type bdry_map;\r
-\r
- void add_boundary(unsigned int bd, std::vector<bc_type>& flags, Function<dim> *bf);\r
-\r
- // An object to store parameter information about the Aztec solver.\r
- typedef struct {\r
- int LIN_OUTPUT;\r
- typedef enum { GMRES = 0, DIRECT = 1} solver_type;\r
- solver_type SOLVER;\r
- typedef enum { QUIET = 0, VERBOSE = 1 } output_type;\r
- output_type OUTPUT;\r
- // Linear residual tolerance.\r
- double RES;\r
- int MAX_ITERS;\r
- // We use the ILUT preconditioner. This is similar\r
- // to the ILU. FILL is the number of extra entries\r
- // to add when forming the ILU decomposition.\r
- double ILUT_FILL;\r
- // When forming the preconditioner, for certain problems\r
- // bad conditioning (or just bad luck) can cause the\r
- // preconditioner to be very poorly conditioned. Hence\r
- // it can help to add diagonal perturbations to the\r
- // original matrix and form the preconditioner for this\r
- // slightly better matrix. ATOL is an absolute perturbation\r
- // that is added to the diagonal before forming the\r
- // prec, and RTOL is a scaling factor $rtol >= 1$.\r
- double ILUT_ATOL;\r
- double ILUT_RTOL;\r
- // The ILUT will drop any values that have magnitude less\r
- // than this value. This is a way to\r
- // manage the amount of memory used by this preconditioner.\r
- double ILUT_DROP;\r
- } solver_params_type;\r
-\r
- solver_params_type solver_params;\r
-\r
- // Some refinement parameters.\r
- typedef struct {\r
- typedef enum { NONE = 0, FIXED_NUMBER = 1, SHOCK = 2} refine_type;\r
- double high_frac;\r
- double low_frac;\r
- refine_type refine;\r
- double high_frac_sav;\r
- double max_cells;\r
- double shock_val;\r
- double shock_levels;\r
- } refinement_params_type;\r
-\r
- refinement_params_type refinement_params;\r
-\r
- // The user can set the stabilization parameter $\alpha$ \r
- // in the Lax-Friedrich's flux. \r
- typedef struct {\r
- typedef enum {CONSTANT=1,MESH=2} LF_stab_type;\r
- LF_stab_type LF_stab;\r
- double LF_stab_value;\r
- } flux_params_type;\r
-\r
- flux_params_type flux_params;\r
-\r
- bool is_stationary;\r
-\r
- // Power for the mesh stabilization term.\r
- double diffusion_power;\r
- double gravity;\r
- // If true, we output the squared gradient of the\r
- // density instead of density. Using this one can\r
- // create shock plots.\r
- bool schlieren_plot;\r
- // How often to create an output file.\r
- double output_step;\r
-\r
- Epetra_CrsMatrix *Matrix;\r
- Epetra_Map *Map;\r
- Vector<double> indicator;\r
- \r
- // Crank-Nicolson value\r
- const double theta; \r
-\r
-};\r
-\r
-\r
- // Asign a row of the conservation law a specified\r
- // boundary type and (possibly) function.\r
-template <int dim>\r
-void ConsLaw<dim>::add_boundary(unsigned int bd,\r
- std::vector<bc_type> &flags, Function<dim> *bf) {\r
-\r
- std::pair<std::vector<bc_type>, Function<dim> *> entry(flags, bf);\r
- bdry_map[bd] = entry;\r
-}\r
-\r
-\r
- // Apply the initialial condition. Simultaneously\r
- // initialize the non-linear solution.\r
-template <int dim>\r
-void ConsLaw<dim>::initialize() {\r
- VectorTools::interpolate(dof_handler,\r
- ic, solution);\r
- VectorTools::interpolate(dof_handler,\r
- ic, nlsolution);\r
-\r
-}\r
-\r
- // @sect3{Assembly}\r
- // @sect4{Function: assemble_cell_term}\r
- //\r
- // Assembles the cell term, adding minus the residual\r
- // to the right hand side, and adding in the Jacobian\r
- // contributions.\r
-template <int dim>\r
-void ConsLaw<dim>::assemble_cell_term(\r
- const FEValues<dim> &fe_v,\r
- std::vector<unsigned int> &dofs,\r
- unsigned int cell_no\r
- ) \r
-{\r
- // The residual for each row (i) will be accumulating \r
- // into this fad variable. At the end of the assembly\r
- // for this row, we will query for the sensitivities\r
- // to this variable and add them into the Jacobian.\r
- fad_double F_i;\r
- int dofs_per_cell = fe_v.dofs_per_cell;\r
- int n_q_points = fe_v.n_quadrature_points;\r
-\r
- // We will define the dofs on this cell in these fad variables.\r
- std::vector<fad_double> DOF(dofs_per_cell);\r
-\r
- // Values of the conservative variables at the quadrature points.\r
- std::vector<std::vector<fad_double > > W (n_q_points,\r
- std::vector<fad_double >(get_n_components()));\r
-\r
- // Values at the last time step of the conservative variables.\r
- // Note that these do not use fad variables, since they do\r
- // not depend on the 'variables to be sought'=DOFS.\r
- std::vector<std::vector<double > > Wl (n_q_points,\r
- std::vector<double >(get_n_components()));\r
-\r
- // Here we will hold the averaged values of the conservative\r
- // variables that we will linearize around (cn=Crank Nicholson).\r
- std::vector<std::vector<fad_double > > Wcn (n_q_points,\r
- std::vector<fad_double >(get_n_components()));\r
-\r
- // Gradients of the current variables. It is a\r
- // bit of a shame that we have to compute these; we almost don't.\r
- // The nice thing about a simple conservation law is that the\r
- // the flux doesn't generally involve any gradients. We do\r
- // need these, however, for the diffusion stabilization. \r
- std::vector<std::vector<std::vector<fad_double> > > Wgrads (n_q_points,\r
- std::vector<std::vector<fad_double> >(get_n_components(),\r
- std::vector<fad_double>(dim)));\r
-\r
-\r
- const std::vector<double> &JxW = fe_v.get_JxW_values ();\r
-\r
- \r
- // Here is the magical point where we declare a subset\r
- // of the fad variables as degrees of freedom. All \r
- // calculations that reference these variables (either\r
- // directly or indirectly) will accumulate sensitivies\r
- // with respect to these dofs.\r
- for (int in = 0; in < dofs_per_cell; in++) {\r
- DOF[in] = nlsolution(dofs[in]);\r
- DOF[in].diff(in, dofs_per_cell);\r
- }\r
-\r
- // Here we compute the shape function values and gradients\r
- // at the quadrature points. Ideally, we could call into \r
- // something like get_function_values, get_function_grads,\r
- // but since we don't want to make the entire solution vector\r
- // fad types, only the local cell variables, we explicitly\r
- // code this loop;\r
- for (int q = 0; q < n_q_points; q++) {\r
- for (int di = 0; di < get_n_components(); di++) {\r
- W[q][di] = 0;\r
- Wl[q][di] = 0;\r
- Wcn[q][di] = 0;\r
- for (int d = 0; d < dim; d++) {\r
- Wgrads[q][di][d] = 0;\r
- }\r
- }\r
- for (int sf = 0; sf < dofs_per_cell; sf++) {\r
- int di = fe_v.get_fe().system_to_component_index(sf).first;\r
- W[q][di] +=\r
- DOF[sf]*fe_v.shape_value_component(sf, q, di);\r
- Wl[q][di] +=\r
- solution(dofs[sf])*fe_v.shape_value_component(sf, q, di);\r
- Wcn[q][di] +=\r
- (theta*DOF[sf]+(1-theta)*solution(dofs[sf]))*fe_v.shape_value_component(sf, q, di);\r
-\r
- for (int d = 0; d < dim; d++) {\r
- Wgrads[q][di][d] += DOF[sf]*\r
- fe_v.shape_grad_component(sf, q, di)[d];\r
- } // for d\r
-\r
- }\r
-\r
- } // for q\r
-\r
- // Gather the flux values for all components at\r
- // all of the quadrature points. This also\r
- // computes the matrix of sensitivities. Perhaps\r
- // this could be done in a better way, since this\r
- // could be a rather large object, but for now it \r
- // seems to work just fine.\r
- std::vector<std::vector<std::vector<fad_double> > > flux(n_q_points, \r
- std::vector<std::vector<fad_double> >(get_n_components(),\r
- std::vector<fad_double>(dim, 0)));\r
-\r
- for (unsigned int q=0; q < n_q_points; ++q) {\r
- Flux<fad_double, dim>(flux[q], fe_v.get_quadrature_points()[q], Wcn[q]);\r
- }\r
-\r
- // We now have all of the function values/grads/fluxes,\r
- // so perform the assembly. We have an outer loop\r
- // through the components of the system, and an\r
- // inner loop over the quadrature points, where we\r
- // accumulate contributions to the ith residual.\r
- //\r
- // We initialy sum all contributions of the residual\r
- // in the positive sense, so that we don't need to\r
- // negative the Jacobian entries. Then, when we sum\r
- // into the <code> right_hand_side </code> vector,\r
- // we negate this residual.\r
- for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) \r
- {\r
- // Find which component this dof contributes to.\r
- const unsigned int\r
- component_i = fe_v.get_fe().system_to_component_index(i).first;\r
-\r
- // Initialize the fad residual to zero (removes\r
- // any previous sensitivities.\r
- F_i = 0;\r
-\r
- // Loop quadrature points.\r
- for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point) {\r
-\r
- fad_double fdotgv = 0;\r
-\r
- // Integrate the flux times gradient of the test function\r
- for (int d = 0; d < dim; d++) \r
- fdotgv += flux[point][component_i][d]*fe_v.shape_grad_component(i, point, component_i)[d];\r
- \r
- F_i -= fdotgv*JxW[point];\r
-\r
- // The mass term (if the simulation is non-stationary).\r
- fad_double delta_t= 1.0/dT*(W[point][component_i] - Wl[point][component_i]);\r
- if (!is_stationary) F_i += delta_t*\r
- fe_v.shape_value_component(i, point, component_i)*JxW[point];\r
-\r
- // Stabilization (cell wise diffusion)\r
- fad_double guv = 0;\r
- for (int d = 0; d < dim; d++) {\r
- guv += fe_v.shape_grad_component(i, point, component_i)[d]*\r
- Wgrads[point][component_i][d];\r
- }\r
-\r
- F_i += 1.0*std::pow(cell_diameter, diffusion_power)*guv*JxW[point];\r
- \r
- // The gravity component only enters into the energy \r
- // equation and into the vertical component of the \r
- // velocity.\r
- if (component_i == dim - 1) {\r
- F_i += gravity*Wcn[point][DENS_IDX]*fe_v.shape_value_component(i,point, component_i)*JxW[point];\r
- } else if (component_i == ENERGY_IDX) {\r
- F_i += gravity*Wcn[point][DENS_IDX]*Wcn[point][dim-1]*\r
- fe_v.shape_value_component(i,point, component_i)*JxW[point];\r
- }\r
- } // for q\r
-\r
- // Here we gain access to the array of sensitivities\r
- // of the residual. We then sum these into the\r
- // Epetra matrix.\r
- double *values = &(F_i.fastAccessDx(0));\r
- Matrix->SumIntoGlobalValues(dofs[i],\r
- dofs_per_cell, &values[0], reinterpret_cast<int*>(&dofs[0]));\r
- \r
- // Add minus the residual to the right hand side.\r
- right_hand_side(dofs[i]) -= F_i.val();\r
-\r
- } // for i\r
-}\r
- // @sect4{Function: assemble_face_term}\r
- // These are either\r
- // boundary terms or terms across differing \r
- // levels of refinement. In the first case,\r
- // fe_v==fe_v_neighbor and dofs==dofs_neighbor.\r
- // The int boundary < 0 if not at a boundary,\r
- // otherwise it is the boundary indicator.\r
-template <int dim>\r
-void ConsLaw<dim>::assemble_face_term(\r
- int face_no,\r
- const FEFaceValuesBase<dim>& fe_v,\r
- const FEFaceValuesBase<dim>& fe_v_neighbor, \r
- std::vector<unsigned int> &dofs,\r
- std::vector<unsigned int> &dofs_neighbor,\r
- int boundary\r
- ) \r
-{\r
- fad_double F_i;\r
- const unsigned int n_q_points = fe_v.n_quadrature_points;\r
- const unsigned int dofs_per_cell = fe_v.get_fe().dofs_per_cell;\r
- const unsigned int ndofs_per_cell = fe_v_neighbor.get_fe().dofs_per_cell;\r
- Assert(dofs_per_cell == ndofs_per_cell,\r
- ExcDimensionMismatch(dofs_per_cell, ndofs_per_cell));\r
-\r
- // As above, the fad degrees of freedom\r
- std::vector<fad_double> DOF(dofs_per_cell+ndofs_per_cell);\r
-\r
- // The conservative variables for this cell,\r
- // and for \r
- std::vector<std::vector<fad_double > > Wplus (n_q_points,\r
- std::vector<fad_double >(get_n_components()));\r
- std::vector<std::vector<fad_double > > Wminus (n_q_points,\r
- std::vector<fad_double >(get_n_components()));\r
-\r
-\r
- const std::vector<double> &JxW = fe_v.get_JxW_values ();\r
- const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();\r
-\r
-\r
- // If we are at a boundary, then dofs_neighbor are\r
- // the same as dofs, so we do not want to duplicate them.\r
- // If there is a neighbor cell, then we want to include \r
- // them.\r
- int ndofs = (boundary < 0 ? dofs_per_cell + ndofs_per_cell : dofs_per_cell);\r
- // Set the local DOFS.\r
- for (int in = 0; in < dofs_per_cell; in++) {\r
- DOF[in] = nlsolution(dofs[in]);\r
- DOF[in].diff(in, ndofs);\r
- }\r
- // If present, set the neighbor dofs.\r
- if (boundary < 0)\r
- for (int in = 0; in < ndofs_per_cell; in++) {\r
- DOF[in+dofs_per_cell] = nlsolution(dofs_neighbor[in]);\r
- DOF[in+dofs_per_cell].diff(in+dofs_per_cell, ndofs);\r
- }\r
-\r
- // Set the values of the local conservative variables.\r
- // Initialize all variables to zero.\r
- for (int q = 0; q < n_q_points; q++) {\r
- for (int di = 0; di < get_n_components(); di++) {\r
- Wplus[q][di] = 0;\r
- Wminus[q][di] = 0;\r
- }\r
- for (int sf = 0; sf < dofs_per_cell; sf++) {\r
- int di = fe_v.get_fe().system_to_component_index(sf).first;\r
- Wplus[q][di] +=\r
- (theta*DOF[sf]+(1.0-theta)*solution(dofs[sf]))*fe_v.shape_value_component(sf, q, di);\r
- }\r
-\r
-\r
- // If there is a cell across, then initialize\r
- // the exterior trace as a function of the other\r
- // cell degrees of freedom.\r
- if (boundary < 0) {\r
- for (int sf = 0; sf < ndofs_per_cell; sf++) {\r
- int di = fe_v_neighbor.get_fe().system_to_component_index(sf).first;\r
- Wminus[q][di] +=\r
- (theta*DOF[sf+dofs_per_cell]+(1.0-theta)*solution(dofs_neighbor[sf]))*\r
- fe_v_neighbor.shape_value_component(sf, q, di);\r
- }\r
- } \r
- } // for q\r
-\r
- // If this is a boundary, then the values of $W^-$ will\r
- // be either functions of $W^+$, or they will be prescribed.\r
- // This switch sets them appropriately. Since we are\r
- // using fad variables here, sensitivities will be updated \r
- // appropriately. These sensitivities would be tremendously\r
- // difficult to manage without fad!!!\r
- if (boundary >= 0) {\r
- // Get the boundary descriptor.\r
- typename bdry_map_type::iterator bme = bdry_map.find(boundary);\r
- assert(bme != bdry_map.end());\r
-\r
- // Evaluate the function object. This is a bit\r
- // tricky; a given boundary might have both prescribed\r
- // and implicit values. If a particular component is not\r
- // prescribed, the values evaluate to zero and are\r
- // ignored, below.\r
- std::vector<Vector<double> > bvals(n_q_points, Vector<double>(N_COMP));\r
- bme->second.second->vector_value_list(fe_v.get_quadrature_points(), bvals);\r
-\r
- // We loop the quadrature points, and we treat each\r
- // component individualy.\r
- for (int q = 0; q < n_q_points; q++) {\r
- for (int di = 0; di < get_n_components(); di++) {\r
-\r
- // An inflow/dirichlet type of boundary condition\r
- if (bme->second.first[di] == INFLOW_BC) {\r
- Wminus[q][di] = bvals[q](di);\r
- } else if (bme->second.first[di] == PRESSURE_BC) {\r
- // A prescribed pressure boundary condition. This boundary\r
- // condition is complicated by the fact that even though\r
- // the pressure is prescribed, we really are setting\r
- // the energy index here, which will depend on velocity\r
- // and pressure. So even though this seems like a dirichlet\r
- // type boundary condition, we get sensitivities of\r
- // energy to velocity and density (unless these\r
- // are also prescribed.\r
- fad_double rho_vel_sqr = 0;\r
- fad_double dens;\r
- \r
- dens = bme->second.first[DENS_IDX] == INFLOW_BC ? bvals[q](DENS_IDX) :\r
- Wplus[q][DENS_IDX];\r
-\r
- for (int d=0; d < dim; d++) {\r
- if (bme->second.first[d] == INFLOW_BC)\r
- rho_vel_sqr += bvals[q](d)*bvals[q](d);\r
- else\r
- rho_vel_sqr += Wplus[q][d]*Wplus[q][d];\r
- }\r
- rho_vel_sqr /= dens;\r
- // Finally set the energy value as determined by the\r
- // prescribed pressure and the other variables.\r
- Wminus[q][di] = bvals[q](di)/(GAMMA-1.0) +\r
- 0.5*rho_vel_sqr;\r
-\r
- } else if (bme->second.first[di] == OUTFLOW_BC) {\r
- // A free/outflow boundary, very simple.\r
- Wminus[q][di] = Wplus[q][di];\r
-\r
- } else { \r
- // We must be at a no-penetration boundary. We\r
- // prescribe the velocity (we are dealing with a\r
- // particular component here so that the average\r
- // of the velocities is orthogonal to the surface\r
- // normal. This creates sensitivies of across\r
- // the velocity components.\r
- fad_double vdotn = 0;\r
- for (int d = 0; d < dim; d++) {\r
- vdotn += Wplus[q][d]*normals[q](d);\r
- }\r
-\r
- Wminus[q][di] = Wplus[q][di] - 2.0*vdotn*normals[q](di);\r
- }\r
- }\r
- } // for q\r
- } // b>= 0\r
- \r
- // Determine the Lax-Friedrich's stability parameter,\r
- // and evaluate the numerical flux function at the quadrature points\r
- std::vector<std::vector<fad_double> > nflux(n_q_points, std::vector<fad_double>(get_n_components(), 0));\r
- double alpha = 1;\r
-\r
- switch(flux_params.LF_stab) {\r
- case flux_params_type::CONSTANT:\r
- alpha = flux_params.LF_stab_value;\r
- break;\r
- case flux_params_type::MESH:\r
- alpha = face_diameter/(2.0*dT);\r
- break;\r
- }\r
-\r
- LFNumFlux<fad_double, dim>(nflux, fe_v.get_quadrature_points(), normals, Wplus, Wminus,\r
- alpha);\r
-\r
- // Now assemble the face term\r
- for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) {\r
- if (!fe_v.get_fe().has_support_on_face(i, face_no)) continue;\r
- F_i = 0;\r
- for (unsigned int point=0; point<n_q_points; ++point)\r
- {\r
- const unsigned int\r
- component_i = fe_v.get_fe().system_to_component_index(i).first;\r
-\r
- F_i += nflux[point][component_i]*fe_v.shape_value_component(i, point, component_i)*JxW[point];\r
-\r
- } \r
-\r
- // Retrieve a pointer to the jacobian.\r
- double *values = &(F_i.fastAccessDx(0));\r
-\r
- // Honestly, I forget why this can happen, but \r
- // for some reason it can!!\r
- if (!values) continue;\r
-\r
- // Update the matrix. Depending on whether there\r
- // is/isn't a neighboring cell, we add more/less\r
- // entries.\r
- Matrix->SumIntoGlobalValues(dofs[i],\r
- dofs_per_cell, &values[0], reinterpret_cast<int*>(&dofs[0]));\r
- if (boundary < 0) {\r
- Matrix->SumIntoGlobalValues(dofs[i],\r
- dofs_per_cell, &values[dofs_per_cell], reinterpret_cast<int*>(&dofs_neighbor[0]));\r
- }\r
-\r
- // And add into the residual\r
- right_hand_side(dofs[i]) -= F_i.val();\r
- } \r
-\r
-}\r
- // @sect4{Assembling the whole system}\r
- // Now we put all of the assembly pieces together\r
- // in a routine that dispatches the correct\r
- // piece for each cell/face. We keep track of\r
- // the norm of the resdual for the Newton iteration.\r
-template <int dim>\r
-void ConsLaw<dim>::assemble_system (double &res_norm) \r
-{\r
- FESystem<dim> &fe = *fe_ptr;\r
- const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;\r
-\r
- // We track the dofs on this cell and (if necessary)\r
- // the adjacent cell.\r
- std::vector<unsigned int> dofs (dofs_per_cell);\r
- std::vector<unsigned int> dofs_neighbor (dofs_per_cell);\r
-\r
- // First we create the\r
- // ``UpdateFlags'' for the\r
- // ``FEValues'' and the\r
- // ``FEFaceValues'' objects.\r
- UpdateFlags update_flags = update_values\r
- | update_gradients\r
- | update_q_points\r
- | update_JxW_values;\r
-\r
- // Note, that on faces we do not\r
- // need gradients but we need\r
- // normal vectors.\r
- UpdateFlags face_update_flags = update_values\r
- | update_q_points\r
- | update_JxW_values\r
- | update_normal_vectors;\r
- \r
- // On the neighboring cell we only\r
- // need the shape values. Given a\r
- // specific face, the quadrature\r
- // points and `JxW values' are the\r
- // same as for the current cells,\r
- // the normal vectors are known to\r
- // be the negative of the normal\r
- // vectors of the current cell.\r
- UpdateFlags neighbor_face_update_flags = update_values;\r
- \r
- // Then we create the ``FEValues''\r
- // object. Note, that since version\r
- // 3.2.0 of deal.II the constructor\r
- // of this class takes a\r
- // ``Mapping'' object as first\r
- // argument. Although the\r
- // constructor without ``Mapping''\r
- // argument is still supported it\r
- // is recommended to use the new\r
- // constructor. This reduces the\r
- // effect of `hidden magic' (the\r
- // old constructor implicitely\r
- // assumes a ``MappingQ1'' mapping)\r
- // and makes it easier to change\r
- // the mapping object later.\r
- FEValues<dim> fe_v (\r
- mapping, fe, quadrature, update_flags);\r
- \r
- // Similarly we create the\r
- // ``FEFaceValues'' and\r
- // ``FESubfaceValues'' objects for\r
- // both, the current and the\r
- // neighboring cell. Within the\r
- // following nested loop over all\r
- // cells and all faces of the cell\r
- // they will be reinited to the\r
- // current cell and the face (and\r
- // subface) number.\r
- FEFaceValues<dim> fe_v_face (\r
- mapping, fe, face_quadrature, face_update_flags);\r
- FESubfaceValues<dim> fe_v_subface (\r
- mapping, fe, face_quadrature, face_update_flags);\r
- FEFaceValues<dim> fe_v_face_neighbor (\r
- mapping, fe, face_quadrature, neighbor_face_update_flags);\r
- FESubfaceValues<dim> fe_v_subface_neighbor (\r
- mapping, fe, face_quadrature, neighbor_face_update_flags);\r
-\r
- // Furthermore we need some cell\r
- // iterators.\r
- typename DoFHandler<dim>::active_cell_iterator\r
- cell = dof_handler.begin_active(),\r
- endc = dof_handler.end();\r
-\r
- // Now we start the loop over all\r
- // active cells.\r
- int fdofs_per_cell = fe_v.dofs_per_cell;\r
- int fn_q_points = face_quadrature.n_quadrature_points;\r
-\r
- unsigned int cell_no = 0;\r
- for (;cell!=endc; ++cell, ++cell_no) \r
- {\r
- \r
- // Now we reinit the ``FEValues''\r
- // object for the current cell\r
- fe_v.reinit (cell);\r
-\r
- // Collect the local dofs and\r
- // asssemble the cell term.\r
- cell->get_dof_indices (dofs);\r
-\r
- cell_diameter = cell->diameter();\r
-\r
- assemble_cell_term(fe_v,\r
- dofs,\r
- cell_no);\r
-\r
- // We use the DG style loop through faces\r
- // to determine if we need to apply a\r
- // 'hanging node' flux calculation or a boundary\r
- // computation.\r
- for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)\r
- {\r
- // First we set the face\r
- // iterator\r
- typename DoFHandler<dim>::face_iterator face=cell->face(face_no);\r
- face_diameter = face->diameter();\r
- \r
- if (face->at_boundary())\r
- {\r
- // We reinit the\r
- // ``FEFaceValues''\r
- // object to the\r
- // current face\r
- fe_v_face.reinit (cell, face_no);\r
-\r
- // and assemble the\r
- // corresponding face\r
- // terms. We send the same\r
- // fe_v and dofs as described\r
- // in the assembly routine.\r
- assemble_face_term(\r
- face_no, fe_v_face,\r
- fe_v_face,\r
- dofs,\r
- dofs,\r
- face->boundary_indicator());\r
- }\r
- else\r
- {\r
- // Now we are not on\r
- // the boundary of the\r
- // domain, therefore\r
- // there must exist a\r
- // neighboring cell.\r
- typename DoFHandler<dim>::cell_iterator neighbor=\r
- cell->neighbor(face_no);;\r
-\r
- if (face->has_children())\r
- {\r
- // case I: This cell refined compared to neighbor\r
-\r
- const unsigned int neighbor2=\r
- cell->neighbor_of_neighbor(face_no);\r
- \r
- \r
- // We loop over\r
- // subfaces\r
- for (unsigned int subface_no=0;\r
- subface_no<GeometryInfo<dim>::subfaces_per_face;\r
- ++subface_no)\r
- {\r
- typename DoFHandler<dim>::active_cell_iterator\r
- neighbor_child\r
- = cell->neighbor_child_on_subface (face_no, subface_no);\r
-\r
- face_diameter = neighbor_child->diameter(); // working on subface\r
- \r
- Assert (neighbor_child->face(neighbor2) == face->child(subface_no),\r
- ExcInternalError());\r
- Assert (!neighbor_child->has_children(), ExcInternalError());\r
-\r
- fe_v_subface.reinit (cell, face_no, subface_no);\r
- fe_v_face_neighbor.reinit (neighbor_child, neighbor2);\r
- neighbor_child->get_dof_indices (dofs_neighbor);\r
-\r
- // Assemble as if we are working with\r
- // a DG element.\r
- assemble_face_term(\r
- face_no, fe_v_subface,\r
- fe_v_face_neighbor,\r
- dofs,\r
- dofs_neighbor);\r
- \r
- }\r
- // End of ``if\r
- // (face->has_children())''\r
- }\r
- else\r
- {\r
- // We have no children, but \r
- // the neighbor cell may be refine\r
- // compared to use\r
- neighbor->get_dof_indices (dofs_neighbor);\r
- if (neighbor->level() != cell->level()) \r
- {\r
- // case II: This is refined compared to neighbor\r
- Assert(neighbor->level() < cell->level(), ExcInternalError());\r
- const std::pair<unsigned int, unsigned int> faceno_subfaceno=\r
- cell->neighbor_of_coarser_neighbor(face_no);\r
- const unsigned int neighbor_face_no=faceno_subfaceno.first,\r
- neighbor_subface_no=faceno_subfaceno.second;\r
-\r
- Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,\r
- neighbor_subface_no)\r
- == cell,\r
- ExcInternalError());\r
-\r
- // Reinit the\r
- // appropriate\r
- // ``FEFaceValues''\r
- // and assemble\r
- // the face\r
- // terms.\r
- fe_v_face.reinit (cell, face_no);\r
- fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no,\r
- neighbor_subface_no);\r
- \r
- assemble_face_term(\r
- face_no, fe_v_face,\r
- fe_v_subface_neighbor,\r
- dofs,\r
- dofs_neighbor);\r
-\r
- }\r
-\r
- } \r
- // End of ``face not at boundary'':\r
- }\r
- // End of loop over all faces:\r
- } \r
- \r
- // End iteration through cells.\r
- } \r
-\r
- // Notify Epetra that the matrix is done.\r
- Matrix->FillComplete();\r
-\r
- // Compute the nonlinear residual.\r
- res_norm = right_hand_side.l2_norm();\r
- \r
-}\r
-\r
- // Create a conservation law with some defaults.\r
-template <int dim>\r
-ConsLaw<dim>::ConsLaw ()\r
- :\r
- mapping (),\r
- fe_ptr(NULL),\r
- dof_handler (triangulation),\r
- quadrature (2),\r
- face_quadrature (2),\r
- T(0),\r
- dT(0.05),\r
- TF(10),\r
- Map(NULL),\r
- Matrix(NULL),\r
- is_stationary(false),\r
- theta(0.5)\r
-{}\r
-\r
- // At one time this example could work for both DG and\r
- // continuous finite elements. The choice was made here.\r
-template <int dim>\r
-void ConsLaw<dim>::build_fe() {\r
- fe_ptr = new FESystem<dim>(FE_Q<dim>(1), N_COMP);\r
-}\r
-\r
- // Bye bye Conservation law.\r
-template <int dim>\r
-ConsLaw<dim>::~ConsLaw () \r
-{\r
- dof_handler.clear ();\r
- delete fe_ptr;\r
-}\r
-\r
- // @sect3{Initialize System}\r
- // Sizes all of the vectors and sets up the\r
- // sparsity patter. This function is called at\r
- // the very beginning of a simulation. The function\r
- // <code> setup_system </code> repeats some of these\r
- // chores and is called after adaptivity in leiu\r
- // of this function.\r
-template <int dim>\r
-void ConsLaw<dim>::initialize_system ()\r
-{\r
- // First we need to distribute the\r
- // DoFs.\r
- dof_handler.clear();\r
- dof_handler.distribute_dofs (*fe_ptr);\r
- \r
- // Size all of the fields.\r
- solution.reinit (dof_handler.n_dofs());\r
- nlsolution.reinit (dof_handler.n_dofs());\r
- predictor.reinit (dof_handler.n_dofs());\r
- ppsolution.reinit (dof_handler.n_dofs());\r
- dsolution.reinit (dof_handler.n_dofs());\r
- right_hand_side.reinit (dof_handler.n_dofs());\r
- indicator.reinit(triangulation.n_active_cells());\r
-}\r
-\r
- // @sect3{Setup System}\r
- // We call this function to build the sparsity\r
- // and the matrix.\r
-template <int dim>\r
-void ConsLaw<dim>::setup_system ()\r
-{\r
-\r
- // The DoFs of a cell are coupled\r
- // with all DoFs of all neighboring\r
- // cells. Therefore the maximum\r
- // number of matrix entries per row\r
- // is needed when all neighbors of\r
- // a cell are once more refined\r
- // than the cell under\r
- // consideration.\r
- sparsity_pattern.reinit (dof_handler.n_dofs(),\r
- dof_handler.n_dofs(),\r
- (GeometryInfo<dim>::faces_per_cell\r
- *GeometryInfo<dim>::subfaces_per_face+1)*fe_ptr->dofs_per_cell);\r
- \r
- // Since the continuous sparsity pattern is\r
- // a subset of the DG one, and since we need\r
- // the DG terms for handling hanging nodes, we use\r
- // the flux pattern.\r
- DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern);\r
- \r
- sparsity_pattern.compress();\r
- \r
- // Rebuild the map. In serial this doesn't do much,\r
- // but is needed. In parallel, this would desribe\r
- // the parallel dof layout.\r
- if (Map) delete Map;\r
- Map = new Epetra_Map(dof_handler.n_dofs(), 0, *Comm);\r
-\r
- // Epetra can build a more efficient matrix if\r
- // one knows ahead of time the maxiumum number of\r
- // columns in any row entry. We traverse the sparsity\r
- // to discover this.\r
- int cur_row = 0;\r
- int cur_col = 0;\r
- int max_df = -1;\r
- for (SparsityPattern::iterator s_i = sparsity_pattern.begin(); \r
- s_i != sparsity_pattern.end(); s_i++) {\r
- if (s_i->row() != cur_row) {\r
- cur_col = 0;\r
- cur_row = s_i->row();\r
- }\r
- cur_col++;\r
- if (cur_col >= max_df) max_df = cur_col;\r
- }\r
-\r
- if (cur_col >= max_df) max_df = cur_col;\r
- std::cout << "max_df:" << max_df << std::endl;\r
-\r
- // Now we build the matrix, using the constructor\r
- // that optimizes with the <code> max_df </code> variable.\r
- if (Matrix) delete Matrix;\r
- Matrix = new Epetra_CrsMatrix(Copy, *Map, max_df+1, true);\r
-\r
- // We add the sparsity pattern to the matrix by\r
- // inserting zeros.\r
- std::vector<double> vals(max_df, 0);\r
- std::vector<int> row_indices(max_df);\r
- \r
- cur_row = 0;\r
- cur_col = 0;\r
- for (SparsityPattern::iterator s_i = sparsity_pattern.begin(); \r
- s_i != sparsity_pattern.end(); s_i++) {\r
- if (s_i->row() != cur_row) {\r
- Matrix->InsertGlobalValues(cur_row, cur_col, &vals[0], &row_indices[0]);\r
- cur_col = 0;\r
- cur_row = s_i->row();\r
- }\r
- row_indices[cur_col++] = s_i->column();\r
- }\r
- // The last row.\r
- Matrix->InsertGlobalValues(cur_row, cur_col, &vals[0], &row_indices[0]);\r
-\r
- // Epetra requires this function after building or\r
- // filling a matrix. It typically does some parallel\r
- // bookeeping; perhaps more.\r
- Matrix->FillComplete();\r
-\r
-}\r
-\r
- // @sect3{Solving the linear system}\r
- // Actually solve the linear system, using either\r
- // Aztec of Amesos.\r
-template <int dim>\r
-void ConsLaw<dim>::solve (Vector<double> &dsolution, int &niter, double &lin_residual) \r
-{\r
-\r
- // We must hand the solvers Epetra vectors.\r
- // Luckily, they support the concept of a \r
- // 'view', so we just send in a pointer to our\r
- // dealii vectors.\r
- Epetra_Vector x(View, *Map, dsolution.begin());\r
- Epetra_Vector b(View, *Map, right_hand_side.begin());\r
-\r
- // The Direct option selects the Amesos solver.\r
- if (solver_params.SOLVER == solver_params_type::DIRECT) {\r
- \r
- // Setup for solving with\r
- // Amesos.\r
- Epetra_LinearProblem prob;\r
- prob.SetOperator(Matrix);\r
- Amesos_BaseSolver *solver;\r
- Amesos Factory;\r
-\r
- // Other solvers are available\r
- // and may be selected by changing this\r
- // string.\r
- char *stype = "Amesos_Klu";\r
-\r
- solver = Factory.Create(stype, prob);\r
-\r
- Assert (solver != NULL, ExcInternalError());\r
-\r
- // There are two parts to the direct solve.\r
- // As I understand, the symbolic part figures\r
- // out the sparsity patterns, and then the\r
- // numerical part actually performs Gaussian\r
- // elimination or whatever the approach is.\r
- if (solver_params.OUTPUT == solver_params_type::VERBOSE)\r
- std::cout << "Starting Symbolic fact\n" << std::flush;\r
-\r
- solver->SymbolicFactorization();\r
-\r
- if (solver_params.OUTPUT == solver_params_type::VERBOSE)\r
- std::cout << "Starting Numeric fact\n" << std::flush;\r
-\r
- solver->NumericFactorization();\r
-\r
- \r
- // Define the linear problem by setting the\r
- // right hand and left hand sides.\r
- prob.SetRHS(&b);\r
- prob.SetLHS(&x);\r
- // And finally solve the problem.\r
- if (solver_params.OUTPUT == solver_params_type::VERBOSE)\r
- std::cout << "Starting solve\n" << std::flush;\r
- solver->Solve();\r
- niter = 0;\r
- lin_residual = 0;\r
-\r
- // We must free the solver that was created\r
- // for us.\r
- delete solver;\r
-\r
- } else if (solver_params.SOLVER == solver_params_type::GMRES) {\r
-\r
- // For the iterative solvers, we use Aztec.\r
- AztecOO Solver;\r
-\r
- // Select the appropriate level of verbosity.\r
- if (solver_params.OUTPUT == solver_params_type::QUIET)\r
- Solver.SetAztecOption(AZ_output, AZ_none);\r
-\r
- if (solver_params.OUTPUT == solver_params_type::VERBOSE)\r
- Solver.SetAztecOption(AZ_output, AZ_all);\r
-\r
- // Select gmres. Other solvers are available.\r
- Solver.SetAztecOption(AZ_solver, AZ_gmres);\r
- Solver.SetRHS(&b);\r
- Solver.SetLHS(&x);\r
-\r
- // Set up the ILUT preconditioner. I do not know\r
- // why, but we must pretend like we are in parallel\r
- // using domain decomposition or the preconditioner\r
- // refuses to activate.\r
- Solver.SetAztecOption(AZ_precond, AZ_dom_decomp);\r
- Solver.SetAztecOption(AZ_subdomain_solve, AZ_ilut);\r
- Solver.SetAztecOption(AZ_overlap, 0);\r
- Solver.SetAztecOption(AZ_reorder, 0);\r
-\r
- // ILUT parameters as described above.\r
- Solver.SetAztecParam(AZ_drop, solver_params.ILUT_DROP);\r
- Solver.SetAztecParam(AZ_ilut_fill, solver_params.ILUT_FILL);\r
- Solver.SetAztecParam(AZ_athresh, solver_params.ILUT_ATOL);\r
- Solver.SetAztecParam(AZ_rthresh, solver_params.ILUT_RTOL);\r
- Solver.SetUserMatrix(Matrix);\r
-\r
- // Run the solver iteration. Collect the number\r
- // of iterations and the residual.\r
- Solver.Iterate(solver_params.MAX_ITERS, solver_params.RES);\r
- niter = Solver.NumIters();\r
- lin_residual = Solver.TrueResidual();\r
- }\r
-}\r
-\r
- // @sect3{Postprocessing and Output}\r
- // Recover the physical variables from the conservative\r
- // variables so that output will be (perhaps) more\r
- // meaningfull.\r
-template <int dim>\r
-void ConsLaw<dim>::postprocess() {\r
- const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;\r
- std::vector<unsigned int> dofs (dofs_per_cell);\r
- UpdateFlags update_flags = update_values\r
- | update_gradients\r
- | update_q_points\r
- | update_JxW_values;\r
- UpdateFlags update_flags1 = update_values\r
- | update_gradients\r
- | update_q_points\r
- | update_JxW_values;\r
-\r
- QGauss<dim> quadrature_formula(4);\r
-\r
- const std::vector<Point<dim> > &us = fe_ptr->base_element(0).get_unit_support_points();\r
-\r
-\r
- Quadrature<dim> unit_support(us);\r
-\r
- int n_q_points = quadrature_formula.n_quadrature_points;\r
- int n_uq_points = unit_support.n_quadrature_points;\r
-\r
- FEValues<dim> fe_v (\r
- mapping, *fe_ptr, quadrature_formula, update_flags);\r
-\r
- FEValues<dim> fe_v_unit (\r
- mapping, *fe_ptr, unit_support, update_flags1);\r
-\r
- std::vector<Vector<double> > U(n_uq_points,\r
- Vector<double>(get_n_components()));\r
- std::vector<Vector<double> > UU(n_q_points,\r
- Vector<double>(get_n_components()));\r
- std::vector<std::vector<Tensor<1,dim> > > dU(n_uq_points,\r
- std::vector<Tensor<1,dim> >(get_n_components()));\r
- \r
- typename DoFHandler<dim>::active_cell_iterator\r
- cell = dof_handler.begin_active(),\r
- endc = dof_handler.end();\r
-\r
- // Loop the cells\r
- for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) {\r
- cell->get_dof_indices (dofs);\r
- fe_v_unit.reinit(cell);\r
- fe_v.reinit(cell);\r
-\r
- fe_v_unit.get_function_values(solution, U);\r
- fe_v_unit.get_function_grads(solution, dU);\r
- fe_v.get_function_values(solution, UU);\r
-\r
- const std::vector<double> &JxW = fe_v.get_JxW_values ();\r
-\r
- for (int q = 0; q < fe_v.get_fe().base_element(0).n_dofs_per_cell(); q++) {\r
- unsigned int didx = fe_v.get_fe().component_to_system_index(DENS_IDX, q);\r
- unsigned int eidx = fe_v.get_fe().component_to_system_index(ENERGY_IDX, q);\r
- double rho_normVsqr = 0;\r
- for (int d = 0; d < dim; d++) {\r
- unsigned int vidx = fe_v.get_fe().component_to_system_index(d, q);\r
- ppsolution(dofs[vidx]) = solution(dofs[vidx])/solution(dofs[didx]);\r
- rho_normVsqr += solution(dofs[vidx])*solution(dofs[vidx]);\r
- }\r
- rho_normVsqr /= solution(dofs[didx]);\r
- // Pressure\r
- ppsolution(dofs[eidx]) = (GAMMA-1.0)*(solution(dofs[eidx]) - 0.5*rho_normVsqr);\r
-\r
- // Either output density or gradient squared of density,\r
- // depending on what the user wants.\r
- if (!schlieren_plot) {\r
- ppsolution(dofs[didx]) = solution(dofs[didx]);\r
- } else {\r
- double ng = 0;\r
- for (int i = 0; i < dim; i++) ng += dU[q][DENS_IDX][i]*dU[q][DENS_IDX][i];\r
- ng = std::sqrt(ng);\r
- ppsolution(dofs[didx]) = ng;\r
- }\r
- }\r
-\r
- } // cell\r
-\r
-}\r
-\r
- // Loop and assign a value for refinement. We\r
- // simply use the density squared, which selects\r
- // shocks with some success.\r
-template <int dim>\r
-void ConsLaw<dim>::estimate() {\r
- \r
- const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;\r
- std::vector<unsigned int> dofs (dofs_per_cell);\r
- UpdateFlags update_flags = update_values\r
- | update_gradients\r
- | update_q_points\r
- | update_JxW_values;\r
-\r
- QGauss<dim> quadrature_formula(1);\r
- int n_q_points = quadrature_formula.n_quadrature_points;\r
-\r
-\r
- FEValues<dim> fe_v (\r
- mapping, *fe_ptr, quadrature_formula, update_flags);\r
-\r
- std::vector<Vector<double> > U(n_q_points,\r
- Vector<double>(get_n_components()));\r
- std::vector<std::vector<Tensor<1,dim> > > dU(n_q_points,\r
- std::vector<Tensor<1,dim> >(get_n_components()));\r
- \r
- typename DoFHandler<dim>::active_cell_iterator\r
- cell = dof_handler.begin_active(),\r
- endc = dof_handler.end();\r
- for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) {\r
- fe_v.reinit(cell);\r
-\r
- fe_v.get_function_values(predictor, U);\r
- fe_v.get_function_grads(predictor, dU);\r
-\r
- indicator(cell_no) = 0;\r
- for (int q = 0; q < n_q_points; q++) {\r
- double ng = 0;\r
- for (int d = 0; d < dim; d++) ng += dU[q][DENS_IDX][d]*dU[q][DENS_IDX][d];\r
-\r
- indicator(cell_no) += std::log(1+std::sqrt(ng));\r
- \r
- } \r
- indicator(cell_no) /= n_q_points;\r
-\r
- } \r
-}\r
-\r
-template <int dim>\r
-void ConsLaw<dim>::refine_grid ()\r
-{\r
-\r
- SolutionTransfer<dim, double> soltrans(dof_handler);\r
-\r
- typename DoFHandler<dim>::active_cell_iterator\r
- cell = dof_handler.begin_active(),\r
- endc = dof_handler.end();\r
-\r
- // Loop cells. If the indicator\r
- // for the cell matches the refinement criterion,\r
- // refine, else unrefine. The unrefinement has\r
- // a slight hysterisis to avoid 'flashing' from refined\r
- // to unrefined.\r
- for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) {\r
- cell->clear_coarsen_flag();\r
- cell->clear_refine_flag();\r
- if (cell->level() < refinement_params.shock_levels &&\r
- std::fabs(indicator(cell_no)) > refinement_params.shock_val ) {\r
- cell->set_refine_flag();\r
- } else {\r
- if (cell->level() > 0 &&\r
- std::fabs(indicator(cell_no)) < 0.75*refinement_params.shock_val)\r
- cell->set_coarsen_flag();\r
- }\r
- }\r
-\r
- // The following code prolongs the solution\r
- // to the new grid and carries out the refinement.\r
- std::vector<Vector<double> > interp_in;\r
- std::vector<Vector<double> > interp_out;\r
-\r
- interp_in.push_back(solution);\r
- interp_in.push_back(predictor);\r
-\r
- triangulation.prepare_coarsening_and_refinement();\r
- soltrans.prepare_for_coarsening_and_refinement(interp_in);\r
-\r
- triangulation.execute_coarsening_and_refinement ();\r
-\r
- dof_handler.clear();\r
- dof_handler.distribute_dofs (*fe_ptr);\r
-\r
- {\r
- Vector<double> new_solution(1);\r
- Vector<double> new_predictor(1);\r
-\r
- interp_out.push_back(new_solution);\r
- interp_out.push_back(new_predictor);\r
- interp_out[0].reinit(dof_handler.n_dofs());\r
- interp_out[1].reinit(dof_handler.n_dofs());\r
- }\r
-\r
- soltrans.interpolate(interp_in, interp_out);\r
- \r
- // Let the vector delete a very small vector\r
- solution.reinit(1);\r
- predictor.reinit(1);\r
- solution.swap(interp_out[0]);\r
- predictor.swap(interp_out[1]);\r
-\r
- // resize these vectors for the new grid.\r
- nlsolution.reinit(dof_handler.n_dofs());\r
- ppsolution.reinit(dof_handler.n_dofs());\r
- nlsolution = solution;\r
- dsolution.reinit (dof_handler.n_dofs());\r
- right_hand_side.reinit (dof_handler.n_dofs());\r
-\r
- indicator.reinit(triangulation.n_active_cells());\r
-\r
-}\r
-\r
-template <int dim>\r
-void ConsLaw<dim>::output_results (const unsigned int cycle) const\r
-{\r
- char filename[512];\r
- std::sprintf(filename, "solution-%03d.vtk", cycle);\r
- std::ofstream output (filename);\r
-\r
- DataOut<dim> data_out;\r
- data_out.attach_dof_handler (dof_handler);\r
- std::vector<std::string> solution_names;\r
-\r
- // Rename the output with the physical variable\r
- // names. Send the post-processed values.\r
- solution_names.clear();\r
- for (int i = 0; i < dim; i++) {\r
- char buf[512];\r
- std::sprintf(buf, "v_%d", i);\r
- solution_names.push_back (buf); \r
- }\r
- solution_names.push_back("density");\r
- solution_names.push_back("pressure");\r
- data_out.add_data_vector (ppsolution, solution_names);\r
-\r
- data_out.add_data_vector (indicator, "error");\r
- data_out.build_patches ();\r
- data_out.write_vtk (output);\r
-\r
- output.close();\r
-}\r
-\r
- // @sect3{Parsing the Input Deck}\r
- // Declare the parameters for the\r
- // input deck. We assume a certain\r
- // maximum number of boundaries and process\r
- // any boundary the user supplies up to\r
- // that maximum number. We\r
- // leave a detailed explanation of these \r
- // parameters to our description of the input\r
- // sample file.\r
-const UInt MAX_BD = 10;\r
-template <int dim>\r
-void ConsLaw<dim>::declare_parameters() {\r
-\r
- // Global scope parameters/\r
- prm.declare_entry("mesh", "grid.inp",\r
- Patterns::Anything(),\r
- "intput file");\r
-\r
- prm.declare_entry("diffusion power", "2.0",\r
- Patterns::Double(),\r
- "power of mesh size for diffusion");\r
-\r
- prm.declare_entry("gravity", "0.0",\r
- Patterns::Double(),\r
- "gravity forcing");\r
-\r
- // Time stepping block\r
- prm.enter_subsection("time stepping");\r
- prm.declare_entry("time step", "0.1",\r
- Patterns::Double(),\r
- "simulation time step");\r
- prm.declare_entry("final time", "10.0",\r
- Patterns::Double(),\r
- "simulation end time");\r
- prm.leave_subsection();\r
-\r
-\r
- // Declare the boundary parameters\r
- for (int b = 0; b < MAX_BD; b++) {\r
- char bd[512];\r
- std::sprintf(bd, "boundary_%d", b);\r
- prm.enter_subsection(bd);\r
- prm.declare_entry("no penetration", "false",\r
- Patterns::Selection("true|false"),\r
- "<true|false>");\r
- // declare a slot for each of the conservative\r
- // variables.\r
- for (int di = 0; di < N_COMP; di++) {\r
- char var[512];\r
- std::sprintf(var, "w_%d", di);\r
- prm.declare_entry(var, "outflow",\r
- Patterns::Selection(\r
- "inflow|outflow|pressure"),\r
- "<inflow|outflow|pressure>");\r
- \r
- // for dirichlet, a function in x,y,z\r
- std::sprintf(var, "w_%d value", di);\r
- prm.declare_entry(var, "0.0",\r
- Patterns::Anything(),\r
- "expression in x,y,z");\r
- }\r
-\r
- prm.leave_subsection();\r
- }\r
-\r
- // Initial condition block.\r
- prm.enter_subsection("initial condition");\r
- for (int di = 0; di < N_COMP; di++) {\r
- char var[512];\r
- std::sprintf(var, "w_%d", di);\r
- \r
- // for dirichlet, a function in x,y,z\r
- std::sprintf(var, "w_%d value", di);\r
- prm.declare_entry(var, "0.0",\r
- Patterns::Anything(),\r
- "expression in x,y,z");\r
- }\r
- prm.leave_subsection();\r
-\r
- // The linear solver block.\r
- prm.enter_subsection("linear solver");\r
- prm.declare_entry("output", "quiet",\r
- Patterns::Selection(\r
- "quiet|verbose"),\r
- "<quiet|verbose>");\r
- prm.declare_entry("method", "gmres",\r
- Patterns::Selection(\r
- "gmres|direct"),\r
- "<gmres|direct>");\r
- prm.declare_entry("residual", "1e-10",\r
- Patterns::Double(),\r
- "linear solver residual");\r
- prm.declare_entry("max iters", "300",\r
- Patterns::Double(),\r
- "maximum solver iterations");\r
- prm.declare_entry("ilut fill", "2",\r
- Patterns::Double(),\r
- "ilut preconditioner fill");\r
- prm.declare_entry("ilut absolute tolerance", "1e-9",\r
- Patterns::Double(),\r
- "ilut preconditioner tolerance");\r
- prm.declare_entry("ilut relative tolerance", "1.1",\r
- Patterns::Double(),\r
- "rel tol");\r
- prm.declare_entry("ilut drop tolerance", "1e-10",\r
- Patterns::Double(),\r
- "ilut drop tol");\r
- prm.leave_subsection();\r
-\r
-\r
- // A refinement controller block.\r
- prm.enter_subsection("refinement");\r
- prm.declare_entry("refinement", "none",\r
- Patterns::Selection(\r
- "none|fixed number|shock"),\r
- "<on|off>");\r
- prm.declare_entry("refinement fraction", "0.1",\r
- Patterns::Double(),\r
- "Fraction of high refinement");\r
- prm.declare_entry("unrefinement fraction", "0.1",\r
- Patterns::Double(),\r
- "Fraction of low unrefinement");\r
- prm.declare_entry("max elements", "1000000",\r
- Patterns::Double(),\r
- "maximum number of elements");\r
- prm.declare_entry("shock value", "4.0",\r
- Patterns::Double(),\r
- "value for shock indicator");\r
- prm.declare_entry("shock levels", "3.0",\r
- Patterns::Double(),\r
- "number of shock refinement levels");\r
- prm.leave_subsection();\r
-\r
- // Output control.\r
- prm.enter_subsection("output");\r
- prm.declare_entry("density", "standard",\r
- Patterns::Selection(\r
- "standard|schlieren"),\r
- "<standard|schlieren>");\r
- prm.declare_entry("step", "-1",\r
- Patterns::Double(),\r
- "output once per this period");\r
- prm.leave_subsection();\r
-\r
- // Flux control\r
- prm.enter_subsection("flux");\r
- prm.declare_entry("stab", "alpha",\r
- Patterns::Selection(\r
- "alpha|constant|mesh"),\r
- "<alpha|constant|mesh>");\r
- prm.declare_entry("stab value", "1",\r
- Patterns::Double(),\r
- "alpha stabilization");\r
- prm.leave_subsection();\r
-\r
-\r
-}\r
-\r
- // Code to actually parse an input file. This function\r
- // matches the declarations above.\r
-template <int dim>\r
-void ConsLaw<dim>::load_parameters(const char *infile){\r
-\r
- prm.read_input(infile);\r
-\r
- // The global parameters.\r
- mesh = prm.get("mesh");\r
-\r
- diffusion_power = prm.get_double("diffusion power");\r
-\r
- gravity = prm.get_double("gravity");\r
-\r
- // The time stepping.\r
- prm.enter_subsection("time stepping");\r
- dT = prm.get_double("time step");\r
- std::cout << "dT=" << dT << std::endl;\r
- if (dT == 0) {\r
- is_stationary = true;\r
- dT = 1.0;\r
- TF = 1.0;\r
- std::cout << "Stationary mode" << std::endl;\r
- }\r
- TF = prm.get_double("final time");\r
- std::cout << "TF=" << TF << std::endl;\r
- prm.leave_subsection();\r
-\r
- // The boundary info\r
- for (int b = 0; b < MAX_BD; b++) {\r
- std::vector<bc_type> flags(N_COMP, OUTFLOW_BC);\r
-\r
- // Define a parser for every boundary, though it may be\r
- // unused.\r
- SideCondition<dim> *sd = new SideCondition<dim>(N_COMP);\r
- char bd[512];\r
- std::sprintf(bd, "boundary_%d", b);\r
- prm.enter_subsection(bd);\r
-\r
- const std::string &nopen = prm.get("no penetration");\r
-\r
- // Determine how each component is handled.\r
- for (int di = 0; di < N_COMP; di++) {\r
- char var[512];\r
- std::sprintf(var, "w_%d", di);\r
- std::string btype = prm.get(var);\r
- std::sprintf(var, "w_%d value", di);\r
- std::string var_value = prm.get(var);\r
-\r
- if (di < dim && nopen == "true") {\r
- flags[di] = NO_PENETRATION_BC;\r
- } else if (btype == "inflow") {\r
- flags[di] = INFLOW_BC;\r
- sd->set_coeff_row(di, var_value); \r
- } else if (btype == "pressure") {\r
- flags[di] = PRESSURE_BC;\r
- sd->set_coeff_row(di, var_value); \r
- }\r
- } \r
- prm.leave_subsection();\r
-\r
- // Add the boundary condition to the law.\r
- sd->Init();\r
- add_boundary(b, flags, sd);\r
- }\r
-\r
- // Initial conditions.\r
- prm.enter_subsection("initial condition");\r
- for (int di = 0; di < N_COMP; di++) {\r
- char var[512];\r
-\r
- std::sprintf(var, "w_%d value", di);\r
- std::string var_value = prm.get(var);\r
- ic.set_ic(di, var_value); \r
- }\r
- ic.Init();\r
- prm.leave_subsection();\r
-\r
- // The linear solver.\r
- prm.enter_subsection("linear solver");\r
- const std::string &op = prm.get("output");\r
- if (op == "verbose") solver_params.OUTPUT = solver_params_type::VERBOSE;\r
- if (op == "quiet") solver_params.OUTPUT = solver_params_type::QUIET;\r
- const std::string &sv = prm.get("method");\r
- if (sv == "direct") {\r
- solver_params.SOLVER = solver_params_type::DIRECT;\r
- } else if (sv == "gmres") {\r
- solver_params.SOLVER = solver_params_type::GMRES;\r
- } \r
-\r
- solver_params.RES = prm.get_double("residual");\r
- solver_params.MAX_ITERS = (int) prm.get_double("max iters");\r
- solver_params.ILUT_FILL = prm.get_double("ilut fill");\r
- solver_params.ILUT_ATOL = prm.get_double("ilut absolute tolerance");\r
- solver_params.ILUT_RTOL = prm.get_double("ilut relative tolerance");\r
- solver_params.ILUT_DROP = prm.get_double("ilut drop tolerance");\r
- solver_params.RES = prm.get_double("residual");\r
- prm.leave_subsection();\r
-\r
-\r
- // And refiement.\r
- prm.enter_subsection("refinement");\r
- const std::string &ref = prm.get("refinement");\r
- if (ref == "none") {\r
- refinement_params.refine = refinement_params_type::NONE;\r
- } else if (ref == "fixed number") {\r
- refinement_params.refine = refinement_params_type::FIXED_NUMBER;\r
- } else if (ref == "shock") {\r
- refinement_params.refine = refinement_params_type::SHOCK;\r
- } else\r
- refinement_params.high_frac = prm.get_double("refinement fraction");\r
- refinement_params.high_frac_sav = refinement_params.high_frac;\r
- refinement_params.low_frac = prm.get_double("unrefinement fraction");\r
- refinement_params.max_cells = prm.get_double("max elements");\r
- refinement_params.shock_val = prm.get_double("shock value");\r
- refinement_params.shock_levels = prm.get_double("shock levels");\r
- prm.leave_subsection();\r
- \r
- // Output control.\r
- prm.enter_subsection("output");\r
- const std::string &dens = prm.get("density");\r
- schlieren_plot = dens == "schlieren" ? true : false;\r
- output_step = prm.get_double("step");\r
- prm.leave_subsection();\r
-\r
- // Flux control.\r
- prm.enter_subsection("flux");\r
- const std::string &stab = prm.get("stab");\r
- if (stab == "constant") {\r
- flux_params.LF_stab = flux_params_type::CONSTANT;\r
- } else if (stab == "mesh ") {\r
- flux_params.LF_stab = flux_params_type::MESH;\r
- }\r
- flux_params.LF_stab_value = prm.get_double("stab value");\r
- prm.leave_subsection();\r
-\r
-\r
-}\r
-\r
-template<int dim>\r
-void ConsLaw<dim>::zero_matrix() {\r
- Matrix->PutScalar(0); Matrix->FillComplete();\r
-}\r
-\r
- // We use a predictor to try and make adaptivity\r
- // work better. The idea is to try and refine ahead\r
- // of a front, rather than stepping into a coarse\r
- // set of elements and smearing the solution. This\r
- // simple time extrapolator does the job.\r
-template<int dim>\r
-void ConsLaw<dim>::compute_predictor() {\r
- predictor = nlsolution;\r
- predictor.sadd(3/2.0, -1/2.0, solution);\r
-}\r
-\r
- // @sect3{Run the simulation} Contains the initialization,\r
- // the time loop, and the inner Newton iteration.\r
-template <int dim>\r
-void ConsLaw<dim>::run () \r
-{\r
-\r
- // Open and load the mesh.\r
- GridIn<dim> grid_in;\r
- grid_in.attach_triangulation(triangulation);\r
- std::cout << "Opening mesh <" << mesh << ">" << std::endl;\r
- std::ifstream input_file(mesh.c_str(), std::ios::in);\r
-\r
- Assert (infile,\r
- ExcFileNotOpen());\r
-\r
- grid_in.read_ucd(input_file); \r
- input_file.close();\r
- \r
- build_fe();\r
-\r
- unsigned int nstep = 0;\r
- \r
- // Initialize fields and matrices.\r
- initialize_system (); \r
- setup_system();\r
- initialize(); \r
- predictor = solution;\r
-\r
- // Initial refinement. We apply the ic,\r
- // estimate, refine, and repeat until\r
- // happy.\r
- if (refinement_params.refine != refinement_params_type::NONE)\r
- for (int i = 0; i < refinement_params.shock_levels; i++) {\r
- estimate();\r
- refine_grid();\r
- setup_system();\r
- initialize(); \r
- predictor = solution;\r
- }\r
- postprocess();\r
- output_results (nstep);\r
-\r
- // Determine when we will output next.\r
- double next_output = T + output_step;\r
-\r
- // @sect4{Main time stepping loop}\r
- predictor = solution;\r
- while(T < TF)\r
- {\r
- std::cout << "T=" << T << ", ";\r
-\r
-\r
- std::cout << " Number of active cells: "\r
- << triangulation.n_active_cells()\r
- << std::endl;\r
-\r
-\r
- std::cout << " Number of degrees of freedom: "\r
- << dof_handler.n_dofs()\r
- << std::endl;\r
-\r
- bool nonlin_done = false;\r
- double res_norm;\r
- int lin_iter;\r
-\r
- // Print some relevant information during the\r
- // Newton iteration.\r
- std::cout << "NonLin Res: Lin Iter Lin Res" << std::endl;\r
- std::cout << "______________________________________" << std::endl;\r
-\r
- int max_nonlin = 7;\r
- int nonlin_iter = 0;\r
- double lin_res;\r
-\r
- // @sect5{Newton iteration}\r
- nlsolution = predictor;\r
- while (!nonlin_done) {\r
- lin_iter = 0;\r
- zero_matrix();\r
- right_hand_side = 0;\r
- assemble_system (res_norm);\r
- // Flash a star to the screen so one can\r
- // know when the assembly has stopped and the linear\r
- // solution is starting.\r
- std::cout << "* " << std::flush;\r
-\r
- // Test against a (hardcoded) nonlinear tolderance.\r
- // Do not solve the linear system at the last step \r
- // (since it would be a waste).\r
- \r
- if (fabs(res_norm) < 1e-10) {\r
- nonlin_done = true;\r
- } else {\r
- // Solve the linear system and update with the\r
- // delta.\r
- dsolution = 0;\r
- solve (dsolution, lin_iter, lin_res);\r
- nlsolution.add(1.0, dsolution);\r
- }\r
-\r
- // Print the residuals.\r
- std::printf("%-16.3e %04d %-5.2e\n",\r
- res_norm, lin_iter, lin_res);\r
-\r
- nonlin_iter++;\r
- } \r
-\r
- // Various post convergence tasks.\r
- compute_predictor();\r
-\r
- solution = nlsolution;\r
-\r
-\r
- estimate();\r
-\r
- postprocess();\r
-\r
- T += dT;\r
-\r
- // Output if it is time.\r
- if (output_step < 0) {\r
- output_results (++nstep);\r
- } else if (T >= next_output) {\r
- output_results (++nstep);\r
- next_output += output_step;\r
- }\r
-\r
- // Refine, if refinement is selected.\r
- if (refinement_params.refine != refinement_params_type::NONE) {\r
- refine_grid();\r
- setup_system();\r
- }\r
- }\r
-}\r
-\r
- // The following ``main'' function is\r
- // similar to previous examples and\r
- // need not to be commented on.\r
-int main (int argc, char *argv[]) \r
-{\r
-\r
- MPI_Init(&argc, &argv);\r
- Comm = new Epetra_MpiComm(MPI_COMM_WORLD);\r
-\r
- if (argc != 2) {\r
- std::cout << "Usage:" << argv[0] << " infile" << std::endl;\r
- std::exit(1);\r
- }\r
- try\r
- {\r
- ConsLaw<DIMENSION> cons;\r
- cons.declare_parameters();\r
- cons.load_parameters(argv[1]);\r
- cons.run ();\r
- }\r
- catch (std::exception &exc)\r
- {\r
- std::cerr << std::endl << std::endl\r
- << "----------------------------------------------------"\r
- << std::endl;\r
- std::cerr << "Exception on processing: " << std::endl\r
- << exc.what() << std::endl\r
- << "Aborting!" << std::endl\r
- << "----------------------------------------------------"\r
- << std::endl;\r
- return 1;\r
- }\r
- catch (...) \r
- {\r
- std::cerr << std::endl << std::endl\r
- << "----------------------------------------------------"\r
- << std::endl;\r
- std::cerr << "Unknown exception!" << std::endl\r
- << "Aborting!" << std::endl\r
- << "----------------------------------------------------"\r
- << std::endl;\r
- return 1;\r
- };\r
- \r
- return 0;\r
-}\r
-\r
+/* step-12.cc,v 1.33 2005/08/08 16:41:40 wolf Exp */
+/* Author: David Neckels, Boulder Colorado 2007 */
+/* step-12.cc,v 1.33 2005/08/08 16:41:40 wolf Exp */
+/* Version: Version-5-2-0 */
+/* */
+/* Copyright (C) 2001, 2002, 2003, 2004, 2005, 2008 by the deal.II authors */
+/* */
+/* This file is subject to QPL and may not be distributed */
+/* without copyright and license information. Please refer */
+/* to the file deal.II/doc/license.html for the text and */
+/* further information on this license. */
+
+ // This program solves the Euler equations
+ // of gas dynamics for a given configuration
+ // file. It uses a standard Galerkin approach
+ // with weakly applied boundary conditions.
+
+ // <h4>Include files</h4>
+
+ // Aztecoo require mpi (even though we run on only
+ // one processor in this example).
+#include <mpi.h>
+
+ // Here we have the necessary TRILINOS includes.
+ //
+ // Epetra is the basic trilinos vector/matrix library.
+#include <Epetra_MpiComm.h>
+#include <Epetra_Map.h>
+#include <Epetra_CrsGraph.h>
+#include <Epetra_CrsMatrix.h>
+#include <Epetra_Vector.h>
+ // Teuchos is a Trilinos utility library that is used
+ // to set parameters within the Aztec solver library.
+#include "Teuchos_ParameterList.hpp"
+ // Aztec is the iterative solver library.
+#include <AztecOO.h>
+#include <AztecOO_Operator.h>
+#define HAVE_IFPACK_TEUCHOS
+#include <Ifpack.h>
+
+ // Amesos is a direct solver package within Trilinos.
+#include <Amesos.h>
+ // Sacado is the automatic differentiation package, which
+ // is used to find the jacobian for a fully implicit Newton
+ // iteration.
+#include <Sacado.hpp>
+
+ // A standard set of dealii includes. Nothing special to
+ // comment on here.
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/parameter_handler.h>
+#include <base/function_parser.h>
+
+#include <lac/vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/vector_memory.h>
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_out.h>
+#include <grid/grid_refinement.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_in.h>
+
+#include <fe/fe_values.h>
+#include <fe/fe_system.h>
+
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+
+#include <numerics/data_out.h>
+#include <numerics/vectors.h>
+#include <numerics/solution_transfer.h>
+
+#include <fe/mapping_q1.h>
+#include <fe/fe_q.h>
+#include <numerics/derivative_approximation.h>
+ // And this again is C++:
+#include <iostream>
+#include <fstream>
+#include <vector>
+
+ // Introduce the dealii library into the current namespace.
+using namespace dealii;
+
+ // We define a shorter name for the automatic differentiation
+ // type.
+typedef Sacado::Fad::DFad<double> fad_double;
+typedef unsigned int UInt;
+ // The Epetra library requires a 'communicator', which describes
+ // the layout of a parallel (or serial) set of processors.
+Epetra_MpiComm *Comm;
+
+ // <h4>Flux function definition</h4>
+ // Here we define the flux function for this system of conservation
+ // laws. Note: it would be terribly difficult to use this example
+ // to solve some other system of conservation laws.
+ //
+ // We define the number of components in the system. Euler's has
+ // one entry for momenta in each spatial direction, plus the energy
+ // and density components.
+#define N_COMP (2 + DIMENSION)
+ // Define a handle to the density and energy indices. We have arrange
+ // the momenta to be first, then density, and, lastly, energy.
+#define DENS_IDX DIMENSION
+#define ENERGY_IDX (DIMENSION+1)
+
+ // The gas constant. This value is representative of air.
+const double GAMMA = 1.4;
+ // We define the flux functions as one large matrix. Each row of this
+ // matrix represents a scalar conservation law for the component in
+ // that row. We template the numerical type of the flux function
+ // so that we may use the automatic differentiation type here.
+ // The flux functions are defined in terms of the
+ // conserved variables $\rho w_0, \dots, \rho w_{d-1}, \rho, E$,
+ // so they do not look exactly like the Euler equations one is
+ // used to seeing. We evaluate the flux at a single quadrature
+ // point.
+template <typename number, int dim>
+void Flux(std::vector<std::vector<number> > &flux,
+ const Point<dim> &point,
+ const std::vector<number> &W)
+{
+
+ // Pressure is a dependent variable: $p =
+ // (\gamma - 1)(E-\frac{1}{2} \rho |v|^2)$.
+ number rho_normVsqr;
+ for (int d0 = 0; d0 < dim; d0++) rho_normVsqr += W[d0]*W[d0];
+ // Since W are $\rho v$, we get a $\rho^2$ in the
+ // numerator, so dividing a $\rho$ out gives the desired $ \rho |v|^2$.
+ rho_normVsqr /= W[DENS_IDX];
+
+ number pressure = (GAMMA-1.0)*(W[ENERGY_IDX] - number(0.5)*(rho_normVsqr));
+
+ // We compute the momentum terms. We divide by the
+ // density here to get $v_i \rho v_j$
+ for (int d = 0; d < dim; d++) {
+ for (int d1 = 0; d1 < dim; d1++) {
+ flux[d][d1] = W[d]*W[d1]/W[DENS_IDX];
+ }
+ // The pressure contribution, along the diagonal:
+ flux[d][d] += pressure;
+ // Advection/conservation of density:
+ flux[DENS_IDX][d] = W[d];
+ // And, lastly, conservation of energy.
+ flux[ENERGY_IDX][d] = W[d]/W[DENS_IDX]*
+ (W[ENERGY_IDX] + pressure); // energy
+ }
+}
+
+ // On the boundaries of the domain and across `hanging nodes` we use
+ // a numerical flux function to enforce boundary conditions. This routine
+ // is the basic Lax-Friedrich's flux with a stabilization parameter
+ // $\alpha$.
+template <typename number, int dim>
+void LFNumFlux(
+ std::vector<std::vector<fad_double> > &nflux,
+ const std::vector<Point<dim> > &points,
+ const std::vector<Point<dim> > &normals,
+ const std::vector<std::vector<number> > &Wplus,
+ const std::vector<std::vector<number> > &Wminus,
+ double alpha)
+{
+ int n_q_points = points.size();
+
+ // We evaluate the flux at each of the quadrature points.
+ for (int q = 0; q < n_q_points; q++) {
+ std::vector<std::vector<fad_double> > iflux(N_COMP,
+ std::vector<fad_double>(dim, 0));
+ std::vector<std::vector<fad_double> > oflux(N_COMP,
+ std::vector<fad_double>(dim, 0));
+
+ Flux<number, dim>(iflux, points[q], Wplus[q]);
+ Flux<number, dim>(oflux, points[q], Wminus[q]);
+
+ for (int di = 0; di < N_COMP; di++) {
+ nflux[q][di] = 0;
+ for (int d = 0; d < dim; d++) {
+ nflux[q][di] += 0.5*(iflux[di][d] + oflux[di][d])*normals[q](d);
+ }
+ nflux[q][di] += 0.5*alpha*(Wplus[q][di] - Wminus[q][di]);
+ }
+ }
+
+}
+
+ // <h4>Initial and side condition parsing</h4>
+ // For the initial condition we use the expression parser function
+ // object.
+template <int dim>
+class InitialCondition : public FunctionParser<dim>
+{
+ public:
+ InitialCondition ();
+
+ // This function should be called after parsing, but before using
+ // the object. It formalizes the expressions and initializes the
+ // function parser with the appropriate expressions.
+ void Init();
+
+ // During parsing we call this function as the initial condition
+ // for one of the $\mathbf{w}$ variables is encountered.
+
+ void set_ic(int _row, std::string &expr) {
+ expressions[_row] = expr;
+ }
+
+ virtual void vector_value_list (const std::vector<Point<dim> > &points,
+ std::vector<Vector<double> > &value_list) const;
+ private:
+ std::vector<std::string> expressions;
+};
+
+template <int dim>
+InitialCondition<dim>::InitialCondition () :
+ FunctionParser<dim> (N_COMP),
+ expressions(N_COMP, "0.0")
+{}
+
+ // Here we set up x,y,z as the variables that one should use in the input
+ // deck to describe their initial condition.
+template<int dim>
+void InitialCondition<dim>::Init() {
+ std::map<std::string, double> constants;
+ constants["M_PI"] = M_PI;
+ std::string variables = (dim == 2 ? "x,y" : "x,y,z");
+
+ FunctionParser<dim>::initialize(variables, expressions, constants);
+
+}
+
+template <int dim>
+void InitialCondition<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+ std::vector<Vector<double> > &value_list) const
+{
+ const unsigned int n_points = points.size();
+
+ Assert (value_list.size() == n_points,
+ ExcDimensionMismatch (value_list.size(), n_points));
+
+ for (unsigned int p=0; p<n_points; ++p)
+ InitialCondition<dim>::vector_value (points[p],
+ value_list[p]);
+}
+
+ // As above, we use the expression function parser for boundary conditions.
+template <int dim>
+class SideCondition : public FunctionParser<dim>
+{
+ public:
+ SideCondition (int ncomp);
+ ~SideCondition ();
+
+ // As above.
+ void Init();
+ // As above.
+ void set_coeff_row(int _row_n, std::string &expr);
+
+ virtual void vector_value_list (const std::vector<Point<dim> > &points,
+ std::vector<Vector<double> > &value_list) const;
+ private:
+ std::vector<std::string> expressions;
+};
+
+template <int dim>
+SideCondition<dim>::SideCondition (int ncomp) :
+ FunctionParser<dim> (ncomp),
+ expressions(ncomp, "0.0")
+{
+}
+template <int dim>
+void SideCondition<dim>::set_coeff_row (int _row_n, std::string &expr)
+{
+ expressions[_row_n] = expr;
+}
+
+template <int dim>
+void SideCondition<dim>::Init() {
+ std::map<std::string, double> constants;
+ constants["M_PI"] = M_PI;
+ std::string variables = (dim == 2 ? "x,y" : "x,y,z");
+
+ FunctionParser<dim>::initialize(variables, expressions, constants);
+
+}
+
+template <int dim>
+SideCondition<dim>::~SideCondition ()
+{
+}
+
+template <int dim>
+void SideCondition<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+ std::vector<Vector<double> > &value_list) const
+{
+ const unsigned int n_points = points.size();
+
+ Assert (value_list.size() == n_points,
+ ExcDimensionMismatch (value_list.size(), n_points));
+
+ for (unsigned int p=0; p<n_points; ++p)
+ SideCondition<dim>::vector_value (points[p],
+ value_list[p]);
+}
+ // <h4>Conservation Law class</h4>
+ // Here we define a Conservation Law class that helps group
+ // operations and data for our Euler equations into a manageable
+ // entity. Functions will be described as their definitions appear.
+template <int dim>
+class ConsLaw
+{
+ public:
+ ConsLaw ();
+ ~ConsLaw ();
+
+ void run ();
+ void declare_parameters();
+ void load_parameters(const char *);
+
+ private:
+ void build_fe();
+ void setup_system ();
+ void initialize_system ();
+ void assemble_system (double &res_norm);
+ void solve (Vector<double> &solution, int &, double &);
+ void refine_grid ();
+ void output_results (const unsigned int cycle) const;
+ void initialize();
+ void zero_matrix();
+ void estimate();
+ void postprocess();
+ void compute_predictor();
+
+ Triangulation<dim> triangulation;
+ const MappingQ1<dim> mapping;
+
+
+ FESystem<dim> *fe_ptr;
+
+ DoFHandler<dim> dof_handler;
+
+ SparsityPattern sparsity_pattern;
+ const QGauss<dim> quadrature;
+ const QGauss<dim-1> face_quadrature;
+
+ // The actual solution to the Euler equation
+ Vector<double> solution;
+ // The current value of the solution during the Newton iteration
+ Vector<double> nlsolution;
+ // An estimate of the next time value; used for adaptivity and as a
+ // guess for the next Newton iteration.
+ Vector<double> predictor;
+ // Values after post-processing (used to output the physical variables).
+ Vector<double> ppsolution;
+ // The solution to the linear problem during the Newton iteration
+ Vector<double> dsolution;
+ Vector<double> right_hand_side;
+
+ public:
+
+ void assemble_cell_term(const FEValues<dim>& fe_v,
+ std::vector<unsigned int> &dofs,
+ unsigned int cell_no
+ );
+
+ void assemble_face_term(
+ int face_no,
+ const FEFaceValuesBase<dim>& fe_v,
+ const FEFaceValuesBase<dim>& fe_v_neighbor,
+ std::vector<unsigned int> &dofs,
+ std::vector<unsigned int> &dofs_neighbor,
+ int boundary = -1
+ );
+
+ unsigned int get_n_components() const { return N_COMP;}
+
+ private:
+ // T = current time, dT = time step, TF = final time.
+ double T, dT, TF;
+ double face_diameter;
+ double cell_diameter;
+ // An object to handle parsing the input deck.
+ ParameterHandler prm;
+ // Name of the mesh to read in.
+ string mesh;
+ InitialCondition<dim> ic;
+
+ // Enums for the various supported boundary conditions.
+ typedef enum {INFLOW_BC = 1, OUTFLOW_BC=2, NO_PENETRATION_BC=3, PRESSURE_BC=4} bc_type;
+
+ // For each boundary we store a map from boundary # to the type
+ // of boundary condition. If the boundary condition is prescribed,
+ // we store a pointer to a function object that will hold the expression
+ // for that boundary condition.
+ typedef typename std::map<unsigned int, std::pair<std::vector<bc_type>, Function<dim>*> > bdry_map_type;
+ bdry_map_type bdry_map;
+
+ void add_boundary(unsigned int bd, std::vector<bc_type>& flags, Function<dim> *bf);
+
+ // An object to store parameter information about the Aztec solver.
+ typedef struct {
+ int LIN_OUTPUT;
+ typedef enum { GMRES = 0, DIRECT = 1} solver_type;
+ solver_type SOLVER;
+ typedef enum { QUIET = 0, VERBOSE = 1 } output_type;
+ output_type OUTPUT;
+ // Linear residual tolerance.
+ double RES;
+ int MAX_ITERS;
+ // We use the ILUT preconditioner. This is similar
+ // to the ILU. FILL is the number of extra entries
+ // to add when forming the ILU decomposition.
+ double ILUT_FILL;
+ // When forming the preconditioner, for certain problems
+ // bad conditioning (or just bad luck) can cause the
+ // preconditioner to be very poorly conditioned. Hence
+ // it can help to add diagonal perturbations to the
+ // original matrix and form the preconditioner for this
+ // slightly better matrix. ATOL is an absolute perturbation
+ // that is added to the diagonal before forming the
+ // prec, and RTOL is a scaling factor $rtol >= 1$.
+ double ILUT_ATOL;
+ double ILUT_RTOL;
+ // The ILUT will drop any values that have magnitude less
+ // than this value. This is a way to
+ // manage the amount of memory used by this preconditioner.
+ double ILUT_DROP;
+ } solver_params_type;
+
+ solver_params_type solver_params;
+
+ // Some refinement parameters.
+ typedef struct {
+ typedef enum { NONE = 0, FIXED_NUMBER = 1, SHOCK = 2} refine_type;
+ double high_frac;
+ double low_frac;
+ refine_type refine;
+ double high_frac_sav;
+ double max_cells;
+ double shock_val;
+ double shock_levels;
+ } refinement_params_type;
+
+ refinement_params_type refinement_params;
+
+ // The user can set the stabilization parameter $\alpha$
+ // in the Lax-Friedrich's flux.
+ typedef struct {
+ typedef enum {CONSTANT=1,MESH=2} LF_stab_type;
+ LF_stab_type LF_stab;
+ double LF_stab_value;
+ } flux_params_type;
+
+ flux_params_type flux_params;
+
+ bool is_stationary;
+
+ // Power for the mesh stabilization term.
+ double diffusion_power;
+ double gravity;
+ // If true, we output the squared gradient of the
+ // density instead of density. Using this one can
+ // create shock plots.
+ bool schlieren_plot;
+ // How often to create an output file.
+ double output_step;
+
+ Epetra_CrsMatrix *Matrix;
+ Epetra_Map *Map;
+ Vector<double> indicator;
+
+ // Crank-Nicolson value
+ const double theta;
+
+};
+
+
+ // Asign a row of the conservation law a specified
+ // boundary type and (possibly) function.
+template <int dim>
+void ConsLaw<dim>::add_boundary(unsigned int bd,
+ std::vector<bc_type> &flags, Function<dim> *bf) {
+
+ std::pair<std::vector<bc_type>, Function<dim> *> entry(flags, bf);
+ bdry_map[bd] = entry;
+}
+
+
+ // Apply the initialial condition. Simultaneously
+ // initialize the non-linear solution.
+template <int dim>
+void ConsLaw<dim>::initialize() {
+ VectorTools::interpolate(dof_handler,
+ ic, solution);
+ VectorTools::interpolate(dof_handler,
+ ic, nlsolution);
+
+}
+
+ // <h4>Assembly</h4>
+ // <h5>Function: assemble_cell_term</h5>
+ //
+ // Assembles the cell term, adding minus the residual
+ // to the right hand side, and adding in the Jacobian
+ // contributions.
+template <int dim>
+void ConsLaw<dim>::assemble_cell_term(
+ const FEValues<dim> &fe_v,
+ std::vector<unsigned int> &dofs,
+ unsigned int cell_no
+ )
+{
+ // 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.
+ fad_double F_i;
+ int dofs_per_cell = fe_v.dofs_per_cell;
+ int n_q_points = fe_v.n_quadrature_points;
+
+ // We will define the dofs on this cell in these fad variables.
+ std::vector<fad_double> DOF(dofs_per_cell);
+
+ // Values of the conservative variables at the quadrature points.
+ std::vector<std::vector<fad_double > > W (n_q_points,
+ std::vector<fad_double >(get_n_components()));
+
+ // Values at the last time step of the conservative variables.
+ // Note that these do not use fad variables, since they do
+ // not depend on the 'variables to be sought'=DOFS.
+ std::vector<std::vector<double > > Wl (n_q_points,
+ std::vector<double >(get_n_components()));
+
+ // Here we will hold the averaged values of the conservative
+ // variables that we will linearize around (cn=Crank Nicholson).
+ std::vector<std::vector<fad_double > > Wcn (n_q_points,
+ std::vector<fad_double >(get_n_components()));
+
+ // Gradients of the current variables. It is a
+ // bit of a shame that we have to compute these; we almost don't.
+ // The nice thing about a simple conservation law is that the
+ // the flux doesn't generally involve any gradients. We do
+ // need these, however, for the diffusion stabilization.
+ std::vector<std::vector<std::vector<fad_double> > > Wgrads (n_q_points,
+ std::vector<std::vector<fad_double> >(get_n_components(),
+ std::vector<fad_double>(dim)));
+
+
+ const std::vector<double> &JxW = fe_v.get_JxW_values ();
+
+
+ // Here is the magical point where we declare a subset
+ // of the fad variables as degrees of freedom. All
+ // calculations that reference these variables (either
+ // directly or indirectly) will accumulate sensitivies
+ // with respect to these dofs.
+ for (int in = 0; in < dofs_per_cell; in++) {
+ DOF[in] = nlsolution(dofs[in]);
+ DOF[in].diff(in, dofs_per_cell);
+ }
+
+ // Here we compute the shape function values and gradients
+ // at the quadrature points. Ideally, we could call into
+ // something like get_function_values, get_function_grads,
+ // but since we don't want to make the entire solution vector
+ // fad types, only the local cell variables, we explicitly
+ // code this loop;
+ for (int q = 0; q < n_q_points; q++) {
+ for (int di = 0; di < get_n_components(); di++) {
+ W[q][di] = 0;
+ Wl[q][di] = 0;
+ Wcn[q][di] = 0;
+ for (int d = 0; d < dim; d++) {
+ Wgrads[q][di][d] = 0;
+ }
+ }
+ for (int sf = 0; sf < dofs_per_cell; sf++) {
+ int di = fe_v.get_fe().system_to_component_index(sf).first;
+ W[q][di] +=
+ DOF[sf]*fe_v.shape_value_component(sf, q, di);
+ Wl[q][di] +=
+ solution(dofs[sf])*fe_v.shape_value_component(sf, q, di);
+ Wcn[q][di] +=
+ (theta*DOF[sf]+(1-theta)*solution(dofs[sf]))*fe_v.shape_value_component(sf, q, di);
+
+ for (int d = 0; d < dim; d++) {
+ Wgrads[q][di][d] += DOF[sf]*
+ fe_v.shape_grad_component(sf, q, di)[d];
+ } // for d
+
+ }
+
+ } // for q
+
+ // Gather the flux values for all components at
+ // all of the quadrature points. This also
+ // computes the matrix of sensitivities. Perhaps
+ // this could be done in a better way, since this
+ // could be a rather large object, but for now it
+ // seems to work just fine.
+ std::vector<std::vector<std::vector<fad_double> > > flux(n_q_points,
+ std::vector<std::vector<fad_double> >(get_n_components(),
+ std::vector<fad_double>(dim, 0)));
+
+ for (unsigned int q=0; q < n_q_points; ++q) {
+ Flux<fad_double, dim>(flux[q], fe_v.get_quadrature_points()[q], Wcn[q]);
+ }
+
+ // We now have all of the function values/grads/fluxes,
+ // so perform the assembly. We have an outer loop
+ // through the components of the system, and an
+ // inner loop over the quadrature points, where we
+ // accumulate contributions to the ith residual.
+ //
+ // We initialy sum all contributions of the residual
+ // in the positive sense, so that we don't need to
+ // negative the Jacobian entries. Then, when we sum
+ // into the <code> right_hand_side </code> vector,
+ // we negate this residual.
+ for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
+ {
+ // Find which component this dof contributes to.
+ const unsigned int
+ component_i = fe_v.get_fe().system_to_component_index(i).first;
+
+ // Initialize the fad residual to zero (removes
+ // any previous sensitivities.
+ F_i = 0;
+
+ // Loop quadrature points.
+ for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point) {
+
+ fad_double fdotgv = 0;
+
+ // Integrate the flux times gradient of the test function
+ for (int d = 0; d < dim; d++)
+ fdotgv += flux[point][component_i][d]*fe_v.shape_grad_component(i, point, component_i)[d];
+
+ F_i -= fdotgv*JxW[point];
+
+ // The mass term (if the simulation is non-stationary).
+ fad_double delta_t= 1.0/dT*(W[point][component_i] - Wl[point][component_i]);
+ if (!is_stationary) F_i += delta_t*
+ fe_v.shape_value_component(i, point, component_i)*JxW[point];
+
+ // Stabilization (cell wise diffusion)
+ fad_double guv = 0;
+ for (int d = 0; d < dim; d++) {
+ guv += fe_v.shape_grad_component(i, point, component_i)[d]*
+ Wgrads[point][component_i][d];
+ }
+
+ F_i += 1.0*std::pow(cell_diameter, diffusion_power)*guv*JxW[point];
+
+ // The gravity component only enters into the energy
+ // equation and into the vertical component of the
+ // velocity.
+ if (component_i == dim - 1) {
+ F_i += gravity*Wcn[point][DENS_IDX]*fe_v.shape_value_component(i,point, component_i)*JxW[point];
+ } else if (component_i == ENERGY_IDX) {
+ F_i += gravity*Wcn[point][DENS_IDX]*Wcn[point][dim-1]*
+ fe_v.shape_value_component(i,point, component_i)*JxW[point];
+ }
+ } // for q
+
+ // Here we gain access to the array of sensitivities
+ // of the residual. We then sum these into the
+ // Epetra matrix.
+ double *values = &(F_i.fastAccessDx(0));
+ Matrix->SumIntoGlobalValues(dofs[i],
+ dofs_per_cell, &values[0], reinterpret_cast<int*>(&dofs[0]));
+
+ // Add minus the residual to the right hand side.
+ right_hand_side(dofs[i]) -= F_i.val();
+
+ } // for i
+}
+ // <h5>Function: assemble_face_term</h5>
+ // 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.
+template <int dim>
+void ConsLaw<dim>::assemble_face_term(
+ int face_no,
+ const FEFaceValuesBase<dim>& fe_v,
+ const FEFaceValuesBase<dim>& fe_v_neighbor,
+ std::vector<unsigned int> &dofs,
+ std::vector<unsigned int> &dofs_neighbor,
+ int boundary
+ )
+{
+ fad_double F_i;
+ 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<fad_double> DOF(dofs_per_cell+ndofs_per_cell);
+
+ // The conservative variables for this cell,
+ // and for
+ std::vector<std::vector<fad_double > > Wplus (n_q_points,
+ std::vector<fad_double >(get_n_components()));
+ std::vector<std::vector<fad_double > > Wminus (n_q_points,
+ std::vector<fad_double >(get_n_components()));
+
+
+ const std::vector<double> &JxW = fe_v.get_JxW_values ();
+ const std::vector<Point<dim> > &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 = (boundary < 0 ? dofs_per_cell + ndofs_per_cell : dofs_per_cell);
+ // Set the local DOFS.
+ for (int in = 0; in < dofs_per_cell; in++) {
+ DOF[in] = nlsolution(dofs[in]);
+ DOF[in].diff(in, ndofs);
+ }
+ // If present, set the neighbor dofs.
+ if (boundary < 0)
+ for (int in = 0; in < ndofs_per_cell; in++) {
+ DOF[in+dofs_per_cell] = nlsolution(dofs_neighbor[in]);
+ DOF[in+dofs_per_cell].diff(in+dofs_per_cell, ndofs);
+ }
+
+ // Set the values of the local conservative variables.
+ // Initialize all variables to zero.
+ for (int q = 0; q < n_q_points; q++) {
+ for (int di = 0; di < get_n_components(); di++) {
+ Wplus[q][di] = 0;
+ Wminus[q][di] = 0;
+ }
+ for (int sf = 0; sf < dofs_per_cell; sf++) {
+ int di = fe_v.get_fe().system_to_component_index(sf).first;
+ Wplus[q][di] +=
+ (theta*DOF[sf]+(1.0-theta)*solution(dofs[sf]))*fe_v.shape_value_component(sf, q, di);
+ }
+
+
+ // If there is a cell across, then initialize
+ // the exterior trace as a function of the other
+ // cell degrees of freedom.
+ if (boundary < 0) {
+ for (int sf = 0; sf < ndofs_per_cell; sf++) {
+ int di = fe_v_neighbor.get_fe().system_to_component_index(sf).first;
+ Wminus[q][di] +=
+ (theta*DOF[sf+dofs_per_cell]+(1.0-theta)*solution(dofs_neighbor[sf]))*
+ fe_v_neighbor.shape_value_component(sf, q, di);
+ }
+ }
+ } // for q
+
+ // If this is a boundary, then the values of $W^-$ will
+ // be either functions of $W^+$, or they will be prescribed.
+ // This switch sets them appropriately. Since we are
+ // using fad variables here, sensitivities will be updated
+ // appropriately. These sensitivities would be tremendously
+ // difficult to manage without fad!!!
+ if (boundary >= 0) {
+ // Get the boundary descriptor.
+ typename bdry_map_type::iterator bme = bdry_map.find(boundary);
+ assert(bme != bdry_map.end());
+
+ // 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<Vector<double> > bvals(n_q_points, Vector<double>(N_COMP));
+ bme->second.second->vector_value_list(fe_v.get_quadrature_points(), bvals);
+
+ // We loop the quadrature points, and we treat each
+ // component individualy.
+ for (int q = 0; q < n_q_points; q++) {
+ for (int di = 0; di < get_n_components(); di++) {
+
+ // An inflow/dirichlet type of boundary condition
+ if (bme->second.first[di] == INFLOW_BC) {
+ Wminus[q][di] = bvals[q](di);
+ } else if (bme->second.first[di] == PRESSURE_BC) {
+ // 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.
+ fad_double rho_vel_sqr = 0;
+ fad_double dens;
+
+ dens = bme->second.first[DENS_IDX] == INFLOW_BC ? bvals[q](DENS_IDX) :
+ Wplus[q][DENS_IDX];
+
+ for (int d=0; d < dim; d++) {
+ if (bme->second.first[d] == INFLOW_BC)
+ 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)/(GAMMA-1.0) +
+ 0.5*rho_vel_sqr;
+
+ } else if (bme->second.first[di] == OUTFLOW_BC) {
+ // 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.
+ fad_double vdotn = 0;
+ for (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
+
+ // Determine the Lax-Friedrich's stability parameter,
+ // and evaluate the numerical flux function at the quadrature points
+ std::vector<std::vector<fad_double> > nflux(n_q_points, std::vector<fad_double>(get_n_components(), 0));
+ double alpha = 1;
+
+ switch(flux_params.LF_stab) {
+ case flux_params_type::CONSTANT:
+ alpha = flux_params.LF_stab_value;
+ break;
+ case flux_params_type::MESH:
+ alpha = face_diameter/(2.0*dT);
+ break;
+ }
+
+ LFNumFlux<fad_double, dim>(nflux, fe_v.get_quadrature_points(), normals, Wplus, Wminus,
+ alpha);
+
+ // Now assemble the face term
+ for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) {
+ if (!fe_v.get_fe().has_support_on_face(i, face_no)) continue;
+ F_i = 0;
+ for (unsigned int point=0; point<n_q_points; ++point)
+ {
+ const unsigned int
+ component_i = fe_v.get_fe().system_to_component_index(i).first;
+
+ F_i += nflux[point][component_i]*fe_v.shape_value_component(i, point, component_i)*JxW[point];
+
+ }
+
+ // Retrieve a pointer to the jacobian.
+ double *values = &(F_i.fastAccessDx(0));
+
+ // Honestly, I forget why this can happen, but
+ // for some reason it can!!
+ if (!values) continue;
+
+ // Update the matrix. Depending on whether there
+ // is/isn't a neighboring cell, we add more/less
+ // entries.
+ Matrix->SumIntoGlobalValues(dofs[i],
+ dofs_per_cell, &values[0], reinterpret_cast<int*>(&dofs[0]));
+ if (boundary < 0) {
+ Matrix->SumIntoGlobalValues(dofs[i],
+ dofs_per_cell, &values[dofs_per_cell], reinterpret_cast<int*>(&dofs_neighbor[0]));
+ }
+
+ // And add into the residual
+ right_hand_side(dofs[i]) -= F_i.val();
+ }
+
+}
+ // <h5>Assembling the whole system</h5>
+ // Now we put all of the assembly pieces together
+ // in a routine that dispatches the correct
+ // piece for each cell/face. We keep track of
+ // the norm of the resdual for the Newton iteration.
+template <int dim>
+void ConsLaw<dim>::assemble_system (double &res_norm)
+{
+ FESystem<dim> &fe = *fe_ptr;
+ const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+
+ // We track the dofs on this cell and (if necessary)
+ // the adjacent cell.
+ std::vector<unsigned int> dofs (dofs_per_cell);
+ std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
+
+ // First we create the
+ // ``UpdateFlags'' for the
+ // ``FEValues'' and the
+ // ``FEFaceValues'' objects.
+ UpdateFlags update_flags = update_values
+ | update_gradients
+ | update_q_points
+ | update_JxW_values;
+
+ // Note, that on faces we do not
+ // need gradients but we need
+ // normal vectors.
+ UpdateFlags face_update_flags = update_values
+ | update_q_points
+ | update_JxW_values
+ | update_normal_vectors;
+
+ // On the neighboring cell we only
+ // need the shape values. Given a
+ // specific face, the quadrature
+ // points and `JxW values' are the
+ // same as for the current cells,
+ // the normal vectors are known to
+ // be the negative of the normal
+ // vectors of the current cell.
+ UpdateFlags neighbor_face_update_flags = update_values;
+
+ // Then we create the ``FEValues''
+ // object. Note, that since version
+ // 3.2.0 of deal.II the constructor
+ // of this class takes a
+ // ``Mapping'' object as first
+ // argument. Although the
+ // constructor without ``Mapping''
+ // argument is still supported it
+ // is recommended to use the new
+ // constructor. This reduces the
+ // effect of `hidden magic' (the
+ // old constructor implicitely
+ // assumes a ``MappingQ1'' mapping)
+ // and makes it easier to change
+ // the mapping object later.
+ FEValues<dim> fe_v (
+ mapping, fe, quadrature, update_flags);
+
+ // Similarly we create the
+ // ``FEFaceValues'' and
+ // ``FESubfaceValues'' objects for
+ // both, the current and the
+ // neighboring cell. Within the
+ // following nested loop over all
+ // cells and all faces of the cell
+ // they will be reinited to the
+ // current cell and the face (and
+ // subface) number.
+ FEFaceValues<dim> fe_v_face (
+ mapping, fe, face_quadrature, face_update_flags);
+ FESubfaceValues<dim> fe_v_subface (
+ mapping, fe, face_quadrature, face_update_flags);
+ FEFaceValues<dim> fe_v_face_neighbor (
+ mapping, fe, face_quadrature, neighbor_face_update_flags);
+ FESubfaceValues<dim> fe_v_subface_neighbor (
+ mapping, fe, face_quadrature, neighbor_face_update_flags);
+
+ // Furthermore we need some cell
+ // iterators.
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+ // Now we start the loop over all
+ // active cells.
+ int fdofs_per_cell = fe_v.dofs_per_cell;
+ int fn_q_points = face_quadrature.n_quadrature_points;
+
+ unsigned int cell_no = 0;
+ for (;cell!=endc; ++cell, ++cell_no)
+ {
+
+ // Now we reinit the ``FEValues''
+ // object for the current cell
+ fe_v.reinit (cell);
+
+ // Collect the local dofs and
+ // asssemble the cell term.
+ cell->get_dof_indices (dofs);
+
+ cell_diameter = cell->diameter();
+
+ assemble_cell_term(fe_v,
+ dofs,
+ cell_no);
+
+ // We use the DG style loop through faces
+ // to determine if we need to apply a
+ // 'hanging node' flux calculation or a boundary
+ // computation.
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ {
+ // First we set the face
+ // iterator
+ typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
+ face_diameter = face->diameter();
+
+ if (face->at_boundary())
+ {
+ // We reinit the
+ // ``FEFaceValues''
+ // object to the
+ // current face
+ fe_v_face.reinit (cell, face_no);
+
+ // and assemble the
+ // corresponding face
+ // terms. We send the same
+ // fe_v and dofs as described
+ // in the assembly routine.
+ assemble_face_term(
+ face_no, fe_v_face,
+ fe_v_face,
+ dofs,
+ dofs,
+ face->boundary_indicator());
+ }
+ else
+ {
+ // Now we are not on
+ // the boundary of the
+ // domain, therefore
+ // there must exist a
+ // neighboring cell.
+ typename DoFHandler<dim>::cell_iterator neighbor=
+ cell->neighbor(face_no);;
+
+ if (face->has_children())
+ {
+ // case I: This cell refined compared to neighbor
+
+ const unsigned int neighbor2=
+ cell->neighbor_of_neighbor(face_no);
+
+
+ // We loop over
+ // subfaces
+ for (unsigned int subface_no=0;
+ subface_no<GeometryInfo<dim>::subfaces_per_face;
+ ++subface_no)
+ {
+ typename DoFHandler<dim>::active_cell_iterator
+ neighbor_child
+ = cell->neighbor_child_on_subface (face_no, subface_no);
+
+ face_diameter = neighbor_child->diameter(); // working on subface
+
+ Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
+ ExcInternalError());
+ Assert (!neighbor_child->has_children(), ExcInternalError());
+
+ fe_v_subface.reinit (cell, face_no, subface_no);
+ fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
+ neighbor_child->get_dof_indices (dofs_neighbor);
+
+ // Assemble as if we are working with
+ // a DG element.
+ assemble_face_term(
+ face_no, fe_v_subface,
+ fe_v_face_neighbor,
+ dofs,
+ dofs_neighbor);
+
+ }
+ // End of ``if
+ // (face->has_children())''
+ }
+ else
+ {
+ // We have no children, but
+ // the neighbor cell may be refine
+ // compared to use
+ neighbor->get_dof_indices (dofs_neighbor);
+ if (neighbor->level() != cell->level())
+ {
+ // case II: This is refined compared to neighbor
+ Assert(neighbor->level() < cell->level(), ExcInternalError());
+ const std::pair<unsigned int, unsigned int> faceno_subfaceno=
+ cell->neighbor_of_coarser_neighbor(face_no);
+ const unsigned int neighbor_face_no=faceno_subfaceno.first,
+ neighbor_subface_no=faceno_subfaceno.second;
+
+ Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,
+ neighbor_subface_no)
+ == cell,
+ ExcInternalError());
+
+ // Reinit the
+ // appropriate
+ // ``FEFaceValues''
+ // and assemble
+ // the face
+ // terms.
+ fe_v_face.reinit (cell, face_no);
+ fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no,
+ neighbor_subface_no);
+
+ assemble_face_term(
+ face_no, fe_v_face,
+ fe_v_subface_neighbor,
+ dofs,
+ dofs_neighbor);
+
+ }
+
+ }
+ // End of ``face not at boundary'':
+ }
+ // End of loop over all faces:
+ }
+
+ // End iteration through cells.
+ }
+
+ // Notify Epetra that the matrix is done.
+ Matrix->FillComplete();
+
+ // Compute the nonlinear residual.
+ res_norm = right_hand_side.l2_norm();
+
+}
+
+ // Create a conservation law with some defaults.
+template <int dim>
+ConsLaw<dim>::ConsLaw ()
+ :
+ mapping (),
+ fe_ptr(NULL),
+ dof_handler (triangulation),
+ quadrature (2),
+ face_quadrature (2),
+ T(0),
+ dT(0.05),
+ TF(10),
+ Map(NULL),
+ Matrix(NULL),
+ is_stationary(false),
+ theta(0.5)
+{}
+
+ // At one time this example could work for both DG and
+ // continuous finite elements. The choice was made here.
+template <int dim>
+void ConsLaw<dim>::build_fe() {
+ fe_ptr = new FESystem<dim>(FE_Q<dim>(1), N_COMP);
+}
+
+ // Bye bye Conservation law.
+template <int dim>
+ConsLaw<dim>::~ConsLaw ()
+{
+ dof_handler.clear ();
+ delete fe_ptr;
+}
+
+ // <h4>Initialize System</h4>
+ // Sizes all of the vectors and sets up the
+ // sparsity patter. This function is called at
+ // the very beginning of a simulation. The function
+ // <code> setup_system </code> repeats some of these
+ // chores and is called after adaptivity in leiu
+ // of this function.
+template <int dim>
+void ConsLaw<dim>::initialize_system ()
+{
+ // First we need to distribute the
+ // DoFs.
+ dof_handler.clear();
+ dof_handler.distribute_dofs (*fe_ptr);
+
+ // Size all of the fields.
+ solution.reinit (dof_handler.n_dofs());
+ nlsolution.reinit (dof_handler.n_dofs());
+ predictor.reinit (dof_handler.n_dofs());
+ ppsolution.reinit (dof_handler.n_dofs());
+ dsolution.reinit (dof_handler.n_dofs());
+ right_hand_side.reinit (dof_handler.n_dofs());
+ indicator.reinit(triangulation.n_active_cells());
+}
+
+ // <h4>Setup System</h4>
+ // We call this function to build the sparsity
+ // and the matrix.
+template <int dim>
+void ConsLaw<dim>::setup_system ()
+{
+
+ // The DoFs of a cell are coupled
+ // with all DoFs of all neighboring
+ // cells. Therefore the maximum
+ // number of matrix entries per row
+ // is needed when all neighbors of
+ // a cell are once more refined
+ // than the cell under
+ // consideration.
+ sparsity_pattern.reinit (dof_handler.n_dofs(),
+ dof_handler.n_dofs(),
+ (GeometryInfo<dim>::faces_per_cell
+ *GeometryInfo<dim>::subfaces_per_face+1)*fe_ptr->dofs_per_cell);
+
+ // Since the continuous sparsity pattern is
+ // a subset of the DG one, and since we need
+ // the DG terms for handling hanging nodes, we use
+ // the flux pattern.
+ DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern);
+
+ sparsity_pattern.compress();
+
+ // Rebuild the map. In serial this doesn't do much,
+ // but is needed. In parallel, this would desribe
+ // the parallel dof layout.
+ if (Map) delete Map;
+ Map = new Epetra_Map(dof_handler.n_dofs(), 0, *Comm);
+
+ // Epetra can build a more efficient matrix if
+ // one knows ahead of time the maxiumum number of
+ // columns in any row entry. We traverse the sparsity
+ // to discover this.
+ int cur_row = 0;
+ int cur_col = 0;
+ int max_df = -1;
+ for (SparsityPattern::iterator s_i = sparsity_pattern.begin();
+ s_i != sparsity_pattern.end(); s_i++) {
+ if (s_i->row() != cur_row) {
+ cur_col = 0;
+ cur_row = s_i->row();
+ }
+ cur_col++;
+ if (cur_col >= max_df) max_df = cur_col;
+ }
+
+ if (cur_col >= max_df) max_df = cur_col;
+ std::cout << "max_df:" << max_df << std::endl;
+
+ // Now we build the matrix, using the constructor
+ // that optimizes with the <code> max_df </code> variable.
+ if (Matrix) delete Matrix;
+ Matrix = new Epetra_CrsMatrix(Copy, *Map, max_df+1, true);
+
+ // We add the sparsity pattern to the matrix by
+ // inserting zeros.
+ std::vector<double> vals(max_df, 0);
+ std::vector<int> row_indices(max_df);
+
+ cur_row = 0;
+ cur_col = 0;
+ for (SparsityPattern::iterator s_i = sparsity_pattern.begin();
+ s_i != sparsity_pattern.end(); s_i++) {
+ if (s_i->row() != cur_row) {
+ Matrix->InsertGlobalValues(cur_row, cur_col, &vals[0], &row_indices[0]);
+ cur_col = 0;
+ cur_row = s_i->row();
+ }
+ row_indices[cur_col++] = s_i->column();
+ }
+ // The last row.
+ Matrix->InsertGlobalValues(cur_row, cur_col, &vals[0], &row_indices[0]);
+
+ // Epetra requires this function after building or
+ // filling a matrix. It typically does some parallel
+ // bookeeping; perhaps more.
+ Matrix->FillComplete();
+
+}
+
+ // <h4>Solving the linear system</h4>
+ // Actually solve the linear system, using either
+ // Aztec of Amesos.
+template <int dim>
+void ConsLaw<dim>::solve (Vector<double> &dsolution, int &niter, double &lin_residual)
+{
+
+ // We must hand the solvers Epetra vectors.
+ // Luckily, they support the concept of a
+ // 'view', so we just send in a pointer to our
+ // dealii vectors.
+ Epetra_Vector x(View, *Map, dsolution.begin());
+ Epetra_Vector b(View, *Map, right_hand_side.begin());
+
+ // The Direct option selects the Amesos solver.
+ if (solver_params.SOLVER == solver_params_type::DIRECT) {
+
+ // Setup for solving with
+ // Amesos.
+ Epetra_LinearProblem prob;
+ prob.SetOperator(Matrix);
+ Amesos_BaseSolver *solver;
+ Amesos Factory;
+
+ // Other solvers are available
+ // and may be selected by changing this
+ // string.
+ char *stype = "Amesos_Klu";
+
+ solver = Factory.Create(stype, prob);
+
+ Assert (solver != NULL, ExcInternalError());
+
+ // There are two parts to the direct solve.
+ // As I understand, the symbolic part figures
+ // out the sparsity patterns, and then the
+ // numerical part actually performs Gaussian
+ // elimination or whatever the approach is.
+ if (solver_params.OUTPUT == solver_params_type::VERBOSE)
+ std::cout << "Starting Symbolic fact\n" << std::flush;
+
+ solver->SymbolicFactorization();
+
+ if (solver_params.OUTPUT == solver_params_type::VERBOSE)
+ std::cout << "Starting Numeric fact\n" << std::flush;
+
+ solver->NumericFactorization();
+
+
+ // Define the linear problem by setting the
+ // right hand and left hand sides.
+ prob.SetRHS(&b);
+ prob.SetLHS(&x);
+ // And finally solve the problem.
+ if (solver_params.OUTPUT == solver_params_type::VERBOSE)
+ std::cout << "Starting solve\n" << std::flush;
+ solver->Solve();
+ niter = 0;
+ lin_residual = 0;
+
+ // We must free the solver that was created
+ // for us.
+ delete solver;
+
+ } else if (solver_params.SOLVER == solver_params_type::GMRES) {
+
+ // For the iterative solvers, we use Aztec.
+ AztecOO Solver;
+
+ // Select the appropriate level of verbosity.
+ if (solver_params.OUTPUT == solver_params_type::QUIET)
+ Solver.SetAztecOption(AZ_output, AZ_none);
+
+ if (solver_params.OUTPUT == solver_params_type::VERBOSE)
+ Solver.SetAztecOption(AZ_output, AZ_all);
+
+ // Select gmres. Other solvers are available.
+ Solver.SetAztecOption(AZ_solver, AZ_gmres);
+ Solver.SetRHS(&b);
+ Solver.SetLHS(&x);
+
+ // Set up the ILUT preconditioner. I do not know
+ // why, but we must pretend like we are in parallel
+ // using domain decomposition or the preconditioner
+ // refuses to activate.
+ Solver.SetAztecOption(AZ_precond, AZ_dom_decomp);
+ Solver.SetAztecOption(AZ_subdomain_solve, AZ_ilut);
+ Solver.SetAztecOption(AZ_overlap, 0);
+ Solver.SetAztecOption(AZ_reorder, 0);
+
+ // ILUT parameters as described above.
+ Solver.SetAztecParam(AZ_drop, solver_params.ILUT_DROP);
+ Solver.SetAztecParam(AZ_ilut_fill, solver_params.ILUT_FILL);
+ Solver.SetAztecParam(AZ_athresh, solver_params.ILUT_ATOL);
+ Solver.SetAztecParam(AZ_rthresh, solver_params.ILUT_RTOL);
+ Solver.SetUserMatrix(Matrix);
+
+ // Run the solver iteration. Collect the number
+ // of iterations and the residual.
+ Solver.Iterate(solver_params.MAX_ITERS, solver_params.RES);
+ niter = Solver.NumIters();
+ lin_residual = Solver.TrueResidual();
+ }
+}
+
+ // <h4>Postprocessing and Output</h4>
+ // Recover the physical variables from the conservative
+ // variables so that output will be (perhaps) more
+ // meaningfull.
+template <int dim>
+void ConsLaw<dim>::postprocess() {
+ const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+ std::vector<unsigned int> dofs (dofs_per_cell);
+ UpdateFlags update_flags = update_values
+ | update_gradients
+ | update_q_points
+ | update_JxW_values;
+ UpdateFlags update_flags1 = update_values
+ | update_gradients
+ | update_q_points
+ | update_JxW_values;
+
+ QGauss<dim> quadrature_formula(4);
+
+ const std::vector<Point<dim> > &us = fe_ptr->base_element(0).get_unit_support_points();
+
+
+ Quadrature<dim> unit_support(us);
+
+ int n_q_points = quadrature_formula.n_quadrature_points;
+ int n_uq_points = unit_support.n_quadrature_points;
+
+ FEValues<dim> fe_v (
+ mapping, *fe_ptr, quadrature_formula, update_flags);
+
+ FEValues<dim> fe_v_unit (
+ mapping, *fe_ptr, unit_support, update_flags1);
+
+ std::vector<Vector<double> > U(n_uq_points,
+ Vector<double>(get_n_components()));
+ std::vector<Vector<double> > UU(n_q_points,
+ Vector<double>(get_n_components()));
+ std::vector<std::vector<Tensor<1,dim> > > dU(n_uq_points,
+ std::vector<Tensor<1,dim> >(get_n_components()));
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+ // Loop the cells
+ for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) {
+ cell->get_dof_indices (dofs);
+ fe_v_unit.reinit(cell);
+ fe_v.reinit(cell);
+
+ fe_v_unit.get_function_values(solution, U);
+ fe_v_unit.get_function_grads(solution, dU);
+ fe_v.get_function_values(solution, UU);
+
+ const std::vector<double> &JxW = fe_v.get_JxW_values ();
+
+ for (int q = 0; q < fe_v.get_fe().base_element(0).n_dofs_per_cell(); q++) {
+ unsigned int didx = fe_v.get_fe().component_to_system_index(DENS_IDX, q);
+ unsigned int eidx = fe_v.get_fe().component_to_system_index(ENERGY_IDX, q);
+ double rho_normVsqr = 0;
+ for (int d = 0; d < dim; d++) {
+ unsigned int vidx = fe_v.get_fe().component_to_system_index(d, q);
+ ppsolution(dofs[vidx]) = solution(dofs[vidx])/solution(dofs[didx]);
+ rho_normVsqr += solution(dofs[vidx])*solution(dofs[vidx]);
+ }
+ rho_normVsqr /= solution(dofs[didx]);
+ // Pressure
+ ppsolution(dofs[eidx]) = (GAMMA-1.0)*(solution(dofs[eidx]) - 0.5*rho_normVsqr);
+
+ // Either output density or gradient squared of density,
+ // depending on what the user wants.
+ if (!schlieren_plot) {
+ ppsolution(dofs[didx]) = solution(dofs[didx]);
+ } else {
+ double ng = 0;
+ for (int i = 0; i < dim; i++) ng += dU[q][DENS_IDX][i]*dU[q][DENS_IDX][i];
+ ng = std::sqrt(ng);
+ ppsolution(dofs[didx]) = ng;
+ }
+ }
+
+ } // cell
+
+}
+
+ // Loop and assign a value for refinement. We
+ // simply use the density squared, which selects
+ // shocks with some success.
+template <int dim>
+void ConsLaw<dim>::estimate() {
+
+ const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+ std::vector<unsigned int> dofs (dofs_per_cell);
+ UpdateFlags update_flags = update_values
+ | update_gradients
+ | update_q_points
+ | update_JxW_values;
+
+ QGauss<dim> quadrature_formula(1);
+ int n_q_points = quadrature_formula.n_quadrature_points;
+
+
+ FEValues<dim> fe_v (
+ mapping, *fe_ptr, quadrature_formula, update_flags);
+
+ std::vector<Vector<double> > U(n_q_points,
+ Vector<double>(get_n_components()));
+ std::vector<std::vector<Tensor<1,dim> > > dU(n_q_points,
+ std::vector<Tensor<1,dim> >(get_n_components()));
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) {
+ fe_v.reinit(cell);
+
+ fe_v.get_function_values(predictor, U);
+ fe_v.get_function_grads(predictor, dU);
+
+ indicator(cell_no) = 0;
+ for (int q = 0; q < n_q_points; q++) {
+ double ng = 0;
+ for (int d = 0; d < dim; d++) ng += dU[q][DENS_IDX][d]*dU[q][DENS_IDX][d];
+
+ indicator(cell_no) += std::log(1+std::sqrt(ng));
+
+ }
+ indicator(cell_no) /= n_q_points;
+
+ }
+}
+
+template <int dim>
+void ConsLaw<dim>::refine_grid ()
+{
+
+ SolutionTransfer<dim, double> soltrans(dof_handler);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+ // Loop cells. If the indicator
+ // for the cell matches the refinement criterion,
+ // refine, else unrefine. The unrefinement has
+ // a slight hysterisis to avoid 'flashing' from refined
+ // to unrefined.
+ for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) {
+ cell->clear_coarsen_flag();
+ cell->clear_refine_flag();
+ if (cell->level() < refinement_params.shock_levels &&
+ std::fabs(indicator(cell_no)) > refinement_params.shock_val ) {
+ cell->set_refine_flag();
+ } else {
+ if (cell->level() > 0 &&
+ std::fabs(indicator(cell_no)) < 0.75*refinement_params.shock_val)
+ cell->set_coarsen_flag();
+ }
+ }
+
+ // The following code prolongs the solution
+ // to the new grid and carries out the refinement.
+ std::vector<Vector<double> > interp_in;
+ std::vector<Vector<double> > interp_out;
+
+ interp_in.push_back(solution);
+ interp_in.push_back(predictor);
+
+ triangulation.prepare_coarsening_and_refinement();
+ soltrans.prepare_for_coarsening_and_refinement(interp_in);
+
+ triangulation.execute_coarsening_and_refinement ();
+
+ dof_handler.clear();
+ dof_handler.distribute_dofs (*fe_ptr);
+
+ {
+ Vector<double> new_solution(1);
+ Vector<double> new_predictor(1);
+
+ interp_out.push_back(new_solution);
+ interp_out.push_back(new_predictor);
+ interp_out[0].reinit(dof_handler.n_dofs());
+ interp_out[1].reinit(dof_handler.n_dofs());
+ }
+
+ soltrans.interpolate(interp_in, interp_out);
+
+ // Let the vector delete a very small vector
+ solution.reinit(1);
+ predictor.reinit(1);
+ solution.swap(interp_out[0]);
+ predictor.swap(interp_out[1]);
+
+ // resize these vectors for the new grid.
+ nlsolution.reinit(dof_handler.n_dofs());
+ ppsolution.reinit(dof_handler.n_dofs());
+ nlsolution = solution;
+ dsolution.reinit (dof_handler.n_dofs());
+ right_hand_side.reinit (dof_handler.n_dofs());
+
+ indicator.reinit(triangulation.n_active_cells());
+
+}
+
+template <int dim>
+void ConsLaw<dim>::output_results (const unsigned int cycle) const
+{
+ char filename[512];
+ std::sprintf(filename, "solution-%03d.vtk", cycle);
+ std::ofstream output (filename);
+
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler (dof_handler);
+ std::vector<std::string> solution_names;
+
+ // Rename the output with the physical variable
+ // names. Send the post-processed values.
+ solution_names.clear();
+ for (int i = 0; i < dim; i++) {
+ char buf[512];
+ std::sprintf(buf, "v_%d", i);
+ solution_names.push_back (buf);
+ }
+ solution_names.push_back("density");
+ solution_names.push_back("pressure");
+ data_out.add_data_vector (ppsolution, solution_names);
+
+ data_out.add_data_vector (indicator, "error");
+ data_out.build_patches ();
+ data_out.write_vtk (output);
+
+ output.close();
+}
+
+ // <h4>Parsing the Input Deck</h4>
+ // Declare the parameters for the
+ // input deck. We assume a certain
+ // maximum number of boundaries and process
+ // any boundary the user supplies up to
+ // that maximum number. We
+ // leave a detailed explanation of these
+ // parameters to our description of the input
+ // sample file.
+const UInt MAX_BD = 10;
+template <int dim>
+void ConsLaw<dim>::declare_parameters() {
+
+ // Global scope parameters/
+ prm.declare_entry("mesh", "grid.inp",
+ Patterns::Anything(),
+ "intput file");
+
+ prm.declare_entry("diffusion power", "2.0",
+ Patterns::Double(),
+ "power of mesh size for diffusion");
+
+ prm.declare_entry("gravity", "0.0",
+ Patterns::Double(),
+ "gravity forcing");
+
+ // Time stepping block
+ prm.enter_subsection("time stepping");
+ prm.declare_entry("time step", "0.1",
+ Patterns::Double(),
+ "simulation time step");
+ prm.declare_entry("final time", "10.0",
+ Patterns::Double(),
+ "simulation end time");
+ prm.leave_subsection();
+
+
+ // Declare the boundary parameters
+ for (int b = 0; b < MAX_BD; b++) {
+ char bd[512];
+ std::sprintf(bd, "boundary_%d", b);
+ prm.enter_subsection(bd);
+ prm.declare_entry("no penetration", "false",
+ Patterns::Selection("true|false"),
+ "<true|false>");
+ // declare a slot for each of the conservative
+ // variables.
+ for (int di = 0; di < N_COMP; di++) {
+ char var[512];
+ std::sprintf(var, "w_%d", di);
+ prm.declare_entry(var, "outflow",
+ Patterns::Selection(
+ "inflow|outflow|pressure"),
+ "<inflow|outflow|pressure>");
+
+ // for dirichlet, a function in x,y,z
+ std::sprintf(var, "w_%d value", di);
+ prm.declare_entry(var, "0.0",
+ Patterns::Anything(),
+ "expression in x,y,z");
+ }
+
+ prm.leave_subsection();
+ }
+
+ // Initial condition block.
+ prm.enter_subsection("initial condition");
+ for (int di = 0; di < N_COMP; di++) {
+ char var[512];
+ std::sprintf(var, "w_%d", di);
+
+ // for dirichlet, a function in x,y,z
+ std::sprintf(var, "w_%d value", di);
+ prm.declare_entry(var, "0.0",
+ Patterns::Anything(),
+ "expression in x,y,z");
+ }
+ prm.leave_subsection();
+
+ // The linear solver block.
+ prm.enter_subsection("linear solver");
+ prm.declare_entry("output", "quiet",
+ Patterns::Selection(
+ "quiet|verbose"),
+ "<quiet|verbose>");
+ prm.declare_entry("method", "gmres",
+ Patterns::Selection(
+ "gmres|direct"),
+ "<gmres|direct>");
+ prm.declare_entry("residual", "1e-10",
+ Patterns::Double(),
+ "linear solver residual");
+ prm.declare_entry("max iters", "300",
+ Patterns::Double(),
+ "maximum solver iterations");
+ prm.declare_entry("ilut fill", "2",
+ Patterns::Double(),
+ "ilut preconditioner fill");
+ prm.declare_entry("ilut absolute tolerance", "1e-9",
+ Patterns::Double(),
+ "ilut preconditioner tolerance");
+ prm.declare_entry("ilut relative tolerance", "1.1",
+ Patterns::Double(),
+ "rel tol");
+ prm.declare_entry("ilut drop tolerance", "1e-10",
+ Patterns::Double(),
+ "ilut drop tol");
+ prm.leave_subsection();
+
+
+ // A refinement controller block.
+ prm.enter_subsection("refinement");
+ prm.declare_entry("refinement", "none",
+ Patterns::Selection(
+ "none|fixed number|shock"),
+ "<on|off>");
+ prm.declare_entry("refinement fraction", "0.1",
+ Patterns::Double(),
+ "Fraction of high refinement");
+ prm.declare_entry("unrefinement fraction", "0.1",
+ Patterns::Double(),
+ "Fraction of low unrefinement");
+ prm.declare_entry("max elements", "1000000",
+ Patterns::Double(),
+ "maximum number of elements");
+ prm.declare_entry("shock value", "4.0",
+ Patterns::Double(),
+ "value for shock indicator");
+ prm.declare_entry("shock levels", "3.0",
+ Patterns::Double(),
+ "number of shock refinement levels");
+ prm.leave_subsection();
+
+ // Output control.
+ prm.enter_subsection("output");
+ prm.declare_entry("density", "standard",
+ Patterns::Selection(
+ "standard|schlieren"),
+ "<standard|schlieren>");
+ prm.declare_entry("step", "-1",
+ Patterns::Double(),
+ "output once per this period");
+ prm.leave_subsection();
+
+ // Flux control
+ prm.enter_subsection("flux");
+ prm.declare_entry("stab", "alpha",
+ Patterns::Selection(
+ "alpha|constant|mesh"),
+ "<alpha|constant|mesh>");
+ prm.declare_entry("stab value", "1",
+ Patterns::Double(),
+ "alpha stabilization");
+ prm.leave_subsection();
+
+
+}
+
+ // Code to actually parse an input file. This function
+ // matches the declarations above.
+template <int dim>
+void ConsLaw<dim>::load_parameters(const char *infile){
+
+ prm.read_input(infile);
+
+ // The global parameters.
+ mesh = prm.get("mesh");
+
+ diffusion_power = prm.get_double("diffusion power");
+
+ gravity = prm.get_double("gravity");
+
+ // The time stepping.
+ prm.enter_subsection("time stepping");
+ dT = prm.get_double("time step");
+ std::cout << "dT=" << dT << std::endl;
+ if (dT == 0) {
+ is_stationary = true;
+ dT = 1.0;
+ TF = 1.0;
+ std::cout << "Stationary mode" << std::endl;
+ }
+ TF = prm.get_double("final time");
+ std::cout << "TF=" << TF << std::endl;
+ prm.leave_subsection();
+
+ // The boundary info
+ for (int b = 0; b < MAX_BD; b++) {
+ std::vector<bc_type> flags(N_COMP, OUTFLOW_BC);
+
+ // Define a parser for every boundary, though it may be
+ // unused.
+ SideCondition<dim> *sd = new SideCondition<dim>(N_COMP);
+ char bd[512];
+ std::sprintf(bd, "boundary_%d", b);
+ prm.enter_subsection(bd);
+
+ const std::string &nopen = prm.get("no penetration");
+
+ // Determine how each component is handled.
+ for (int di = 0; di < N_COMP; di++) {
+ char var[512];
+ std::sprintf(var, "w_%d", di);
+ std::string btype = prm.get(var);
+ std::sprintf(var, "w_%d value", di);
+ std::string var_value = prm.get(var);
+
+ if (di < dim && nopen == "true") {
+ flags[di] = NO_PENETRATION_BC;
+ } else if (btype == "inflow") {
+ flags[di] = INFLOW_BC;
+ sd->set_coeff_row(di, var_value);
+ } else if (btype == "pressure") {
+ flags[di] = PRESSURE_BC;
+ sd->set_coeff_row(di, var_value);
+ }
+ }
+ prm.leave_subsection();
+
+ // Add the boundary condition to the law.
+ sd->Init();
+ add_boundary(b, flags, sd);
+ }
+
+ // Initial conditions.
+ prm.enter_subsection("initial condition");
+ for (int di = 0; di < N_COMP; di++) {
+ char var[512];
+
+ std::sprintf(var, "w_%d value", di);
+ std::string var_value = prm.get(var);
+ ic.set_ic(di, var_value);
+ }
+ ic.Init();
+ prm.leave_subsection();
+
+ // The linear solver.
+ prm.enter_subsection("linear solver");
+ const std::string &op = prm.get("output");
+ if (op == "verbose") solver_params.OUTPUT = solver_params_type::VERBOSE;
+ if (op == "quiet") solver_params.OUTPUT = solver_params_type::QUIET;
+ const std::string &sv = prm.get("method");
+ if (sv == "direct") {
+ solver_params.SOLVER = solver_params_type::DIRECT;
+ } else if (sv == "gmres") {
+ solver_params.SOLVER = solver_params_type::GMRES;
+ }
+
+ solver_params.RES = prm.get_double("residual");
+ solver_params.MAX_ITERS = (int) prm.get_double("max iters");
+ solver_params.ILUT_FILL = prm.get_double("ilut fill");
+ solver_params.ILUT_ATOL = prm.get_double("ilut absolute tolerance");
+ solver_params.ILUT_RTOL = prm.get_double("ilut relative tolerance");
+ solver_params.ILUT_DROP = prm.get_double("ilut drop tolerance");
+ solver_params.RES = prm.get_double("residual");
+ prm.leave_subsection();
+
+
+ // And refiement.
+ prm.enter_subsection("refinement");
+ const std::string &ref = prm.get("refinement");
+ if (ref == "none") {
+ refinement_params.refine = refinement_params_type::NONE;
+ } else if (ref == "fixed number") {
+ refinement_params.refine = refinement_params_type::FIXED_NUMBER;
+ } else if (ref == "shock") {
+ refinement_params.refine = refinement_params_type::SHOCK;
+ } else
+ refinement_params.high_frac = prm.get_double("refinement fraction");
+ refinement_params.high_frac_sav = refinement_params.high_frac;
+ refinement_params.low_frac = prm.get_double("unrefinement fraction");
+ refinement_params.max_cells = prm.get_double("max elements");
+ refinement_params.shock_val = prm.get_double("shock value");
+ refinement_params.shock_levels = prm.get_double("shock levels");
+ prm.leave_subsection();
+
+ // Output control.
+ prm.enter_subsection("output");
+ const std::string &dens = prm.get("density");
+ schlieren_plot = dens == "schlieren" ? true : false;
+ output_step = prm.get_double("step");
+ prm.leave_subsection();
+
+ // Flux control.
+ prm.enter_subsection("flux");
+ const std::string &stab = prm.get("stab");
+ if (stab == "constant") {
+ flux_params.LF_stab = flux_params_type::CONSTANT;
+ } else if (stab == "mesh ") {
+ flux_params.LF_stab = flux_params_type::MESH;
+ }
+ flux_params.LF_stab_value = prm.get_double("stab value");
+ prm.leave_subsection();
+
+
+}
+
+template<int dim>
+void ConsLaw<dim>::zero_matrix() {
+ Matrix->PutScalar(0); Matrix->FillComplete();
+}
+
+ // We use a predictor to try and make adaptivity
+ // work better. The idea is to try and refine ahead
+ // of a front, rather than stepping into a coarse
+ // set of elements and smearing the solution. This
+ // simple time extrapolator does the job.
+template<int dim>
+void ConsLaw<dim>::compute_predictor() {
+ predictor = nlsolution;
+ predictor.sadd(3/2.0, -1/2.0, solution);
+}
+
+ // <h4>Run the simulation</h4>
+ // Contains the initialization
+ // the time loop, and the inner Newton iteration.
+template <int dim>
+void ConsLaw<dim>::run ()
+{
+
+ // Open and load the mesh.
+ GridIn<dim> grid_in;
+ grid_in.attach_triangulation(triangulation);
+ std::cout << "Opening mesh <" << mesh << ">" << std::endl;
+ std::ifstream input_file(mesh.c_str(), std::ios::in);
+
+ Assert (infile,
+ ExcFileNotOpen());
+
+ grid_in.read_ucd(input_file);
+ input_file.close();
+
+ build_fe();
+
+ unsigned int nstep = 0;
+
+ // Initialize fields and matrices.
+ initialize_system ();
+ setup_system();
+ initialize();
+ predictor = solution;
+
+ // Initial refinement. We apply the ic,
+ // estimate, refine, and repeat until
+ // happy.
+ if (refinement_params.refine != refinement_params_type::NONE)
+ for (int i = 0; i < refinement_params.shock_levels; i++) {
+ estimate();
+ refine_grid();
+ setup_system();
+ initialize();
+ predictor = solution;
+ }
+ postprocess();
+ output_results (nstep);
+
+ // Determine when we will output next.
+ double next_output = T + output_step;
+
+ // <h5>Main time stepping loop</h5>
+ predictor = solution;
+ while(T < TF)
+ {
+ std::cout << "T=" << T << ", ";
+
+
+ std::cout << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << std::endl;
+
+
+ std::cout << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl;
+
+ bool nonlin_done = false;
+ double res_norm;
+ int lin_iter;
+
+ // Print some relevant information during the
+ // Newton iteration.
+ std::cout << "NonLin Res: Lin Iter Lin Res" << std::endl;
+ std::cout << "______________________________________" << std::endl;
+
+ int max_nonlin = 7;
+ int nonlin_iter = 0;
+ double lin_res;
+
+ // <h6>Newton iteration</h6>
+ nlsolution = predictor;
+ while (!nonlin_done) {
+ lin_iter = 0;
+ zero_matrix();
+ right_hand_side = 0;
+ assemble_system (res_norm);
+ // Flash a star to the screen so one can
+ // know when the assembly has stopped and the linear
+ // solution is starting.
+ std::cout << "* " << std::flush;
+
+ // Test against a (hardcoded) nonlinear tolderance.
+ // Do not solve the linear system at the last step
+ // (since it would be a waste).
+
+ if (fabs(res_norm) < 1e-10) {
+ nonlin_done = true;
+ } else {
+ // Solve the linear system and update with the
+ // delta.
+ dsolution = 0;
+ solve (dsolution, lin_iter, lin_res);
+ nlsolution.add(1.0, dsolution);
+ }
+
+ // Print the residuals.
+ std::printf("%-16.3e %04d %-5.2e\n",
+ res_norm, lin_iter, lin_res);
+
+ nonlin_iter++;
+ }
+
+ // Various post convergence tasks.
+ compute_predictor();
+
+ solution = nlsolution;
+
+
+ estimate();
+
+ postprocess();
+
+ T += dT;
+
+ // Output if it is time.
+ if (output_step < 0) {
+ output_results (++nstep);
+ } else if (T >= next_output) {
+ output_results (++nstep);
+ next_output += output_step;
+ }
+
+ // Refine, if refinement is selected.
+ if (refinement_params.refine != refinement_params_type::NONE) {
+ refine_grid();
+ setup_system();
+ }
+ }
+}
+
+ // The following ``main'' function is
+ // similar to previous examples and
+ // need not to be commented on.
+int main (int argc, char *argv[])
+{
+
+ MPI_Init(&argc, &argv);
+ Comm = new Epetra_MpiComm(MPI_COMM_WORLD);
+
+ if (argc != 2) {
+ std::cout << "Usage:" << argv[0] << " infile" << std::endl;
+ std::exit(1);
+ }
+ try
+ {
+ ConsLaw<DIMENSION> cons;
+ cons.declare_parameters();
+ cons.load_parameters(argv[1]);
+ cons.run ();
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+
+ return 0;
+}
+