From a4fdef1d1e8c32d2b4b0106ce3ca5e4fb3a875d5 Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 19 May 2008 19:51:11 +0000 Subject: [PATCH] Remove a few class variables and instead pass them as arguments. git-svn-id: https://svn.dealii.org/trunk@16124 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-33/doc/results.dox | 1 + deal.II/examples/step-33/step-33.cc | 515 ++++++++++++----------- 2 files changed, 275 insertions(+), 241 deletions(-) diff --git a/deal.II/examples/step-33/doc/results.dox b/deal.II/examples/step-33/doc/results.dox index c426e7b7aa..b214cdedb7 100644 --- a/deal.II/examples/step-33/doc/results.dox +++ b/deal.II/examples/step-33/doc/results.dox @@ -59,6 +59,7 @@ end subsection time stepping set final time = 10.0 # simulation end time set time step = 0.02 # simulation time step + set theta scheme value = 0.5 end subsection linear solver diff --git a/deal.II/examples/step-33/step-33.cc b/deal.II/examples/step-33/step-33.cc index 7aaeb15511..1ed175886b 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -1097,6 +1097,7 @@ namespace Parameters double gravity; double time_step, final_time; + double theta; bool is_stationary; std::string mesh_filename; @@ -1143,11 +1144,16 @@ namespace Parameters prm.enter_subsection("time stepping"); { prm.declare_entry("time step", "0.1", - Patterns::Double(), + Patterns::Double(0), "simulation time step"); prm.declare_entry("final time", "10.0", - Patterns::Double(), + Patterns::Double(0), "simulation end time"); + prm.declare_entry("theta scheme value", "0.5", + Patterns::Double(0,1), + "value for theta that interpolated between explicit " + "Euler (theta=0), Crank-Nicolson (theta=0.5), and " + "implicit Euler (theta=1)."); } prm.leave_subsection(); @@ -1216,6 +1222,7 @@ namespace Parameters is_stationary = false; final_time = prm.get_double("final time"); + theta = prm.get_double("theta scheme value"); } prm.leave_subsection(); @@ -1283,101 +1290,143 @@ namespace Parameters // @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. Member - // functions will be described as - // their definitions appear. + // Here finally comes the class that + // actually does something with all + // the Euler equation and parameter + // specifics we've defined above. The + // public interface is pretty much + // the same as always (the + // constructor now takes the name of + // a file from which to read + // parameters, which is passed on the + // command line). The private + // function interface is also pretty + // similar to the usual arrangement, + // with the + // assemble_system + // function split into three parts: + // one that contains the main loop + // over all cells and that then calls + // the other two for integrals over + // cells and faces, respectively. template -class ConsLaw +class ConservationLaw { public: - ConsLaw (const char *input_filename); - ~ConsLaw (); + ConservationLaw (const char *input_filename); + ~ConservationLaw (); void run (); private: void setup_system (); - void initialize_system (); - void assemble_system (double &res_norm); - void solve (Vector &solution, int &, double &); + + void assemble_system (); + + void assemble_cell_term (const FEValues &fe_v, + const std::vector &dofs); + + void assemble_face_term(const unsigned int face_no, + const FEFaceValuesBase &fe_v, + const FEFaceValuesBase &fe_v_neighbor, + const std::vector &dofs, + const std::vector &dofs_neighbor, + const bool external_face, + const unsigned int boundary_id, + const double face_diameter); + + std::pair solve (Vector &solution); + void refine_grid (); - void output_results (const unsigned int cycle) const; - void initialize(); - void estimate(); - void compute_predictor(); + void output_results () const; + + void compute_refinement_indicator (); + + + // The first few member variables + // are also rather standard. Note + // that we define a mapping + // object to be used throughout + // the program when assembling + // terms (we will hand it to + // every FEValues and + // FEFaceValues object); the + // mapping we use is just the + // standard $Q_1$ mapping -- + // nothing fancy, in other words + // -- but declaring one here and + // using it throughout the + // program will make it simpler + // later on to change it if that + // should become necessary. This + // is, in fact, rather pertinent: + // it is known that for + // transsonic simulations with + // the Euler equations, + // computations do not converge + // even as $h\rightarrow 0$ if + // the boundary approximation is + // not of sufficiently high + // order. Triangulation triangulation; - const MappingQ1 mapping; + const MappingQ1 mapping; - - FESystem fe; - + const FESystem fe; DoFHandler dof_handler; SparsityPattern sparsity_pattern; - const QGauss quadrature; - const QGauss face_quadrature; + const QGauss quadrature; + const QGauss face_quadrature; - // The actual solution to the Euler equation - Vector solution; - // The current value of the solution during the Newton iteration - Vector nlsolution; - // An estimate of the next time value; used for adaptivity and as a - // guess for the next Newton iteration. + // Next come a number of data + // vectors that correspond to the + // solution of the previous time + // step + // (old_solution), + // the best guess of the current + // solution + // (current_solution; + // we say guess because + // the Newton iteration to + // compute it may not have + // converged yet, whereas + // old_solution + // refers to the fully converged + // final result of the previous + // time step), and a predictor + // for the solution at the next + // time step, computed by + // extrapolating the current and + // previous solution one time + // step into the future: + Vector old_solution; + Vector current_solution; Vector predictor; - // The solution to the linear problem during the Newton iteration - Vector dsolution; + Vector right_hand_side; Epetra_SerialComm communicator; - - public: - - void assemble_cell_term (const FEValues &fe_v, - const std::vector &dofs); - - void assemble_face_term( - int face_no, - const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - std::vector &dofs, - std::vector &dofs_neighbor, - const bool external_face, - const unsigned int boundary_id); - - private: - double T; - double face_diameter; - double cell_diameter; - - Parameters::AllParameters parameters; Epetra_Map *Map; Epetra_CrsMatrix *Matrix; Vector indicator; - // Crank-Nicolson value - const double theta; - + Parameters::AllParameters parameters; }; // Create a conservation law with some defaults. template -ConsLaw::ConsLaw (const char *input_filename) +ConservationLaw::ConservationLaw (const char *input_filename) : mapping (), fe (FE_Q(1), EulerEquations::n_components), dof_handler (triangulation), quadrature (2), face_quadrature (2), - T(0), Map(NULL), - Matrix(NULL), - theta(0.5) + Matrix(NULL) { ParameterHandler prm; Parameters::AllParameters::declare_parameters (prm); @@ -1389,21 +1438,12 @@ ConsLaw::ConsLaw (const char *input_filename) // Bye bye Conservation law. template -ConsLaw::~ConsLaw () +ConservationLaw::~ConservationLaw () { dof_handler.clear (); } - // Apply the initialial condition. Simultaneously - // initialize the non-linear solution. -template -void ConsLaw::initialize() { - VectorTools::interpolate(dof_handler, - parameters.initial_conditions, solution); - nlsolution = solution; -} - // @sect3{Assembly} // @sect4{%Function: assemble_cell_term} // @@ -1411,8 +1451,8 @@ void ConsLaw::initialize() { // to the right hand side, and adding in the Jacobian // contributions. template -void ConsLaw::assemble_cell_term (const FEValues &fe_v, - const std::vector &dofs) +void ConservationLaw::assemble_cell_term (const FEValues &fe_v, + const std::vector &dofs) { unsigned int dofs_per_cell = fe_v.dofs_per_cell; unsigned int n_q_points = fe_v.n_quadrature_points; @@ -1451,14 +1491,14 @@ void ConsLaw::assemble_cell_term (const FEValues &fe_v, // 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] = current_solution(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 + // but since we don't want to make the entire old_solution vector // fad types, only the local cell variables, we explicitly // code this loop; for (unsigned int q = 0; q < n_q_points; q++) { @@ -1475,9 +1515,9 @@ void ConsLaw::assemble_cell_term (const FEValues &fe_v, 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); + old_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); + (parameters.theta*DOF[sf]+(1-parameters.theta)*old_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]* @@ -1541,7 +1581,7 @@ void ConsLaw::assemble_cell_term (const FEValues &fe_v, // Stabilization (cell wise diffusion) for (unsigned int d = 0; d < dim; d++) - F_i += 1.0*std::pow(cell_diameter, parameters.diffusion_power) * + F_i += 1.0*std::pow(fe_v.get_cell()->diameter(), parameters.diffusion_power) * fe_v.shape_grad_component(i, point, component_i)[d] * Wgrads[point][component_i][d] * fe_v.JxW(point); @@ -1583,15 +1623,15 @@ void ConsLaw::assemble_cell_term (const FEValues &fe_v, // The int boundary < 0 if not at a boundary, // otherwise it is the boundary indicator. template -void ConsLaw::assemble_face_term( - int face_no, - const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - std::vector &dofs, - std::vector &dofs_neighbor, - const bool external_face, - const unsigned int boundary_id -) +void +ConservationLaw::assemble_face_term(const unsigned int face_no, + const FEFaceValuesBase &fe_v, + const FEFaceValuesBase &fe_v_neighbor, + const std::vector &dofs, + const std::vector &dofs_neighbor, + const bool external_face, + const unsigned int boundary_id, + const double face_diameter) { Sacado::Fad::DFad F_i; const unsigned int n_q_points = fe_v.n_quadrature_points; @@ -1622,13 +1662,13 @@ void ConsLaw::assemble_face_term( int ndofs = (external_face == false ? dofs_per_cell + ndofs_per_cell : dofs_per_cell); // Set the local DOFS. for (unsigned int in = 0; in < dofs_per_cell; in++) { - DOF[in] = nlsolution(dofs[in]); + DOF[in] = current_solution(dofs[in]); DOF[in].diff(in, ndofs); } // If present, set the neighbor dofs. if (external_face == false) for (unsigned int in = 0; in < ndofs_per_cell; in++) { - DOF[in+dofs_per_cell] = nlsolution(dofs_neighbor[in]); + DOF[in+dofs_per_cell] = current_solution(dofs_neighbor[in]); DOF[in+dofs_per_cell].diff(in+dofs_per_cell, ndofs); } @@ -1642,7 +1682,7 @@ void ConsLaw::assemble_face_term( 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); + (parameters.theta*DOF[sf]+(1.0-parameters.theta)*old_solution(dofs[sf]))*fe_v.shape_value_component(sf, q, di); } @@ -1653,7 +1693,7 @@ void ConsLaw::assemble_face_term( 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]))* + (parameters.theta*DOF[sf+dofs_per_cell]+(1.0-parameters.theta)*old_solution(dofs_neighbor[sf]))* fe_v_neighbor.shape_value_component(sf, q, di); } } @@ -1798,13 +1838,15 @@ void ConsLaw::assemble_face_term( // is/isn't a neighboring cell, we add more/less // entries. Matrix->SumIntoGlobalValues(dofs[i], - dofs_per_cell, &values[0], reinterpret_cast(&dofs[0])); + dofs_per_cell, + &values[0], + reinterpret_cast(const_cast(&dofs[0]))); if (external_face == false) Matrix->SumIntoGlobalValues(dofs[i], dofs_per_cell, &values[dofs_per_cell], - reinterpret_cast(&dofs_neighbor[0])); + reinterpret_cast(const_cast(&dofs_neighbor[0]))); // And add into the residual @@ -1819,7 +1861,7 @@ void ConsLaw::assemble_face_term( // piece for each cell/face. We keep track of // the norm of the resdual for the Newton iteration. template -void ConsLaw::assemble_system (double &res_norm) +void ConservationLaw::assemble_system () { const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; @@ -1912,10 +1954,7 @@ void ConsLaw::assemble_system (double &res_norm) // asssemble the cell term. cell->get_dof_indices (dofs); - cell_diameter = cell->diameter(); - - assemble_cell_term(fe_v, - dofs); + assemble_cell_term(fe_v, dofs); // We use the DG style loop through faces // to determine if we need to apply a @@ -1926,7 +1965,6 @@ void ConsLaw::assemble_system (double &res_norm) // First we set the face // iterator typename DoFHandler::face_iterator face=cell->face(face_no); - face_diameter = face->diameter(); if (face->at_boundary()) { @@ -1941,13 +1979,13 @@ void ConsLaw::assemble_system (double &res_norm) // 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()); + assemble_face_term(face_no, fe_v_face, + fe_v_face, + dofs, + dofs, + true, + face->boundary_indicator(), + face->diameter()); } else { @@ -1977,8 +2015,6 @@ void ConsLaw::assemble_system (double &res_norm) neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no); - face_diameter = neighbor_child->diameter(); // working on subface - Assert (neighbor_child->face(neighbor2) == face->child(subface_no), ExcInternalError()); Assert (!neighbor_child->has_children(), ExcInternalError()); @@ -1989,14 +2025,13 @@ void ConsLaw::assemble_system (double &res_norm) // 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); - + 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())'' @@ -2031,14 +2066,13 @@ void ConsLaw::assemble_system (double &res_norm) 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); - + assemble_face_term(face_no, fe_v_face, + fe_v_subface_neighbor, + dofs, + dofs_neighbor, + false, + numbers::invalid_unsigned_int, + face->diameter()); } } @@ -2052,42 +2086,13 @@ void ConsLaw::assemble_system (double &res_norm) // Notify Epetra that the matrix is done. Matrix->FillComplete(); - - - // Compute the nonlinear residual. - res_norm = right_hand_side.l2_norm(); - -} - - // @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 () -{ - // First we need to distribute the - // DoFs. - dof_handler.clear(); - dof_handler.distribute_dofs (fe); - - // Size all of the fields. - solution.reinit (dof_handler.n_dofs()); - nlsolution.reinit (dof_handler.n_dofs()); - predictor.reinit (dof_handler.n_dofs()); - dsolution.reinit (dof_handler.n_dofs()); - right_hand_side.reinit (dof_handler.n_dofs()); - indicator.reinit(triangulation.n_active_cells()); } // @sect3{Setup System} // We call this function to build the sparsity // and the matrix. template -void ConsLaw::setup_system () +void ConservationLaw::setup_system () { // The DoFs of a cell are coupled @@ -2163,14 +2168,15 @@ void ConsLaw::setup_system () // Actually solve the linear system, using either // Aztec or Amesos. template -void ConsLaw::solve (Vector &dsolution, int &niter, double &lin_residual) +std::pair +ConservationLaw::solve (Vector &newton_update) { // 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 x(View, *Map, newton_update.begin()); Epetra_Vector b(View, *Map, right_hand_side.begin()); // The Direct option selects the Amesos solver. @@ -2211,14 +2217,14 @@ void ConsLaw::solve (Vector &dsolution, int &niter, double &lin_res if (parameters.output == Parameters::Solver::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 (parameters.solver == Parameters::Solver::gmres) { + return std::make_pair (0, 0); + } + else if (parameters.solver == Parameters::Solver::gmres) + { // For the iterative solvers, we use Aztec. AztecOO Solver; @@ -2254,17 +2260,21 @@ void ConsLaw::solve (Vector &dsolution, int &niter, double &lin_res // Run the solver iteration. Collect the number // of iterations and the residual. Solver.Iterate(parameters.max_iterations, parameters.linear_residual); - niter = Solver.NumIters(); - lin_residual = Solver.TrueResidual(); + + return std::make_pair (Solver.NumIters(), + Solver.TrueResidual()); } + + Assert (false, ExcNotImplemented()); + return std::make_pair (0,0); } // Loop and assign a value for refinement. We // simply use the density squared, which selects // shocks with some success. template -void ConsLaw::estimate() { - +void ConservationLaw::compute_refinement_indicator () +{ const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; std::vector dofs (dofs_per_cell); UpdateFlags update_flags = update_values @@ -2307,7 +2317,7 @@ void ConsLaw::estimate() { } template -void ConsLaw::refine_grid () +void ConservationLaw::refine_grid () { SolutionTransfer soltrans(dof_handler); @@ -2334,12 +2344,12 @@ void ConsLaw::refine_grid () } } - // The following code prolongs the solution + // The following code prolongs the old_solution // to the new grid and carries out the refinement. std::vector > interp_in; std::vector > interp_out; - interp_in.push_back(solution); + interp_in.push_back(old_solution); interp_in.push_back(predictor); triangulation.prepare_coarsening_and_refinement(); @@ -2351,44 +2361,45 @@ void ConsLaw::refine_grid () dof_handler.distribute_dofs (fe); { - Vector new_solution(1); + Vector new_old_solution(1); Vector new_predictor(1); - interp_out.push_back(new_solution); + interp_out.push_back(new_old_solution); interp_out.push_back(new_predictor); interp_out[0].reinit(dof_handler.n_dofs()); interp_out[1].reinit(dof_handler.n_dofs()); } soltrans.interpolate(interp_in, interp_out); - - // Let the vector delete a very small vector - solution.reinit(1); - predictor.reinit(1); - solution.swap(interp_out[0]); - predictor.swap(interp_out[1]); + + old_solution.reinit (interp_out[0].size()); + old_solution = interp_out[0]; + + predictor.reinit (interp_out[1].size()); + predictor = interp_out[1]; // resize these vectors for the new grid. - nlsolution.reinit(dof_handler.n_dofs()); - nlsolution = solution; - dsolution.reinit (dof_handler.n_dofs()); + current_solution.reinit(dof_handler.n_dofs()); + current_solution = old_solution; right_hand_side.reinit (dof_handler.n_dofs()); indicator.reinit(triangulation.n_active_cells()); } + + template -void ConsLaw::output_results (const unsigned int cycle) const +void ConservationLaw::output_results () const { typename EulerEquations::Postprocessor postprocessor (parameters.schlieren_plot); DataOut data_out; data_out.attach_dof_handler (dof_handler); - std::vector solution_names (dim, "momentum"); - solution_names.push_back ("density"); - solution_names.push_back ("energy_density"); + std::vector old_solution_names (dim, "momentum"); + old_solution_names.push_back ("density"); + old_solution_names.push_back ("energy_density"); std::vector data_component_interpretation @@ -2398,43 +2409,34 @@ void ConsLaw::output_results (const unsigned int cycle) const data_component_interpretation .push_back (DataComponentInterpretation::component_is_scalar); - data_out.add_data_vector (solution, solution_names, + data_out.add_data_vector (old_solution, old_solution_names, DataOut::type_dof_data, data_component_interpretation); - data_out.add_data_vector (solution, postprocessor); + data_out.add_data_vector (old_solution, postprocessor); data_out.add_data_vector (indicator, "error"); data_out.build_patches (); - std::string filename = "solution-" + - Utilities::int_to_string (cycle, 3) + + static unsigned int output_file_number = 0; + std::string filename = "old_solution-" + + Utilities::int_to_string (output_file_number, 3) + ".vtk"; std::ofstream output (filename.c_str()); data_out.write_vtk (output); + + ++output_file_number; } - // 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. template -void ConsLaw::run () +void ConservationLaw::run () { // Open and load the mesh. @@ -2448,13 +2450,25 @@ void ConsLaw::run () grid_in.read_ucd(input_file); } - unsigned int nstep = 0; - // Initialize fields and matrices. - initialize_system (); + // First we need to distribute the + // DoFs. + dof_handler.clear(); + dof_handler.distribute_dofs (fe); + + // Size all of the fields. + old_solution.reinit (dof_handler.n_dofs()); + current_solution.reinit (dof_handler.n_dofs()); + predictor.reinit (dof_handler.n_dofs()); + right_hand_side.reinit (dof_handler.n_dofs()); + indicator.reinit(triangulation.n_active_cells()); + setup_system(); - initialize(); - predictor = solution; + + VectorTools::interpolate(dof_handler, + parameters.initial_conditions, old_solution); + current_solution = old_solution; + predictor = old_solution; // Initial refinement. We apply the ic, // estimate, refine, and repeat until @@ -2462,23 +2476,28 @@ void ConsLaw::run () if (parameters.do_refine == true) for (unsigned int i = 0; i < parameters.shock_levels; i++) { - estimate(); + compute_refinement_indicator(); refine_grid(); setup_system(); - initialize(); - predictor = solution; + + VectorTools::interpolate(dof_handler, + parameters.initial_conditions, old_solution); + current_solution = old_solution; + predictor = old_solution; } - output_results (nstep); + output_results (); // Determine when we will output next. - double next_output = T + parameters.output_step; + double time = 0; + double next_output = time + parameters.output_step; // @sect4{Main time stepping loop} - predictor = solution; - while (T < parameters.final_time) + predictor = old_solution; + Vector newton_update (dof_handler.n_dofs()); + while (time < parameters.final_time) { - std::cout << "T=" << T << ", "; + std::cout << "T=" << time << ", "; std::cout << " Number of active cells: " @@ -2491,8 +2510,6 @@ void ConsLaw::run () << std::endl; bool nonlin_done = false; - double res_norm; - int lin_iter; // Print some relevant information during the // Newton iteration. @@ -2501,40 +2518,45 @@ void ConsLaw::run () const unsigned int max_nonlin = 7; unsigned int nonlin_iter = 0; - double lin_res; // @sect5{Newton iteration} - nlsolution = predictor; + current_solution = predictor; while (!nonlin_done) { - lin_iter = 0; - Matrix->PutScalar(0); Matrix->FillComplete(); right_hand_side = 0; - assemble_system (res_norm); + assemble_system (); + // Flash a star to the screen so one can // know when the assembly has stopped and the linear - // solution is starting. + // old_solution is starting. std::cout << "* " << std::flush; // Test against a (hardcoded) nonlinear tolderance. // Do not solve the linear system at the last step // (since it would be a waste). - if (fabs(res_norm) < 1e-10) { - nonlin_done = true; - } else { + const double res_norm = right_hand_side.l2_norm(); + if (std::fabs(res_norm) < 1e-10) + { + nonlin_done = true; + std::printf("%-16.3e (converged)\n", res_norm); + } + else + { // Solve the linear system and update with the // delta. - dsolution = 0; - solve (dsolution, lin_iter, lin_res); - nlsolution.add(1.0, dsolution); - } - - // Print the residuals. - std::printf("%-16.3e %04d %-5.2e\n", - res_norm, lin_iter, lin_res); + newton_update = 0; + + std::pair convergence + = solve (newton_update); + + current_solution.add(1.0, newton_update); + + std::printf("%-16.3e %04d %-5.2e\n", + res_norm, convergence.first, convergence.second); + } ++nonlin_iter; @@ -2543,27 +2565,38 @@ void ConsLaw::run () } // Various post convergence tasks. - compute_predictor(); - solution = nlsolution; + // 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 old_solution. This + // simple time extrapolator does the job. + predictor = current_solution; + predictor.sadd(3/2.0, -1/2.0, old_solution); - estimate(); + old_solution = current_solution; - T += parameters.time_step; + compute_refinement_indicator(); + + time += parameters.time_step; // Output if it is time. - if (parameters.output_step < 0) { - output_results (++nstep); - } else if (T >= next_output) { - output_results (++nstep); - next_output += parameters.output_step; - } + if (parameters.output_step < 0) + output_results (); + else if (time >= next_output) + { + output_results (); + next_output += parameters.output_step; + } // Refine, if refinement is selected. if (parameters.do_refine == true) { refine_grid(); setup_system(); + + newton_update.reinit (dof_handler.n_dofs()); } } } @@ -2584,7 +2617,7 @@ int main (int argc, char *argv[]) try { - ConsLaw<2> cons (argv[1]); + ConservationLaw<2> cons (argv[1]); cons.run (); } catch (std::exception &exc) -- 2.39.5