From e09b338e6c50f1054d6773068a283a18cd21469b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 19 Aug 2009 18:08:53 +0000 Subject: [PATCH] Fix determination of c_R in the program, and document our value. git-svn-id: https://svn.dealii.org/trunk@19314 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-31/doc/results.dox | 30 +++- deal.II/examples/step-31/step-31.cc | 191 ++++++++++++----------- 2 files changed, 124 insertions(+), 97 deletions(-) diff --git a/deal.II/examples/step-31/doc/results.dox b/deal.II/examples/step-31/doc/results.dox index e991bb03b1..757a9415ae 100644 --- a/deal.II/examples/step-31/doc/results.dox +++ b/deal.II/examples/step-31/doc/results.dox @@ -325,7 +325,29 @@ In all of these cases, we will have to expect that the correct choice of each value depends on that of the others, and most likely also on the space dimension and polynomial degree of the finite element used for the temperature. Below we'll discuss a few numerical experiments to choose -constants. +constants $c_k$ and $\beta$. + +Below, we will not discuss the choice of $c_R$. In the program, we set +it to $c_R=2^{\frac{4-2\alpha}{d}}$. The reason for this value is a +bit complicated and has more to do with the history of the program +than reasoning: while the correct formula for the global scaling +parameter $c(\mathbf{u},T)$ is shown above, the program (including the +version shipped with deal.II 6.2) initially had a bug in that we +computed +$c(\mathbf{u},T) = + \|\mathbf{u}\|_{L^\infty(\Omega)} \ \mathrm{var}(T) + \ \frac{1}{|\mathrm{diam}(\Omega)|^{\alpha-2}}$ instead, where +we had set the scaling parameter to one. Since we only computed on the +unit square/cube where $\mathrm{diam}(\Omega)=2^{1/d}$, this was +entirely equivalent to using the correct formula with +$c_R=\left(2^{1/d}\right)^{4-2\alpha}=2^{\frac{4-2\alpha}{d}}$. Since +this value for $c_R$ appears to work just fine for the current +program, we corrected the formula in the program and set $c_R$ to a +value that reproduces exactly the results we had before. We will, +however, revisit this issue again in @ref step_32 "step-32". + +Now, however, back to the discussion of what values of $c_k$ and +$\beta$ to choose:

Choosing ck and beta

@@ -367,7 +389,11 @@ $ = \beta \|\mathbf{u}\|_{L^\infty(K)} h_K -$ to eliminate the effect of of the constant $c_R$. We then run the program +$ to eliminate the effect of of the constant $c_R$ (we know that +solutions are stable by using above $\nu(T)$ as an artificial +viscosity, but that we can improve things -- i.e. make the solution +sharper -- by using the more complicated formula for this artificial +viscosity). We then run the program for different values $c_k,\beta$ and observe maximal and minimal temperatures in the domain. What we expect to see is this: If we choose the time step too big (i.e. choose a $c_k$ bigger than theoretically allowed) then we will get diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 6125b21223..1d6e147a8e 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -127,7 +127,7 @@ using namespace dealii; // circles (or spheres in 3d) somewhere at // the bottom of the domain, as explained in // the introduction, and zero outside. - // + // // Finally, or maybe firstly, at the top of // this namespace, we define the various // material constants we need ($\eta,\kappa$, @@ -195,16 +195,16 @@ namespace EquationData { Assert (component == 0, ExcMessage ("Invalid operation for a scalar function.")); - + Assert ((dim==2) || (dim==3), ExcNotImplemented()); - + static const Point source_centers[3] = { (dim == 2 ? Point(.3,.1) : Point(.3,.5,.1)), (dim == 2 ? Point(.45,.1) : Point(.45,.5,.1)), (dim == 2 ? Point(.75,.1) : Point(.75,.5,.1)) }; static const double source_radius = (dim == 2 ? 1./32 : 1./8); - + return ((source_centers[0].distance (p) < source_radius) || (source_centers[1].distance (p) < source_radius) @@ -392,7 +392,7 @@ namespace LinearSolvers // \end{array}\right) = // \left(\begin{array}{cc} I & // A^{-1} B^T \\ 0 & I - // \end{array}\right), + // \end{array}\right), // @f} // which indeed is very simple. A GMRES // solver based on exact matrices would @@ -408,7 +408,7 @@ namespace LinearSolvers // Using general block preconditioners", // SIAM J. Numer. Anal., 31 (1994), // pp. 1352-1367). - // + // // Replacing P by $\tilde{P}$ // keeps that spirit alive: the product // $P^{-1} A$ will still be close to a @@ -444,7 +444,7 @@ namespace LinearSolvers // iterations for each outer // solver step (using the AMG // preconditioner). - // + // // Having the above explanations in mind, // we define a preconditioner class with // a vmult functionality, @@ -452,7 +452,7 @@ namespace LinearSolvers // interaction with the usual solver // functions further below in the program // code. - // + // // First the declarations. These are // similar to the definition of the Schur // complement in step-20, with the @@ -627,7 +627,7 @@ class BoussinesqFlowProblem TrilinosWrappers::BlockVector stokes_rhs; - const unsigned int temperature_degree; + const unsigned int temperature_degree; FE_Q temperature_fe; DoFHandler temperature_dof_handler; ConstraintMatrix temperature_constraints; @@ -658,7 +658,7 @@ class BoussinesqFlowProblem // @sect3{BoussinesqFlowProblem class implementation} // @sect4{BoussinesqFlowProblem::BoussinesqFlowProblem} - // + // // The constructor of this class is an // extension of the constructor in // step-22. We need to add the various @@ -785,7 +785,7 @@ double BoussinesqFlowProblem::get_maximal_velocity () const double max_velocity = 0; const FEValuesExtractors::Vector velocities (0); - + typename DoFHandler::active_cell_iterator cell = stokes_dof_handler.begin_active(), endc = stokes_dof_handler.end(); @@ -895,7 +895,7 @@ BoussinesqFlowProblem::get_extrapolated_temperature_range () const for (unsigned int q=0; q::get_extrapolated_temperature_range () const max_temperature = std::max (max_temperature, temperature); } } - + return std::make_pair(min_temperature, max_temperature); - } + } } @@ -991,25 +991,25 @@ compute_viscosity (const std::vector &old_temperature, { const double beta = 0.015 * dim; const double alpha = 1; - + if (global_u_infty == 0) return 5e-3 * cell_diameter; - + const unsigned int n_q_points = old_temperature.size(); - + double max_residual = 0; double max_velocity = 0; - + for (unsigned int q=0; q < n_q_points; ++q) { const Tensor<1,dim> u = (old_velocity_values[q] + old_old_velocity_values[q]) / 2; - + const double dT_dt = (old_temperature[q] - old_old_temperature[q]) / old_time_step; const double u_grad_T = u * (old_temperature_grads[q] + old_old_temperature_grads[q]) / 2; - + const double kappa_Delta_T = EquationData::kappa * (old_temperature_laplacians[q] + old_old_temperature_laplacians[q]) / 2; @@ -1022,8 +1022,9 @@ compute_viscosity (const std::vector &old_temperature, max_residual = std::max (residual, max_residual); max_velocity = std::max (std::sqrt (u*u), max_velocity); } - - const double global_scaling = global_u_infty * global_T_variation * + + const double c_R = std::pow (2., (4.-2*alpha)/dim); + const double global_scaling = c_R * global_u_infty * global_T_variation * std::pow(global_Omega_diameter, alpha - 2.); return (beta * @@ -1036,7 +1037,7 @@ compute_viscosity (const std::vector &old_temperature, // @sect4{BoussinesqFlowProblem::setup_dofs} - // + // // This is the function that sets up the // DoFHandler objects we have here (one for // the Stokes part and one for the @@ -1045,7 +1046,7 @@ compute_viscosity (const std::vector &old_temperature, // for the linear algebra in this // program. Its basic operations are similar // to what we do in step-22. - // + // // The body of the function first // enumerates all degrees of freedom for // the Stokes and temperature systems. For @@ -1066,7 +1067,7 @@ compute_viscosity (const std::vector &old_temperature, // renumbering, but its costs are low // compared to the velocity portion, so the // additional work does not pay off. - // + // // We then proceed with the generation of the // hanging node constraints that arise from // adaptive grid refinement for both @@ -1092,11 +1093,11 @@ void BoussinesqFlowProblem::setup_dofs () { std::vector stokes_sub_blocks (dim+1,0); stokes_sub_blocks[dim] = 1; - + { stokes_dof_handler.distribute_dofs (stokes_fe); DoFRenumbering::component_wise (stokes_dof_handler, stokes_sub_blocks); - + stokes_constraints.clear (); DoFTools::make_hanging_node_constraints (stokes_dof_handler, stokes_constraints); @@ -1119,7 +1120,7 @@ void BoussinesqFlowProblem::setup_dofs () std::vector stokes_dofs_per_block (2); DoFTools::count_dofs_per_block (stokes_dof_handler, stokes_dofs_per_block, stokes_sub_blocks); - + const unsigned int n_u = stokes_dofs_per_block[0], n_p = stokes_dofs_per_block[1], n_T = temperature_dof_handler.n_dofs(); @@ -1135,7 +1136,7 @@ void BoussinesqFlowProblem::setup_dofs () << " (" << n_u << '+' << n_p << '+'<< n_T <<')' << std::endl << std::endl; - + // The next step is to create the sparsity // pattern for the Stokes and temperature // system matrices as well as the @@ -1150,7 +1151,7 @@ void BoussinesqFlowProblem::setup_dofs () // class would consume too much memory when // used in three spatial dimensions as we // intend to do for this program. - // + // // So, we first release the memory stored // in the matrices, then set up an object // of type @@ -1215,7 +1216,7 @@ void BoussinesqFlowProblem::setup_dofs () stokes_matrix.clear (); BlockCompressedSimpleSparsityPattern csp (2,2); - + csp.block(0,0).reinit (n_u, n_u); csp.block(0,1).reinit (n_u, n_p); csp.block(1,0).reinit (n_p, n_u); @@ -1244,12 +1245,12 @@ void BoussinesqFlowProblem::setup_dofs () stokes_preconditioner_matrix.clear (); BlockCompressedSimpleSparsityPattern csp (2,2); - + csp.block(0,0).reinit (n_u, n_u); csp.block(0,1).reinit (n_u, n_p); csp.block(1,0).reinit (n_p, n_u); csp.block(1,1).reinit (n_p, n_p); - + csp.collect_sizes (); Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); @@ -1293,7 +1294,7 @@ void BoussinesqFlowProblem::setup_dofs () temperature_stiffness_matrix.clear (); temperature_matrix.clear (); - CompressedSimpleSparsityPattern csp (n_T, n_T); + CompressedSimpleSparsityPattern csp (n_T, n_T); DoFTools::make_sparsity_pattern (temperature_dof_handler, csp, temperature_constraints, false); @@ -1324,7 +1325,7 @@ void BoussinesqFlowProblem::setup_dofs () // @sect4{BoussinesqFlowProblem::assemble_stokes_preconditioner} - // + // // This function assembles the matrix we use // for preconditioning the Stokes // system. What we need are a vector Laplace @@ -1405,7 +1406,7 @@ BoussinesqFlowProblem::assemble_stokes_preconditioner () phi_grad_u[k] = stokes_fe_values[velocities].gradient(k,q); phi_p[k] = stokes_fe_values[pressure].value (k, q); } - + for (unsigned int i=0; i::assemble_stokes_preconditioner () // @sect4{BoussinesqFlowProblem::build_stokes_preconditioner} - // + // // This function generates the inner // preconditioners that are going to be used // for the Schur complement block @@ -1441,7 +1442,7 @@ BoussinesqFlowProblem::assemble_stokes_preconditioner () // task is to call // assemble_stokes_preconditioner // to generate the preconditioner matrices. - // + // // Next, we set up the preconditioner for // the velocity-velocity matrix // A. As explained in the @@ -1475,18 +1476,18 @@ BoussinesqFlowProblem::build_stokes_preconditioner () { if (rebuild_stokes_preconditioner == false) return; - + std::cout << " Rebuilding Stokes preconditioner..." << std::flush; - + assemble_stokes_preconditioner (); - + Amg_preconditioner = std_cxx1x::shared_ptr (new TrilinosWrappers::PreconditionAMG()); std::vector > constant_modes; std::vector velocity_components (dim+1,true); velocity_components[dim] = false; - DoFTools::extract_constant_modes (stokes_dof_handler, velocity_components, + DoFTools::extract_constant_modes (stokes_dof_handler, velocity_components, constant_modes); TrilinosWrappers::PreconditionAMG::AdditionalData amg_data; amg_data.constant_modes = constant_modes; @@ -1526,7 +1527,7 @@ BoussinesqFlowProblem::build_stokes_preconditioner () // this data set, we then initialize the // preconditioner with the matrix we want // it to apply to. - // + // // Finally, we also initialize the // preconditioner for the inversion of // the pressure mass matrix. This matrix @@ -1564,7 +1565,7 @@ BoussinesqFlowProblem::build_stokes_preconditioner () // @sect4{BoussinesqFlowProblem::assemble_stokes_system} - // + // // The time lag scheme we use for advancing // the coupled Stokes-temperature system // forces us to split up the assembly (and @@ -1575,7 +1576,7 @@ BoussinesqFlowProblem::build_stokes_preconditioner () // and right hand sides for the temperature // dofs, which depends on the result of the // linear system for the velocity. - // + // // This function is called at the beginning // of each time step. In the first time step // or if the mesh has changed, indicated by @@ -1586,7 +1587,7 @@ BoussinesqFlowProblem::build_stokes_preconditioner () // not necessary and all we need to do is // assemble the right hand side vector which // changes in each time step. - // + // // Regarding the technical details of // implementation, not much has changed from // step-22. We reset matrix and vector, @@ -1617,7 +1618,7 @@ BoussinesqFlowProblem::build_stokes_preconditioner () // FEValues objects to ensure that we get // matching information when we loop over the // quadrature points of the two objects. - // + // // The declarations proceed with some // shortcuts for array sizes, the creation // of the local matrix and right hand side @@ -1644,7 +1645,7 @@ void BoussinesqFlowProblem::assemble_stokes_system () update_gradients : UpdateFlags(0))); - + FEValues temperature_fe_values (temperature_fe, quadrature_formula, update_values); @@ -1663,7 +1664,7 @@ void BoussinesqFlowProblem::assemble_stokes_system () // term in the right hand side of the // momentum equation. Let's call this vector // old_solution_values. - // + // // The set of vectors we create next hold // the evaluations of the basis functions // as well as their gradients and @@ -1674,7 +1675,7 @@ void BoussinesqFlowProblem::assemble_stokes_system () // each time it is needed is an // optimization to accelerate the assembly // process, see step-22 for details. - // + // // The last two declarations are used to // extract the individual blocks // (velocity, pressure, temperature) from @@ -1713,16 +1714,16 @@ void BoussinesqFlowProblem::assemble_stokes_system () endc = stokes_dof_handler.end(); typename DoFHandler::active_cell_iterator temperature_cell = temperature_dof_handler.begin_active(); - + for (; cell!=endc; ++cell, ++temperature_cell) { stokes_fe_values.reinit (cell); temperature_fe_values.reinit (temperature_cell); - + local_matrix = 0; local_rhs = 0; - temperature_fe_values.get_function_values (old_temperature_solution, + temperature_fe_values.get_function_values (old_temperature_solution, old_temperature_values); for (unsigned int q=0; q::assemble_stokes_system () // inner products. As shown in // step-22 this helps accelerate // assembly. - // + // // Once this is done, we start the // loop over the rows and columns // of the local matrix and feed the @@ -1772,7 +1773,7 @@ void BoussinesqFlowProblem::assemble_stokes_system () - phi_p[i] * div_phi_u[j]) * stokes_fe_values.JxW(q); - const Point gravity = -( (dim == 2) ? (Point (0,1)) : + const Point gravity = -( (dim == 2) ? (Point (0,1)) : (Point (0,0,1)) ); for (unsigned int i=0; i::assemble_stokes_system () // @sect4{BoussinesqFlowProblem::assemble_temperature_matrix} - // + // // This function assembles the matrix in // the temperature equation. The // temperature matrix consists of two @@ -1838,7 +1839,7 @@ void BoussinesqFlowProblem::assemble_stokes_system () // will then sum up the matrix plus the // stiffness matrix times the time step // size once we know the actual time step. - // + // // So the details for this first step are // very simple. In case we need to // rebuild the matrix (i.e., the mesh has @@ -1853,10 +1854,10 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () { if (rebuild_temperature_matrices == false) return; - + temperature_mass_matrix = 0; temperature_stiffness_matrix = 0; - + QGauss quadrature_formula (temperature_degree+2); FEValues temperature_fe_values (temperature_fe, quadrature_formula, update_values | update_gradients | @@ -1896,7 +1897,7 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () local_stiffness_matrix = 0; temperature_fe_values.reinit (cell); - + for (unsigned int q=0; q::assemble_temperature_matrix () grad_phi_T[k] = temperature_fe_values.shape_grad (k,q); phi_T[k] = temperature_fe_values.shape_value (k, q); } - + for (unsigned int i=0; i::assemble_temperature_matrix () temperature_fe_values.JxW(q)); } } - + cell->get_dof_indices (local_dof_indices); temperature_constraints.distribute_local_to_global (local_mass_matrix, @@ -1928,7 +1929,7 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () local_dof_indices, temperature_stiffness_matrix); } - + rebuild_temperature_matrices = false; } @@ -1980,9 +1981,9 @@ void BoussinesqFlowProblem:: temperature_matrix.copy_from (temperature_mass_matrix); temperature_matrix.add (time_step, temperature_stiffness_matrix); } - + temperature_rhs = 0; - + const QGauss quadrature_formula(temperature_degree+2); FEValues temperature_fe_values (temperature_fe, quadrature_formula, update_values | @@ -2033,12 +2034,12 @@ void BoussinesqFlowProblem:: std::vector phi_T (dofs_per_cell); std::vector > grad_phi_T (dofs_per_cell); - + const std::pair global_T_range = get_extrapolated_temperature_range(); const FEValuesExtractors::Vector velocities (0); - + // Now, let's start the loop over all cells // in the triangulation. Again, we need two // cell iterators that walk in parallel @@ -2084,12 +2085,12 @@ void BoussinesqFlowProblem:: old_temperature_grads); temperature_fe_values.get_function_gradients (old_old_temperature_solution, old_old_temperature_grads); - + temperature_fe_values.get_function_laplacians (old_temperature_solution, old_temperature_laplacians); temperature_fe_values.get_function_laplacians (old_old_temperature_solution, old_old_temperature_laplacians); - + temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(), gamma_values); @@ -2133,7 +2134,7 @@ void BoussinesqFlowProblem:: maximal_velocity, global_T_range.second - global_T_range.first, cell->diameter()); - + for (unsigned int q=0; q:: old_temperature_values[q]); const Tensor<1,dim> ext_grad_T - = (use_bdf2_scheme ? + = (use_bdf2_scheme ? (old_temperature_grads[q] * - (1+time_step/old_time_step) + (1+time_step/old_time_step) - old_old_temperature_grads[q] * - time_step / old_time_step) + time_step / old_time_step) : old_temperature_grads[q]); - + const Tensor<1,dim> extrapolated_u - = (use_bdf2_scheme ? - (old_velocity_values[q] * (1+time_step/old_time_step) - + = (use_bdf2_scheme ? + (old_velocity_values[q] * (1+time_step/old_time_step) - old_old_velocity_values[q] * time_step/old_time_step) : old_velocity_values[q]); @@ -2184,7 +2185,7 @@ void BoussinesqFlowProblem:: * temperature_fe_values.JxW(q); } - + cell->get_dof_indices (local_dof_indices); temperature_constraints.distribute_local_to_global (local_rhs, local_dof_indices, @@ -2219,7 +2220,7 @@ void BoussinesqFlowProblem:: // only pointers, so we use * // to pass down the actual preconditioner // objects. - // + // // Once the preconditioner is ready, we // create a GMRES solver for the block // system. Since we are working with @@ -2331,7 +2332,7 @@ void BoussinesqFlowProblem::solve () std::cout << " " << "Time step: " << time_step << std::endl; - + temperature_solution = old_temperature_solution; // Next we set up the temperature system @@ -2404,7 +2405,7 @@ void BoussinesqFlowProblem::solve () max_temperature = std::max (max_temperature, temperature_solution(i)); } - + std::cout << " Temperature range: " << min_temperature << ' ' << max_temperature << std::endl; @@ -2414,7 +2415,7 @@ void BoussinesqFlowProblem::solve () // @sect4{BoussinesqFlowProblem::output_results} - // + // // This function writes the solution to a VTK // output file for visualization, which is // done every tenth time step. This is @@ -2441,7 +2442,7 @@ void BoussinesqFlowProblem::solve () // that ensures that we got all the dofs from // both Stokes and temperature even in the // combined system. - // + // // Next, we create a vector that will collect // the actual solution values. Since this // vector is only going to be used for @@ -2462,7 +2463,7 @@ void BoussinesqFlowProblem::output_results () const Assert (joint_dof_handler.n_dofs() == stokes_dof_handler.n_dofs() + temperature_dof_handler.n_dofs(), ExcInternalError()); - + Vector joint_solution (joint_dof_handler.n_dofs()); // Unfortunately, there is no @@ -2495,7 +2496,7 @@ void BoussinesqFlowProblem::output_results () const // between global vector and local dofs // looks like on the present cell, which // concludes this tedious work. - // + // // There's one thing worth remembering // when looking at the output: In our // algorithm, we first solve for the @@ -2515,7 +2516,7 @@ void BoussinesqFlowProblem::output_results () const std::vector local_joint_dof_indices (joint_fe.dofs_per_cell); std::vector local_stokes_dof_indices (stokes_fe.dofs_per_cell); std::vector local_temperature_dof_indices (temperature_fe.dofs_per_cell); - + typename DoFHandler::active_cell_iterator joint_cell = joint_dof_handler.begin_active(), joint_endc = joint_dof_handler.end(), @@ -2550,7 +2551,7 @@ void BoussinesqFlowProblem::output_results () const } } } - + // Next, we proceed as we've done in // step-22. We create solution names // (that are going to appear in the @@ -2603,7 +2604,7 @@ void BoussinesqFlowProblem::output_results () const // @sect4{BoussinesqFlowProblem::refine_mesh} - // + // // This function takes care of the adaptive // mesh refinement. The three tasks this // function performs is to first find out @@ -2669,7 +2670,7 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) GridRefinement::refine_and_coarsen_fixed_fraction (triangulation, estimated_error_per_cell, 0.8, 0.1); - if (triangulation.n_levels() > max_grid_level) + if (triangulation.n_levels() > max_grid_level) for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(max_grid_level); cell != triangulation.end(); ++cell) @@ -2766,7 +2767,7 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) // @sect4{BoussinesqFlowProblem::run} - // + // // This function performs all the // essential steps in the Boussinesq // program. It starts by setting up a @@ -2792,7 +2793,7 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) // beginning at the // start_time_iteration // label. - // + // // Before we start, we project the // initial values to the grid and // obtain the first data for the @@ -2815,7 +2816,7 @@ void BoussinesqFlowProblem::run () setup_dofs(); unsigned int pre_refinement_step = 0; - + start_time_iteration: VectorTools::project (temperature_dof_handler, @@ -2823,10 +2824,10 @@ void BoussinesqFlowProblem::run () QGauss(temperature_degree+2), EquationData::TemperatureInitialValues(), old_temperature_solution); - + timestep_number = 0; time_step = old_time_step = 0; - + double time = 0; do @@ -2864,7 +2865,7 @@ void BoussinesqFlowProblem::run () output_results (); std::cout << std::endl; - + if ((timestep_number == 0) && (pre_refinement_step < n_pre_refinement_steps)) { @@ -2891,7 +2892,7 @@ void BoussinesqFlowProblem::run () // @sect3{The main function} - // + // // The main function looks almost the same // as in all other programs. The only // difference is that Trilinos wants to get -- 2.39.5