From 6d151e9328f17515b826e923400911199a0d087d Mon Sep 17 00:00:00 2001 From: bangerth Date: Sat, 10 May 2008 23:50:22 +0000 Subject: [PATCH] Reindent everything. git-svn-id: https://svn.dealii.org/trunk@16071 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-33/step-33.cc | 1519 +++++++++++++-------------- 1 file changed, 757 insertions(+), 762 deletions(-) diff --git a/deal.II/examples/step-33/step-33.cc b/deal.II/examples/step-33/step-33.cc index 7bd3a0e930..234d66ea64 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -11,12 +11,7 @@ /* 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. - - // @sect3{Include files} + // @sect3{Include files} // Here we have the necessary TRILINOS includes. // @@ -80,98 +75,98 @@ #include #include - // Introduce the dealii library into the current namespace. + // Introduce the dealii library into the current namespace. using namespace dealii; #define DIMENSION 2 - // We define a shorter name for the automatic differentiation - // type. + // We define a shorter name for the automatic differentiation + // type. typedef Sacado::Fad::DFad fad_double; typedef unsigned int UInt; - // The Epetra library requires a 'communicator', which describes - // the layout of a parallel (or serial) set of processors. + // The Epetra library requires a 'communicator', which describes + // the layout of a parallel (or serial) set of processors. Epetra_SerialComm *Comm; - // @sect3{Flux function definition} - // 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. + // @sect3{Flux function definition} + // 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 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. + // 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. + // 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 void Flux(std::vector > &flux, const Point &/*point*/, const std::vector &W) { - // Pressure is a dependent variable: $p = - // (\gamma - 1)(E-\frac{1}{2} \rho |v|^2)$. - number rho_normVsqr; - for (unsigned 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 (unsigned int d = 0; d < dim; d++) { - for (unsigned 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 + // Pressure is a dependent variable: $p = + // (\gamma - 1)(E-\frac{1}{2} \rho |v|^2)$. + number rho_normVsqr; + for (unsigned 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 (unsigned int d = 0; d < dim; d++) { + for (unsigned 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$. + // 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 void LFNumFlux( - std::vector > &nflux, - const std::vector > &points, - const std::vector > &normals, - const std::vector > &Wplus, - const std::vector > &Wminus, - double alpha) + std::vector > &nflux, + const std::vector > &points, + const std::vector > &normals, + const std::vector > &Wplus, + const std::vector > &Wminus, + double alpha) { const unsigned int n_q_points = points.size(); - // We evaluate the flux at each of the quadrature points. + // We evaluate the flux at each of the quadrature points. for (unsigned int q = 0; q < n_q_points; q++) { std::vector > iflux(N_COMP, - std::vector(dim, 0)); + std::vector(dim, 0)); std::vector > oflux(N_COMP, - std::vector(dim, 0)); + std::vector(dim, 0)); Flux(iflux, points[q], Wplus[q]); Flux(oflux, points[q], Wminus[q]); @@ -181,28 +176,28 @@ void LFNumFlux( for (unsigned 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]); + nflux[q][di] += 0.5*alpha*(Wplus[q][di] - Wminus[q][di]); } } } - // @sect3{Initial and side condition parsing} - // For the initial condition we use the expression parser function - // object. + // @sect3{Initial and side condition parsing} + // For the initial condition we use the expression parser function + // object. template class InitialCondition : public FunctionParser { 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. + // 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. + // 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; @@ -210,8 +205,8 @@ class InitialCondition : public FunctionParser virtual void vector_value_list (const std::vector > &points, std::vector > &value_list) const; - private: - std::vector expressions; + private: + std::vector expressions; }; template @@ -220,21 +215,21 @@ InitialCondition::InitialCondition () : 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. + // Here we set up x,y,z as the variables that one should use in the input + // deck to describe their initial condition. template void InitialCondition::Init() { - std::map constants; - constants["M_PI"] = M_PI; - std::string variables = (dim == 2 ? "x,y" : "x,y,z"); + std::map constants; + constants["M_PI"] = M_PI; + std::string variables = (dim == 2 ? "x,y" : "x,y,z"); - FunctionParser::initialize(variables, expressions, constants); + FunctionParser::initialize(variables, expressions, constants); } template void InitialCondition::vector_value_list (const std::vector > &points, - std::vector > &value_list) const + std::vector > &value_list) const { const unsigned int n_points = points.size(); @@ -243,10 +238,10 @@ void InitialCondition::vector_value_list (const std::vector > &p for (unsigned int p=0; p::vector_value (points[p], - value_list[p]); + value_list[p]); } - // As above, we use the expression function parser for boundary conditions. + // As above, we use the expression function parser for boundary conditions. template class SideCondition : public FunctionParser { @@ -254,15 +249,15 @@ class SideCondition : public FunctionParser SideCondition (int ncomp); ~SideCondition (); - // As above. + // As above. void Init(); - // As above. + // As above. void set_coeff_row(int _row_n, std::string &expr); virtual void vector_value_list (const std::vector > &points, std::vector > &value_list) const; private: - std::vector expressions; + std::vector expressions; }; template @@ -274,16 +269,16 @@ SideCondition::SideCondition (int ncomp) : template void SideCondition::set_coeff_row (int _row_n, std::string &expr) { - expressions[_row_n] = expr; + expressions[_row_n] = expr; } template void SideCondition::Init() { - std::map constants; - constants["M_PI"] = M_PI; - std::string variables = (dim == 2 ? "x,y" : "x,y,z"); + std::map constants; + constants["M_PI"] = M_PI; + std::string variables = (dim == 2 ? "x,y" : "x,y,z"); - FunctionParser::initialize(variables, expressions, constants); + FunctionParser::initialize(variables, expressions, constants); } @@ -306,9 +301,9 @@ void SideCondition::vector_value_list (const std::vector > &poin value_list[p]); } // @sect3{Conservation Law class} - // 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. + // 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 class ConsLaw { @@ -364,138 +359,138 @@ class ConsLaw void assemble_cell_term(const FEValues& fe_v, std::vector &dofs, unsigned int cell_no - ); + ); void assemble_face_term( - int face_no, - const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - std::vector &dofs, - std::vector &dofs_neighbor, - int boundary = -1 - ); + int face_no, + const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + std::vector &dofs, + std::vector &dofs_neighbor, + int boundary = -1 + ); unsigned int get_n_components() const { return N_COMP;} private: - // T = current time, dT = time step, TF = final time. + // 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. + // An object to handle parsing the input deck. ParameterHandler prm; - // Name of the mesh to read in. + // Name of the mesh to read in. string mesh; InitialCondition ic; - // Enums for the various supported boundary conditions. + // 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. + // 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, Function*> > bdry_map_type; bdry_map_type bdry_map; void add_boundary(unsigned int bd, std::vector& flags, Function *bf); - // An object to store parameter information about the Aztec solver. + // 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; + 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. + // 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; + 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. + // 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; + 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. + // 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. + // 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. + // How often to create an output file. double output_step; Epetra_Map *Map; Epetra_CrsMatrix *Matrix; Vector indicator; - // Crank-Nicolson value + // Crank-Nicolson value const double theta; }; - // Asign a row of the conservation law a specified - // boundary type and (possibly) function. + // Asign a row of the conservation law a specified + // boundary type and (possibly) function. template void ConsLaw::add_boundary(unsigned int bd, - std::vector &flags, Function *bf) { + std::vector &flags, Function *bf) { std::pair, Function *> entry(flags, bf); bdry_map[bd] = entry; } - // Apply the initialial condition. Simultaneously - // initialize the non-linear solution. + // Apply the initialial condition. Simultaneously + // initialize the non-linear solution. template void ConsLaw::initialize() { - VectorTools::interpolate(dof_handler, + VectorTools::interpolate(dof_handler, ic, solution); - VectorTools::interpolate(dof_handler, + VectorTools::interpolate(dof_handler, ic, nlsolution); } @@ -511,63 +506,63 @@ void ConsLaw::assemble_cell_term( const FEValues &fe_v, std::vector &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. + // 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; unsigned int dofs_per_cell = fe_v.dofs_per_cell; unsigned int n_q_points = fe_v.n_quadrature_points; - // We will define the dofs on this cell in these fad variables. + // We will define the dofs on this cell in these fad variables. std::vector DOF(dofs_per_cell); - // Values of the conservative variables at the quadrature points. + // Values of the conservative variables at the quadrature points. std::vector > W (n_q_points, - std::vector(get_n_components())); + std::vector(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. + // 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 > Wl (n_q_points, - std::vector(get_n_components())); + std::vector(get_n_components())); - // Here we will hold the averaged values of the conservative - // variables that we will linearize around (cn=Crank Nicholson). + // Here we will hold the averaged values of the conservative + // variables that we will linearize around (cn=Crank Nicholson). std::vector > Wcn (n_q_points, - std::vector(get_n_components())); + std::vector(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 > > Wgrads (n_q_points, - std::vector >(get_n_components(), - std::vector(dim))); + // 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 > > Wgrads (n_q_points, + std::vector >(get_n_components(), + std::vector(dim))); const std::vector &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. + // 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 (unsigned int in = 0; in < dofs_per_cell; in++) { - DOF[in] = nlsolution(dofs[in]); - DOF[in].diff(in, dofs_per_cell); + 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; + // 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 (unsigned int q = 0; q < n_q_points; q++) { for (unsigned int di = 0; di < get_n_components(); di++) { W[q][di] = 0; @@ -578,18 +573,18 @@ void ConsLaw::assemble_cell_term( } } for (unsigned 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 (unsigned int d = 0; d < dim; d++) { - Wgrads[q][di][d] += DOF[sf]* - fe_v.shape_grad_component(sf, q, di)[d]; - } // for d + 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 (unsigned int d = 0; d < dim; d++) { + Wgrads[q][di][d] += DOF[sf]* + fe_v.shape_grad_component(sf, q, di)[d]; + } // for d } @@ -601,90 +596,90 @@ void ConsLaw::assemble_cell_term( // 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 > > flux(n_q_points, - std::vector >(get_n_components(), - std::vector(dim, 0))); + std::vector > > flux(n_q_points, + std::vector >(get_n_components(), + std::vector(dim, 0))); - for (unsigned int q=0; q < n_q_points; ++q) { - Flux(flux[q], fe_v.get_quadrature_points()[q], Wcn[q]); - } + for (unsigned int q=0; q < n_q_points; ++q) { + Flux(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 right_hand_side vector, - // we negate this residual. - for (unsigned int i=0; i right_hand_side vector, + // we negate this residual. + for (unsigned int i=0; iSumIntoGlobalValues(dofs[i], - dofs_per_cell, &values[0], reinterpret_cast(&dofs[0])); + // 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(&dofs[0])); - // Add minus the residual to the right hand side. - right_hand_side(dofs[i]) -= F_i.val(); + // Add minus the residual to the right hand side. + right_hand_side(dofs[i]) -= F_i.val(); - } // for i + } // for i } - // @sect4{%Function: assemble_face_term} - // These are either - // boundary terms or terms across differing - // levels of refinement. In the first case, - // fe_v==fe_v_neighbor and dofs==dofs_neighbor. - // The int boundary < 0 if not at a boundary, - // otherwise it is the boundary indicator. + // @sect4{%Function: assemble_face_term} + // These are either + // boundary terms or terms across differing + // levels of refinement. In the first case, + // fe_v==fe_v_neighbor and dofs==dofs_neighbor. + // The int boundary < 0 if not at a boundary, + // otherwise it is the boundary indicator. template void ConsLaw::assemble_face_term( int face_no, @@ -693,7 +688,7 @@ void ConsLaw::assemble_face_term( std::vector &dofs, std::vector &dofs_neighbor, int boundary - ) +) { fad_double F_i; const unsigned int n_q_points = fe_v.n_quadrature_points; @@ -702,101 +697,101 @@ void ConsLaw::assemble_face_term( Assert(dofs_per_cell == ndofs_per_cell, ExcDimensionMismatch(dofs_per_cell, ndofs_per_cell)); - // As above, the fad degrees of freedom + // As above, the fad degrees of freedom std::vector DOF(dofs_per_cell+ndofs_per_cell); - // The conservative variables for this cell, - // and for + // The conservative variables for this cell, + // and for std::vector > Wplus (n_q_points, - std::vector(get_n_components())); + std::vector(get_n_components())); std::vector > Wminus (n_q_points, - std::vector(get_n_components())); + std::vector(get_n_components())); const std::vector &JxW = fe_v.get_JxW_values (); const std::vector > &normals = fe_v.get_normal_vectors (); - // If we are at a boundary, then dofs_neighbor are - // the same as dofs, so we do not want to duplicate them. - // If there is a neighbor cell, then we want to include - // them. + // 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. + // Set the local DOFS. for (unsigned int in = 0; in < dofs_per_cell; in++) { - DOF[in] = nlsolution(dofs[in]); - DOF[in].diff(in, ndofs); + DOF[in] = nlsolution(dofs[in]); + DOF[in].diff(in, ndofs); } - // If present, set the neighbor dofs. + // If present, set the neighbor dofs. if (boundary < 0) - for (unsigned int in = 0; in < ndofs_per_cell; in++) { + for (unsigned 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. + // Set the values of the local conservative variables. + // Initialize all variables to zero. for (unsigned int q = 0; q < n_q_points; q++) { for (unsigned int di = 0; di < get_n_components(); di++) { - Wplus[q][di] = 0; - Wminus[q][di] = 0; + Wplus[q][di] = 0; + Wminus[q][di] = 0; } for (unsigned int sf = 0; sf < dofs_per_cell; sf++) { - int di = fe_v.get_fe().system_to_component_index(sf).first; - Wplus[q][di] += - (theta*DOF[sf]+(1.0-theta)*solution(dofs[sf]))*fe_v.shape_value_component(sf, q, di); + 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 there is a cell across, then initialize + // the exterior trace as a function of the other + // cell degrees of freedom. if (boundary < 0) { for (unsigned int sf = 0; sf < ndofs_per_cell; sf++) { - int di = fe_v_neighbor.get_fe().system_to_component_index(sf).first; - Wminus[q][di] += - (theta*DOF[sf+dofs_per_cell]+(1.0-theta)*solution(dofs_neighbor[sf]))* - fe_v_neighbor.shape_value_component(sf, q, di); + 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 > bvals(n_q_points, Vector(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 (unsigned int q = 0; q < n_q_points; q++) { + } // 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 > bvals(n_q_points, Vector(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 (unsigned int q = 0; q < n_q_points; q++) { for (unsigned int di = 0; di < get_n_components(); di++) { - // An inflow/dirichlet type of boundary condition + // 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. + // 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; @@ -810,22 +805,22 @@ void ConsLaw::assemble_face_term( 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. + // 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; + 0.5*rho_vel_sqr; } else if (bme->second.first[di] == OUTFLOW_BC) { - // A free/outflow boundary, very simple. + // 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. + // 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 (unsigned int d = 0; d < dim; d++) { vdotn += Wplus[q][d]*normals[q](d); @@ -834,60 +829,60 @@ void ConsLaw::assemble_face_term( Wminus[q][di] = Wplus[q][di] - 2.0*vdotn*normals[q](di); } } - } // for q - } // b>= 0 + } // for q + } // b>= 0 - // Determine the Lax-Friedrich's stability parameter, - // and evaluate the numerical flux function at the quadrature points - std::vector > nflux(n_q_points, std::vector(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(nflux, fe_v.get_quadrature_points(), normals, Wplus, Wminus, - alpha); - - // Now assemble the face term - for (unsigned int i=0; iSumIntoGlobalValues(dofs[i], - dofs_per_cell, &values[0], reinterpret_cast(&dofs[0])); - if (boundary < 0) { - Matrix->SumIntoGlobalValues(dofs[i], - dofs_per_cell, &values[dofs_per_cell], reinterpret_cast(&dofs_neighbor[0])); - } + // Determine the Lax-Friedrich's stability parameter, + // and evaluate the numerical flux function at the quadrature points + std::vector > nflux(n_q_points, std::vector(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(nflux, fe_v.get_quadrature_points(), normals, Wplus, Wminus, + alpha); + + // Now assemble the face term + for (unsigned int i=0; iSumIntoGlobalValues(dofs[i], + dofs_per_cell, &values[0], reinterpret_cast(&dofs[0])); + if (boundary < 0) { + Matrix->SumIntoGlobalValues(dofs[i], + dofs_per_cell, &values[dofs_per_cell], reinterpret_cast(&dofs_neighbor[0])); + } + + // And add into the residual + right_hand_side(dofs[i]) -= F_i.val(); + } + } // @sect4{Assembling the whole system} // Now we put all of the assembly pieces together @@ -900,8 +895,8 @@ void ConsLaw::assemble_system (double &res_norm) FESystem &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. + // We track the dofs on this cell and (if necessary) + // the adjacent cell. std::vector dofs (dofs_per_cell); std::vector dofs_neighbor (dofs_per_cell); @@ -1020,11 +1015,11 @@ void ConsLaw::assemble_system (double &res_norm) // 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()); + face_no, fe_v_face, + fe_v_face, + dofs, + dofs, + face->boundary_indicator()); } else { @@ -1038,7 +1033,7 @@ void ConsLaw::assemble_system (double &res_norm) if (face->has_children()) { - // case I: This cell refined compared to neighbor + // case I: This cell refined compared to neighbor const unsigned int neighbor2= cell->neighbor_of_neighbor(face_no); @@ -1064,13 +1059,13 @@ void ConsLaw::assemble_system (double &res_norm) 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 as if we are working with + // a DG element. assemble_face_term( - face_no, fe_v_subface, - fe_v_face_neighbor, - dofs, - dofs_neighbor); + face_no, fe_v_subface, + fe_v_face_neighbor, + dofs, + dofs_neighbor); } // End of ``if @@ -1078,13 +1073,13 @@ void ConsLaw::assemble_system (double &res_norm) } else { - // We have no children, but - // the neighbor cell may be refine - // compared to use + // 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 + // case II: This is refined compared to neighbor Assert(neighbor->level() < cell->level(), ExcInternalError()); const std::pair faceno_subfaceno= cell->neighbor_of_coarser_neighbor(face_no); @@ -1107,31 +1102,31 @@ void ConsLaw::assemble_system (double &res_norm) neighbor_subface_no); assemble_face_term( - face_no, fe_v_face, - fe_v_subface_neighbor, - dofs, - dofs_neighbor); + face_no, fe_v_face, + fe_v_subface_neighbor, + dofs, + dofs_neighbor); } } - // End of ``face not at boundary'': + // End of ``face not at boundary'': } - // End of loop over all faces: + // End of loop over all faces: } - // End iteration through cells. + // End iteration through cells. } - // Notify Epetra that the matrix is done. - Matrix->FillComplete(); + // Notify Epetra that the matrix is done. + Matrix->FillComplete(); - // Compute the nonlinear residual. - res_norm = right_hand_side.l2_norm(); + // Compute the nonlinear residual. + res_norm = right_hand_side.l2_norm(); } - // Create a conservation law with some defaults. + // Create a conservation law with some defaults. template ConsLaw::ConsLaw () : @@ -1149,14 +1144,14 @@ ConsLaw::ConsLaw () theta(0.5) {} - // At one time this example could work for both DG and - // continuous finite elements. The choice was made here. + // At one time this example could work for both DG and + // continuous finite elements. The choice was made here. template void ConsLaw::build_fe() { fe_ptr = new FESystem(FE_Q(1), N_COMP); } - // Bye bye Conservation law. + // Bye bye Conservation law. template ConsLaw::~ConsLaw () { @@ -1164,13 +1159,13 @@ ConsLaw::~ConsLaw () delete fe_ptr; } - // @sect3{Initialize System} - // Sizes all of the vectors and sets up the - // sparsity patter. This function is called at - // the very beginning of a simulation. The function - // setup_system repeats some of these - // chores and is called after adaptivity in leiu - // of this function. + // @sect3{Initialize System} + // Sizes all of the vectors and sets up the + // sparsity patter. This function is called at + // the very beginning of a simulation. The function + // setup_system repeats some of these + // chores and is called after adaptivity in leiu + // of this function. template void ConsLaw::initialize_system () { @@ -1189,9 +1184,9 @@ void ConsLaw::initialize_system () indicator.reinit(triangulation.n_active_cells()); } - // @sect3{Setup System} - // We call this function to build the sparsity - // and the matrix. + // @sect3{Setup System} + // We call this function to build the sparsity + // and the matrix. template void ConsLaw::setup_system () { @@ -1230,15 +1225,15 @@ void ConsLaw::setup_system () for (unsigned int i=0; i vals(max_nonzero_entries, 0); @@ -1247,20 +1242,20 @@ void ConsLaw::setup_system () unsigned int cur_row = 0; unsigned int cur_col = 0; for (SparsityPattern::iterator s_i = sparsity_pattern.begin(); - s_i != sparsity_pattern.end(); s_i++) { + 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(); + row_indices[cur_col++] = s_i->column(); } - // The last row. + // 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. + // Epetra requires this function after building or + // filling a matrix. It typically does some parallel + // bookeeping; perhaps more. Matrix->FillComplete(); } @@ -1272,93 +1267,93 @@ template void ConsLaw::solve (Vector &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()); + // 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. + // The Direct option selects the Amesos solver. if (solver_params.SOLVER == solver_params_type::DIRECT) { - // Setup for solving with - // Amesos. Other solvers are - // available and may be selected by - // changing th string given to the - // Create function. - Epetra_LinearProblem prob; - prob.SetOperator(Matrix); - Amesos_BaseSolver *solver = Amesos().Create ("Amesos_Klu", 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; + // Setup for solving with + // Amesos. Other solvers are + // available and may be selected by + // changing th string given to the + // Create function. + Epetra_LinearProblem prob; + prob.SetOperator(Matrix); + Amesos_BaseSolver *solver = Amesos().Create ("Amesos_Klu", 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(); + solver->SymbolicFactorization(); - if (solver_params.OUTPUT == solver_params_type::VERBOSE) - std::cout << "Starting Numeric fact\n" << std::flush; + if (solver_params.OUTPUT == solver_params_type::VERBOSE) + std::cout << "Starting Numeric fact\n" << std::flush; - solver->NumericFactorization(); + 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; + // 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. + // For the iterative solvers, we use Aztec. AztecOO Solver; - // Select the appropriate level of verbosity. + // 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. + // 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. + // 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. + // 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. + // 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(); @@ -1378,19 +1373,19 @@ void ConsLaw::postprocess() { | update_q_points | update_JxW_values; UpdateFlags update_flags1 = update_values - | update_gradients - | update_q_points - | update_JxW_values; + | update_gradients + | update_q_points + | update_JxW_values; - QGauss quadrature_formula(4); + QGauss quadrature_formula(4); - const std::vector > &us = fe_ptr->base_element(0).get_unit_support_points(); + const std::vector > &us = fe_ptr->base_element(0).get_unit_support_points(); - Quadrature unit_support(us); + Quadrature unit_support(us); - int n_q_points = quadrature_formula.n_quadrature_points; - int n_uq_points = unit_support.n_quadrature_points; + int n_q_points = quadrature_formula.n_quadrature_points; + int n_uq_points = unit_support.n_quadrature_points; FEValues fe_v ( mapping, *fe_ptr, quadrature_formula, update_flags); @@ -1401,15 +1396,15 @@ void ConsLaw::postprocess() { std::vector > U(n_uq_points, Vector(get_n_components())); std::vector > UU(n_q_points, - Vector(get_n_components())); + Vector(get_n_components())); std::vector > > dU(n_uq_points, - std::vector >(get_n_components())); + std::vector >(get_n_components())); typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); - // Loop the cells + // Loop the cells for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) { cell->get_dof_indices (dofs); fe_v_unit.reinit(cell); @@ -1429,11 +1424,11 @@ void ConsLaw::postprocess() { rho_normVsqr += solution(dofs[vidx])*solution(dofs[vidx]); } rho_normVsqr /= solution(dofs[didx]); - // Pressure + // 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. + // Either output density or gradient squared of density, + // depending on what the user wants. if (!schlieren_plot) { ppsolution(dofs[didx]) = solution(dofs[didx]); } else { @@ -1448,9 +1443,9 @@ void ConsLaw::postprocess() { } - // Loop and assign a value for refinement. We - // simply use the density squared, which selects - // shocks with some success. + // Loop and assign a value for refinement. We + // simply use the density squared, which selects + // shocks with some success. template void ConsLaw::estimate() { @@ -1471,7 +1466,7 @@ void ConsLaw::estimate() { std::vector > U(n_q_points, Vector(get_n_components())); std::vector > > dU(n_q_points, - std::vector >(get_n_components())); + std::vector >(get_n_components())); typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), @@ -1505,11 +1500,11 @@ void ConsLaw::refine_grid () 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. + // 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(); @@ -1518,13 +1513,13 @@ void ConsLaw::refine_grid () cell->set_refine_flag(); } else { if (cell->level() > 0 && - std::fabs(indicator(cell_no)) < 0.75*refinement_params.shock_val) - cell->set_coarsen_flag(); + 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. + // The following code prolongs the solution + // to the new grid and carries out the refinement. std::vector > interp_in; std::vector > interp_out; @@ -1540,24 +1535,24 @@ void ConsLaw::refine_grid () dof_handler.distribute_dofs (*fe_ptr); { - Vector new_solution(1); - Vector new_predictor(1); + Vector new_solution(1); + Vector 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()); + 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 + // 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. + // resize these vectors for the new grid. nlsolution.reinit(dof_handler.n_dofs()); ppsolution.reinit(dof_handler.n_dofs()); nlsolution = solution; @@ -1600,15 +1595,15 @@ void ConsLaw::output_results (const unsigned int cycle) const output.close(); } - // @sect3{Parsing the Input Deck} - // 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. + // @sect3{Parsing the Input Deck} + // 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 void ConsLaw::declare_parameters() { @@ -1619,160 +1614,160 @@ void ConsLaw::declare_parameters() { "intput file"); prm.declare_entry("diffusion power", "2.0", - Patterns::Double(), - "power of mesh size for diffusion"); + Patterns::Double(), + "power of mesh size for diffusion"); prm.declare_entry("gravity", "0.0", - Patterns::Double(), - "gravity forcing"); + 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.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 + // Declare the boundary parameters for (unsigned 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"), - ""); - // declare a slot for each of the conservative - // variables. + Patterns::Selection("true|false"), + ""); + // declare a slot for each of the conservative + // variables. for (unsigned 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"), - ""); + Patterns::Selection( + "inflow|outflow|pressure"), + ""); - // for dirichlet, a function in x,y,z + // 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"); + Patterns::Anything(), + "expression in x,y,z"); } prm.leave_subsection(); } - // Initial condition block. + // Initial condition block. prm.enter_subsection("initial condition"); - for (unsigned int di = 0; di < N_COMP; di++) { - char var[512]; - std::sprintf(var, "w_%d", di); + for (unsigned 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"); - } + // 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. + // The linear solver block. prm.enter_subsection("linear solver"); - prm.declare_entry("output", "quiet", - Patterns::Selection( - "quiet|verbose"), - ""); - prm.declare_entry("method", "gmres", - Patterns::Selection( - "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.declare_entry("output", "quiet", + Patterns::Selection( + "quiet|verbose"), + ""); + prm.declare_entry("method", "gmres", + Patterns::Selection( + "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. + // A refinement controller block. prm.enter_subsection("refinement"); - prm.declare_entry("refinement", "none", - Patterns::Selection( - "none|fixed number|shock"), - ""); - 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.declare_entry("refinement", "none", + Patterns::Selection( + "none|fixed number|shock"), + ""); + 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. + // Output control. prm.enter_subsection("output"); - prm.declare_entry("density", "standard", - Patterns::Selection( - "standard|schlieren"), - ""); - prm.declare_entry("step", "-1", - Patterns::Double(), - "output once per this period"); + prm.declare_entry("density", "standard", + Patterns::Selection( + "standard|schlieren"), + ""); + prm.declare_entry("step", "-1", + Patterns::Double(), + "output once per this period"); prm.leave_subsection(); - // Flux control + // Flux control prm.enter_subsection("flux"); - prm.declare_entry("stab", "alpha", - Patterns::Selection( - "alpha|constant|mesh"), - ""); - prm.declare_entry("stab value", "1", - Patterns::Double(), - "alpha stabilization"); + prm.declare_entry("stab", "alpha", + Patterns::Selection( + "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. + // Code to actually parse an input file. This function + // matches the declarations above. template void ConsLaw::load_parameters(const char *infile){ prm.read_input(infile); - // The global parameters. + // The global parameters. mesh = prm.get("mesh"); diffusion_power = prm.get_double("diffusion power"); gravity = prm.get_double("gravity"); - // The time stepping. + // The time stepping. prm.enter_subsection("time stepping"); dT = prm.get_double("time step"); std::cout << "dT=" << dT << std::endl; @@ -1786,12 +1781,12 @@ void ConsLaw::load_parameters(const char *infile){ std::cout << "TF=" << TF << std::endl; prm.leave_subsection(); - // The boundary info + // The boundary info for (unsigned int b = 0; b < MAX_BD; b++) { std::vector flags(N_COMP, OUTFLOW_BC); - // Define a parser for every boundary, though it may be - // unused. + // Define a parser for every boundary, though it may be + // unused. SideCondition *sd = new SideCondition(N_COMP); char bd[512]; std::sprintf(bd, "boundary_%d", b); @@ -1799,7 +1794,7 @@ void ConsLaw::load_parameters(const char *infile){ const std::string &nopen = prm.get("no penetration"); - // Determine how each component is handled. + // Determine how each component is handled. for (unsigned int di = 0; di < N_COMP; di++) { char var[512]; std::sprintf(var, "w_%d", di); @@ -1819,79 +1814,79 @@ void ConsLaw::load_parameters(const char *infile){ } prm.leave_subsection(); - // Add the boundary condition to the law. + // Add the boundary condition to the law. sd->Init(); add_boundary(b, flags, sd); - } + } - // Initial conditions. - prm.enter_subsection("initial condition"); - for (unsigned int di = 0; di < N_COMP; di++) { - char var[512]; + // Initial conditions. + prm.enter_subsection("initial condition"); + for (unsigned 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; - } + 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(); - 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"); + // 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. + // 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 + 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"); + 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. + // 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"); + const std::string &dens = prm.get("density"); + schlieren_plot = dens == "schlieren" ? true : false; + output_step = prm.get_double("step"); prm.leave_subsection(); - // Flux control. + // 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"); + 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(); @@ -1902,25 +1897,25 @@ void ConsLaw::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. + // 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 void ConsLaw::compute_predictor() { predictor = nlsolution; predictor.sadd(3/2.0, -1/2.0, solution); } - // @sect3{Run the simulation} - // Contains the initialization - // the time loop, and the inner Newton iteration. + // @sect3{Run the simulation} + // Contains the initialization + // the time loop, and the inner Newton iteration. template void ConsLaw::run () { - // Open and load the mesh. + // Open and load the mesh. GridIn grid_in; grid_in.attach_triangulation(triangulation); std::cout << "Opening mesh <" << mesh << ">" << std::endl; @@ -1935,30 +1930,30 @@ void ConsLaw::run () unsigned int nstep = 0; - // Initialize fields and matrices. + // Initialize fields and matrices. initialize_system (); setup_system(); initialize(); predictor = solution; - // Initial refinement. We apply the ic, - // estimate, refine, and repeat until - // happy. + // Initial refinement. We apply the ic, + // estimate, refine, and repeat until + // happy. if (refinement_params.refine != refinement_params_type::NONE) - for (unsigned int i = 0; i < refinement_params.shock_levels; i++) { - estimate(); - refine_grid(); - setup_system(); - initialize(); - predictor = solution; - } + for (unsigned 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. + // Determine when we will output next. double next_output = T + output_step; - // @sect4{Main time stepping loop} + // @sect4{Main time stepping loop} predictor = solution; while(T < TF) { @@ -1978,8 +1973,8 @@ void ConsLaw::run () double res_norm; int lin_iter; - // Print some relevant information during the - // Newton iteration. + // Print some relevant information during the + // Newton iteration. std::cout << "NonLin Res: Lin Iter Lin Res" << std::endl; std::cout << "______________________________________" << std::endl; @@ -1987,35 +1982,35 @@ void ConsLaw::run () unsigned int nonlin_iter = 0; double lin_res; - // @sect5{Newton iteration} + // @sect5{Newton iteration} 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. + // 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). + // 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); + // 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. + // Print the residuals. std::printf("%-16.3e %04d %-5.2e\n", - res_norm, lin_iter, lin_res); + res_norm, lin_iter, lin_res); ++nonlin_iter; @@ -2023,7 +2018,7 @@ void ConsLaw::run () ExcMessage ("No convergence in nonlinear solver")); } - // Various post convergence tasks. + // Various post convergence tasks. compute_predictor(); solution = nlsolution; @@ -2035,7 +2030,7 @@ void ConsLaw::run () T += dT; - // Output if it is time. + // Output if it is time. if (output_step < 0) { output_results (++nstep); } else if (T >= next_output) { @@ -2043,7 +2038,7 @@ void ConsLaw::run () next_output += output_step; } - // Refine, if refinement is selected. + // Refine, if refinement is selected. if (refinement_params.refine != refinement_params_type::NONE) { refine_grid(); setup_system(); -- 2.39.5