Ifpack_ex_ScalarLaplacian_FEM.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                IFPACK: Robust Algebraic Preconditioning Package
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "Ifpack_ConfigDefs.h"
00030 #include "Galeri_ConfigDefs.h"
00031 #ifdef HAVE_MPI
00032 #include "mpi.h"
00033 #include "Epetra_MpiComm.h"
00034 #else
00035 #include "Epetra_SerialComm.h"
00036 #endif
00037 #include "Epetra_FECrsMatrix.h"
00038 #include "Epetra_FEVector.h"
00039 
00040 #include "Galeri_core_Object.h"
00041 #include "Galeri_core_Workspace.h"
00042 #include "Galeri_grid_Loadable.h"
00043 #include "Galeri_grid_Generator.h"
00044 #include "Galeri_quadrature_Hex.h"
00045 #include "Galeri_problem_ScalarLaplacian.h"
00046 #include "Galeri_viz_MEDIT.h"
00047 
00048 #include "AztecOO.h"
00049 #include "Ifpack.h"
00050 #include "Ifpack_AdditiveSchwarz.h"
00051 
00052 using namespace Galeri;
00053 
00054 class Laplacian 
00055 {
00056   public:
00057     static inline double 
00058     getElementLHS(const double& x, 
00059                   const double& y, 
00060                   const double& z,
00061                   const double& phi_i,
00062                   const double& phi_i_derx, 
00063                   const double& phi_i_dery,
00064                   const double& phi_i_derz,
00065                   const double& phi_j,
00066                   const double& phi_j_derx,
00067                   const double& phi_j_dery,
00068                   const double& phi_j_derz)
00069     {
00070       return(phi_i_derx * phi_j_derx + 
00071              phi_i_dery * phi_j_dery + 
00072              phi_i_derz * phi_j_derz);
00073     }
00074 
00075     static inline double
00076     getElementRHS(const double& x,
00077                   const double& y,
00078                   const double& z,
00079                   const double& phi_i)
00080 
00081     {
00082       return(-getExactSolution('f', x, y, z) * phi_i);
00083     }
00084 
00085     static inline double
00086     getBoundaryValue(const double& x, const double& y, const double& z)
00087     {
00088       return(getExactSolution('f', x, y, z));
00089     }
00090 
00091     static inline char
00092     getBoundaryType(const int ID, const double& x, const double& y, const double& z)
00093     {
00094       return('d');
00095     }
00096 
00097     static inline double 
00098     getExactSolution(const char& what, const double& x, 
00099                      const double& y, const double& z)
00100     {
00101       if (what == 'f')
00102         return(exp(x) + exp(y) + exp(z));
00103       else if (what == 'x')
00104         return(exp(x));
00105       else if (what == 'y')
00106         return(exp(y));
00107       else if (what == 'z')
00108         return(exp(z));
00109       else
00110         return(0.0);
00111     }
00112 };
00113 
00114 // =========== //
00115 // main driver //
00116 // =========== //
00117 
00118 int main(int argc, char *argv[])
00119 {
00120 #ifdef HAVE_MPI
00121   MPI_Init(&argc,&argv);
00122   Epetra_MpiComm comm(MPI_COMM_WORLD);
00123 #else
00124   Epetra_SerialComm comm;
00125 #endif
00126 
00127   Galeri::core::Workspace::setNumDimensions(3);
00128 
00129   Galeri::grid::Loadable domain, boundary;
00130 
00131   int numGlobalElementsX = 2 * comm.NumProc();
00132   int numGlobalElementsY = 2;
00133   int numGlobalElementsZ = 2;
00134 
00135   int mx = comm.NumProc();
00136   int my = 1;
00137   int mz = 1;
00138 
00139   Galeri::grid::Generator::
00140   getCubeWithHexs(comm, numGlobalElementsX, numGlobalElementsY, numGlobalElementsZ,
00141                   mx, my, mz, domain, boundary);
00142 
00143   Epetra_Map matrixMap(domain.getNumGlobalVertices(), 0, comm);
00144 
00145   Epetra_FECrsMatrix A(Copy, matrixMap, 0);
00146   Epetra_FEVector    LHS(matrixMap);
00147   Epetra_FEVector    RHS(matrixMap);
00148 
00149   Galeri::problem::ScalarLaplacian<Laplacian> problem("Hex", 1, 8);
00150 
00151   problem.integrate(domain, A, RHS);
00152 
00153   LHS.PutScalar(0.0);
00154 
00155   problem.imposeDirichletBoundaryConditions(boundary, A, RHS, LHS);
00156 
00157   // ============================================================ //
00158   // Solving the linear system is the next step, using the IFPACK //
00159   // factory. This is done by using the IFPACK factory, then      //
00160   // asking for IC preconditioner, and setting few parameters     //
00161   // using a Teuchos::ParameterList.                              //
00162   // ============================================================ //
00163   
00164   Ifpack Factory;
00165   Ifpack_Preconditioner* Prec = Factory.Create("IC", &A, 0);
00166 
00167   Teuchos::ParameterList list;
00168   
00169   list.set("fact: level-of-fill", 1);
00170   IFPACK_CHK_ERR(Prec->SetParameters(list));
00171   IFPACK_CHK_ERR(Prec->Initialize());
00172   IFPACK_CHK_ERR(Prec->Compute());
00173 
00174   Epetra_LinearProblem linearProblem(&A, &LHS, &RHS);
00175 
00176   AztecOO solver(linearProblem);
00177   solver.SetAztecOption(AZ_solver, AZ_cg);
00178   solver.SetPrecOperator(Prec);
00179   solver.Iterate(1550, 1e-9);
00180 
00181   // visualization using MEDIT -- a VTK module is avaiable as well
00182   Galeri::viz::MEDIT::write(domain, "sol", LHS);
00183 
00184   // now compute the norm of the solution
00185   problem.computeNorms(domain, LHS);
00186 
00187 #ifdef HAVE_MPI
00188   MPI_Finalize();
00189 #endif
00190 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:34 2011 for Ifpack Package Browser (Single Doxygen Collection) by  doxygen 1.6.3