From c41d65fcecb67c8d50aeb8c2a9bd24a44d097206 Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 19 May 2008 21:18:46 +0000 Subject: [PATCH] Some more. git-svn-id: https://svn.dealii.org/trunk@16126 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-33/step-33.cc | 994 ++++++++++++++-------------- 1 file changed, 504 insertions(+), 490 deletions(-) diff --git a/deal.II/examples/step-33/step-33.cc b/deal.II/examples/step-33/step-33.cc index 66db71abf1..6336501c78 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -23,8 +23,7 @@ #include #include -#include -#include +#include #include #include @@ -118,24 +117,24 @@ template struct EulerEquations { // First a few variables that - // describe the various components of our - // solution vector in a generic way. This - // includes the number of components in the - // system (Euler's equations have one entry - // for momenta in each spatial direction, - // plus the energy and density components, - // for a total of dim+2 - // components), as well as functions that - // describe the index within the solution - // vector of the first momentum component, - // the density component, and the energy - // density component. Note that all these - // %numbers depend on the space dimension; - // defining them in a generic way (rather - // than by implicit convention) makes our - // code more flexible and makes it easier - // to later extend it, for example by - // adding more components to the equations. + // describe the various components of our + // solution vector in a generic way. This + // includes the number of components in the + // system (Euler's equations have one entry + // for momenta in each spatial direction, + // plus the energy and density components, + // for a total of dim+2 + // components), as well as functions that + // describe the index within the solution + // vector of the first momentum component, + // the density component, and the energy + // density component. Note that all these + // %numbers depend on the space dimension; + // defining them in a generic way (rather + // than by implicit convention) makes our + // code more flexible and makes it easier + // to later extend it, for example by + // adding more components to the equations. static const unsigned int n_components = dim + 2; static const unsigned int first_momentum_component = 0; static const unsigned int density_component = dim; @@ -185,22 +184,22 @@ struct EulerEquations } - // Next, we define the gas - // constant. We will set it to 1.4 - // in its definition immediately - // following the declaration of - // this class (unlike integer - // variables, like the ones above, - // static const floating point - // member variables cannot be - // initialized within the class - // declaration in C++). This value - // of 1.4 is representative of a - // gas that consists of molecules - // composed of two atoms, such as - // air which consists up to small - // traces almost entirely of $N_2$ - // and $O_2$. + // Next, we define the gas + // constant. We will set it to 1.4 + // in its definition immediately + // following the declaration of + // this class (unlike integer + // variables, like the ones above, + // static const floating point + // member variables cannot be + // initialized within the class + // declaration in C++). This value + // of 1.4 is representative of a + // gas that consists of molecules + // composed of two atoms, such as + // air which consists up to small + // traces almost entirely of $N_2$ + // and $O_2$. static const double gas_gamma; @@ -368,7 +367,7 @@ struct EulerEquations normal_flux[di] += 0.5*alpha*(Wplus[di] - Wminus[di]); } - } + } // Finally, we declare a class that @@ -522,9 +521,9 @@ compute_derived_quantities_vector (const std::vector > if (do_schlieren_plot == true) Assert (duh.size() == n_quadrature_points, ExcInternalError()) - else - Assert (duh.size() == 0, - ExcInternalError()); + else + Assert (duh.size() == 0, + ExcInternalError()); Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError()); @@ -535,9 +534,9 @@ compute_derived_quantities_vector (const std::vector > if (do_schlieren_plot == true) Assert (computed_quantities[0].size() == dim+2, ExcInternalError()) - else - Assert (computed_quantities[0].size() == dim+1, - ExcInternalError()); + else + Assert (computed_quantities[0].size() == dim+1, + ExcInternalError()); // Then loop over all quadrature points and // do our work there. The code should be @@ -1349,18 +1348,14 @@ class ConservationLaw { public: ConservationLaw (const char *input_filename); - ~ConservationLaw (); - void run (); private: void setup_system (); void assemble_system (); - void assemble_cell_term (const FEValues &fe_v, const std::vector &dofs); - void assemble_face_term(const unsigned int face_no, const FEFaceValuesBase &fe_v, const FEFaceValuesBase &fe_v_neighbor, @@ -1410,7 +1405,6 @@ class ConservationLaw const FESystem fe; DoFHandler dof_handler; - SparsityPattern sparsity_pattern; const QGauss quadrature; const QGauss face_quadrature; @@ -1491,7 +1485,7 @@ class ConservationLaw // std::auto_ptr // instead of a plain pointer for // this. - Epetra_SerialComm communicator; + Epetra_SerialComm communicator; std::auto_ptr Map; std::auto_ptr Matrix; @@ -1499,7 +1493,13 @@ class ConservationLaw }; - // Create a conservation law with some defaults. + // @sect4{ConservationLaw::ConservationLaw} + // + // There is nothing much to say about + // the constructor. Essentially, it + // reads the input file and fills the + // parameter object with the parsed + // values: template ConservationLaw::ConservationLaw (const char *input_filename) : @@ -1517,16 +1517,328 @@ ConservationLaw::ConservationLaw (const char *input_filename) } - // Bye bye Conservation law. + + // @sect4{ConservationLaw::setup_system} + // + // The following function is called + // each time the mesh is + // changed. Essentially what it does + // is to resize the Trilinos + // matrix. In addition to just + // resizing it, it also builds a + // sparsity pattern, initializes the + // row lengths of the matrix with the + // ones from this sparsity pattern, + // and finally puts zero entries into + // the places where nonzero entries + // will later be found. This will + // make subsequent operations on the + // matrix faster, because no new + // memory will need to be allocated: template -ConservationLaw::~ConservationLaw () +void ConservationLaw::setup_system () { - dof_handler.clear (); + Map.reset (new Epetra_Map(dof_handler.n_dofs(), 0, communicator)); + + + // Now create a sparsity pattern, + // condense it, and count the + // number of nonzero entries per + // row: + CompressedSparsityPattern sparsity_pattern (dof_handler.n_dofs(), + dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + + ConstraintMatrix hanging_node_constraints; + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + hanging_node_constraints.close (); + + hanging_node_constraints.condense (sparsity_pattern); + sparsity_pattern.compress(); + + std::vector row_lengths (dof_handler.n_dofs()); + for (unsigned int i=0; iFillComplete() + // function that indicates that no + // other matrix entries will be + // needed: + Matrix.reset (new Epetra_CrsMatrix(Copy, *Map, &row_lengths[0], true)); + + const unsigned int max_nonzero_entries + = *std::max_element (row_lengths.begin(), row_lengths.end()); + + std::vector values(max_nonzero_entries, 0); + std::vector row_indices(max_nonzero_entries); + + for (unsigned int row=0; rowInsertGlobalValues(row, row_lengths[row], + &values[0], &row_indices[0]); + } + + Matrix->FillComplete(); } - // @sect3{Assembly} - // @sect4{%Function: assemble_cell_term} + // @sect4{ConservationLaw::assemble_system} + // + // This and the following two + // functions are the meat of this + // program: They assemble the linear + // system that results from applying + // Newton's method to the nonlinear + // system of conservation + // equations. + // + // This first function puts all of + // the assembly pieces together in a + // routine that dispatches the + // correct piece for each cell/face. + // The actual implementation of the + // assembly on these objects is done + // in the following functions. +template +void ConservationLaw::assemble_system () +{ + 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 dofs (dofs_per_cell); + std::vector dofs_neighbor (dofs_per_cell); + + // First we create the + // ``UpdateFlags'' for the + // ``FEValues'' and the + // ``FEFaceValues'' objects. + const 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. + 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. + 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 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 fe_v_face (mapping, fe, face_quadrature, + face_update_flags); + FESubfaceValues fe_v_subface (mapping, fe, face_quadrature, + face_update_flags); + FEFaceValues fe_v_face_neighbor (mapping, fe, face_quadrature, + neighbor_face_update_flags); + FESubfaceValues fe_v_subface_neighbor (mapping, fe, face_quadrature, + neighbor_face_update_flags); + + // Furthermore we need some cell + // iterators. + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + // Now we start the loop over all + // active cells. + for (; cell!=endc; ++cell) + { + + // 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); + + assemble_cell_term(fe_v, dofs); + + // 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::faces_per_cell; ++face_no) + { + if (cell->at_boundary(face_no)) + { + // 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, + true, + cell->face(face_no)->boundary_indicator(), + cell->face(face_no)->diameter()); + } + else + { + // Now we are not on + // the boundary of the + // domain, therefore + // there must exist a + // neighboring cell. + typename DoFHandler::cell_iterator neighbor= + cell->neighbor(face_no);; + + if (cell->face(face_no)->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::subfaces_per_face; + ++subface_no) + { + typename DoFHandler::active_cell_iterator + neighbor_child + = cell->neighbor_child_on_subface (face_no, subface_no); + + 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, + false, + numbers::invalid_unsigned_int, + neighbor_child->diameter()); + } + // 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 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, + false, + numbers::invalid_unsigned_int, + cell->face(face_no)->diameter()); + } + + } + // 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(); +} + + + // @sect4{ConservationLaw::assemble_cell_term} // // Assembles the cell term, adding minus the residual // to the right hand side, and adding in the Jacobian @@ -1543,7 +1855,7 @@ void ConservationLaw::assemble_cell_term (const FEValues & // Values of the conservative variables at the quadrature points. std::vector > > W (n_q_points, - std::vector >(EulerEquations::n_components)); + std::vector >(EulerEquations::n_components)); // Values at the last time step of the conservative variables. // Note that these do not use fad variables, since they do @@ -1554,7 +1866,7 @@ void ConservationLaw::assemble_cell_term (const FEValues & // 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 >(EulerEquations::n_components)); + std::vector >(EulerEquations::n_components)); // Gradients of the current variables. It is a // bit of a shame that we have to compute these; we almost don't. @@ -1562,8 +1874,8 @@ void ConservationLaw::assemble_cell_term (const FEValues & // 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 > >(EulerEquations::n_components, - std::vector >(dim))); + std::vector > >(EulerEquations::n_components, + std::vector >(dim))); // Here is the magical point where we declare a subset @@ -1696,7 +2008,10 @@ void ConservationLaw::assemble_cell_term (const FEValues & delete[] flux; } - // @sect4{%Function: assemble_face_term} + + + // @sect4{ConservationLaw::assemble_face_term} + // // These are either // boundary terms or terms across differing // levels of refinement. In the first case, @@ -1727,9 +2042,9 @@ ConservationLaw::assemble_face_term(const unsigned int face_no, // The conservative variables for this cell, // and for std::vector > > Wplus (n_q_points, - std::vector >(EulerEquations::n_components)); + std::vector >(EulerEquations::n_components)); std::vector > > Wminus (n_q_points, - std::vector >(EulerEquations::n_components)); + std::vector >(EulerEquations::n_components)); const std::vector > &normals = fe_v.get_normal_vectors (); @@ -1795,80 +2110,80 @@ ConservationLaw::assemble_face_term(const unsigned int face_no, ExcIndexRange (boundary_id, 0, Parameters::AllParameters::max_n_boundaries)); - // Evaluate the function object. This is - // a bit tricky; a given boundary might - // have both prescribed and implicit - // values. If a particular component is - // not prescribed, the values evaluate to - // zero and are ignored, below. - std::vector > bvals(n_q_points, Vector(EulerEquations::n_components)); - parameters.boundary_conditions[boundary_id].values.vector_value_list(fe_v.get_quadrature_points(), bvals); - - // We loop the quadrature points, and we treat each - // component individualy. - for (unsigned int q = 0; q < n_q_points; q++) { - for (unsigned int di = 0; di < EulerEquations::n_components; di++) { - - // An inflow/dirichlet type of boundary condition - if (parameters.boundary_conditions[boundary_id].kind[di] == Parameters::AllParameters::inflow_boundary) { - Wminus[q][di] = bvals[q](di); - } else if (parameters.boundary_conditions[boundary_id].kind[di] == Parameters::AllParameters::pressure_boundary) { - // A prescribed pressure boundary - // condition. This boundary - // condition is complicated by the - // fact that even though the - // pressure is prescribed, we - // really are setting the energy - // index here, which will depend on - // velocity and pressure. So even - // though this seems like a - // dirichlet type boundary - // condition, we get sensitivities - // of energy to velocity and - // density (unless these are also - // prescribed. - Sacado::Fad::DFad rho_vel_sqr = 0; - Sacado::Fad::DFad dens; + // Evaluate the function object. This is + // a bit tricky; a given boundary might + // have both prescribed and implicit + // values. If a particular component is + // not prescribed, the values evaluate to + // zero and are ignored, below. + std::vector > bvals(n_q_points, Vector(EulerEquations::n_components)); + parameters.boundary_conditions[boundary_id].values.vector_value_list(fe_v.get_quadrature_points(), bvals); + + // We loop the quadrature points, and we treat each + // component individualy. + for (unsigned int q = 0; q < n_q_points; q++) { + for (unsigned int di = 0; di < EulerEquations::n_components; di++) { + + // An inflow/dirichlet type of boundary condition + if (parameters.boundary_conditions[boundary_id].kind[di] == Parameters::AllParameters::inflow_boundary) { + Wminus[q][di] = bvals[q](di); + } else if (parameters.boundary_conditions[boundary_id].kind[di] == Parameters::AllParameters::pressure_boundary) { + // A prescribed pressure boundary + // condition. This boundary + // condition is complicated by the + // fact that even though the + // pressure is prescribed, we + // really are setting the energy + // index here, which will depend on + // velocity and pressure. So even + // though this seems like a + // dirichlet type boundary + // condition, we get sensitivities + // of energy to velocity and + // density (unless these are also + // prescribed. + Sacado::Fad::DFad rho_vel_sqr = 0; + Sacado::Fad::DFad dens; - dens = parameters.boundary_conditions[boundary_id].kind[EulerEquations::density_component] == Parameters::AllParameters::inflow_boundary ? bvals[q](EulerEquations::density_component) : - Wplus[q][EulerEquations::density_component]; - - for (unsigned int d=0; d < dim; d++) { - if (parameters.boundary_conditions[boundary_id].kind[d] == Parameters::AllParameters::inflow_boundary) - rho_vel_sqr += bvals[q](d)*bvals[q](d); - else - rho_vel_sqr += Wplus[q][d]*Wplus[q][d]; - } - rho_vel_sqr /= dens; - // Finally set the energy value as determined by the - // prescribed pressure and the other variables. - Wminus[q][di] = bvals[q](di)/(EulerEquations::gas_gamma-1.0) + - 0.5*rho_vel_sqr; - - } else if (parameters.boundary_conditions[boundary_id].kind[di] == Parameters::AllParameters::outflow_boundary) { - // A free/outflow boundary, very simple. - Wminus[q][di] = Wplus[q][di]; - - } else { - // We must be at a no-penetration - // boundary. We prescribe the - // velocity (we are dealing with a - // particular component here so - // that the average of the - // velocities is orthogonal to the - // surface normal. This creates - // sensitivies of across the - // velocity components. - Sacado::Fad::DFad vdotn = 0; - for (unsigned int d = 0; d < dim; d++) { - vdotn += Wplus[q][d]*normals[q](d); - } - - Wminus[q][di] = Wplus[q][di] - 2.0*vdotn*normals[q](di); - } - } - } // for q - } // b>= 0 + dens = parameters.boundary_conditions[boundary_id].kind[EulerEquations::density_component] == Parameters::AllParameters::inflow_boundary ? bvals[q](EulerEquations::density_component) : + Wplus[q][EulerEquations::density_component]; + + for (unsigned int d=0; d < dim; d++) { + if (parameters.boundary_conditions[boundary_id].kind[d] == Parameters::AllParameters::inflow_boundary) + rho_vel_sqr += bvals[q](d)*bvals[q](d); + else + rho_vel_sqr += Wplus[q][d]*Wplus[q][d]; + } + rho_vel_sqr /= dens; + // Finally set the energy value as determined by the + // prescribed pressure and the other variables. + Wminus[q][di] = bvals[q](di)/(EulerEquations::gas_gamma-1.0) + + 0.5*rho_vel_sqr; + + } else if (parameters.boundary_conditions[boundary_id].kind[di] == Parameters::AllParameters::outflow_boundary) { + // A free/outflow boundary, very simple. + Wminus[q][di] = Wplus[q][di]; + + } else { + // We must be at a no-penetration + // boundary. We prescribe the + // velocity (we are dealing with a + // particular component here so + // that the average of the + // velocities is orthogonal to the + // surface normal. This creates + // sensitivies of across the + // velocity components. + Sacado::Fad::DFad vdotn = 0; + for (unsigned int d = 0; d < dim; d++) { + vdotn += Wplus[q][d]*normals[q](d); + } + + Wminus[q][di] = Wplus[q][di] - 2.0*vdotn*normals[q](di); + } + } + } // for q + } // b>= 0 // Determine the Lax-Friedrich's stability parameter, // and evaluate the numerical flux function at the quadrature points @@ -1936,311 +2251,8 @@ ConservationLaw::assemble_face_term(const unsigned int face_no, delete[] normal_fluxes; } - // @sect4{Assembling the whole system} - // 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 -void ConservationLaw::assemble_system () -{ - 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 dofs (dofs_per_cell); - std::vector 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 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 fe_v_face ( - mapping, fe, face_quadrature, face_update_flags); - FESubfaceValues fe_v_subface ( - mapping, fe, face_quadrature, face_update_flags); - FEFaceValues fe_v_face_neighbor ( - mapping, fe, face_quadrature, neighbor_face_update_flags); - FESubfaceValues fe_v_subface_neighbor ( - mapping, fe, face_quadrature, neighbor_face_update_flags); - - // Furthermore we need some cell - // iterators. - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - // Now we start the loop over all - // active cells. - 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); - - assemble_cell_term(fe_v, dofs); - - // 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::faces_per_cell; ++face_no) - { - // First we set the face - // iterator - typename DoFHandler::face_iterator face=cell->face(face_no); - - 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, - true, - face->boundary_indicator(), - face->diameter()); - } - else - { - // Now we are not on - // the boundary of the - // domain, therefore - // there must exist a - // neighboring cell. - typename DoFHandler::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::subfaces_per_face; - ++subface_no) - { - typename DoFHandler::active_cell_iterator - neighbor_child - = cell->neighbor_child_on_subface (face_no, subface_no); - - 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, - false, - numbers::invalid_unsigned_int, - neighbor_child->diameter()); - } - // 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 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, - false, - numbers::invalid_unsigned_int, - face->diameter()); - } - - } - // 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(); -} - - // @sect3{Setup System} - // We call this function to build the sparsity - // and the matrix. -template -void ConservationLaw::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::faces_per_cell - *GeometryInfo::subfaces_per_face+1)*fe.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. - Map.reset (new Epetra_Map(dof_handler.n_dofs(), 0, communicator)); - - // Epetra can build a more efficient matrix if - // one knows ahead of time the maximum number of - // columns in any row entry - std::vector row_lengths (dof_handler.n_dofs()); - for (unsigned int i=0; i vals(max_nonzero_entries, 0); - std::vector row_indices(max_nonzero_entries); - - 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++) { - 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(); -} // @sect3{Solving the linear system} // Actually solve the linear system, using either @@ -2304,44 +2316,44 @@ ConservationLaw::solve (Vector &newton_update) else if (parameters.solver == Parameters::Solver::gmres) { - // For the iterative solvers, we use Aztec. - AztecOO Solver; - - // Select the appropriate level of verbosity. - if (parameters.output == Parameters::Solver::quiet) - Solver.SetAztecOption(AZ_output, AZ_none); - - if (parameters.output == Parameters::Solver::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, parameters.ilut_drop); - Solver.SetAztecParam(AZ_ilut_fill, parameters.ilut_fill); - Solver.SetAztecParam(AZ_athresh, parameters.ilut_atol); - Solver.SetAztecParam(AZ_rthresh, parameters.ilut_rtol); - Solver.SetUserMatrix(Matrix.get()); - - // Run the solver iteration. Collect the number - // of iterations and the residual. - Solver.Iterate(parameters.max_iterations, parameters.linear_residual); - - return std::make_pair (Solver.NumIters(), - Solver.TrueResidual()); - } + // For the iterative solvers, we use Aztec. + AztecOO Solver; + + // Select the appropriate level of verbosity. + if (parameters.output == Parameters::Solver::quiet) + Solver.SetAztecOption(AZ_output, AZ_none); + + if (parameters.output == Parameters::Solver::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, parameters.ilut_drop); + Solver.SetAztecParam(AZ_ilut_fill, parameters.ilut_fill); + Solver.SetAztecParam(AZ_athresh, parameters.ilut_atol); + Solver.SetAztecParam(AZ_rthresh, parameters.ilut_rtol); + Solver.SetUserMatrix(Matrix.get()); + + // Run the solver iteration. Collect the number + // of iterations and the residual. + Solver.Iterate(parameters.max_iterations, parameters.linear_residual); + + return std::make_pair (Solver.NumIters(), + Solver.TrueResidual()); + } Assert (false, ExcNotImplemented()); return std::make_pair (0,0); @@ -2377,22 +2389,24 @@ compute_refinement_indicators (Vector &refinement_indicators) const typename DoFHandler::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); + 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); + fe_v.get_function_values(predictor, U); + fe_v.get_function_grads(predictor, dU); - refinement_indicators(cell_no) = 0; - for (unsigned int q = 0; q < n_q_points; q++) { - double ng = 0; - for (unsigned int d = 0; d < dim; d++) ng += dU[q][EulerEquations::density_component][d]*dU[q][EulerEquations::density_component][d]; + refinement_indicators(cell_no) = 0; + for (unsigned int q = 0; q < n_q_points; q++) { + double ng = 0; + for (unsigned int d = 0; d < dim; d++) + ng += dU[q][EulerEquations::density_component][d] * + dU[q][EulerEquations::density_component][d]; - refinement_indicators(cell_no) += std::log(1+std::sqrt(ng)); + refinement_indicators(cell_no) += std::log(1+std::sqrt(ng)); - } - refinement_indicators(cell_no) /= n_q_points; - + } + refinement_indicators(cell_no) /= n_q_points; } } @@ -2610,8 +2624,8 @@ void ConservationLaw::run () } else { - // Solve the linear system and update with the - // delta. + // Solve the linear system and update with the + // delta. newton_update = 0; std::pair convergence -- 2.39.5