Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_ex_ScalarLaplacian_FEM.cpp
Go to the documentation of this file.
1//@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41
42#include "Ifpack_ConfigDefs.h"
43#include "Galeri_ConfigDefs.h"
44#ifdef HAVE_MPI
45#include "mpi.h"
46#include "Epetra_MpiComm.h"
47#else
48#include "Epetra_SerialComm.h"
49#endif
50#include "Epetra_FECrsMatrix.h"
51#include "Epetra_FEVector.h"
52
53#include "Galeri_core_Object.h"
54#include "Galeri_core_Workspace.h"
55#include "Galeri_grid_Loadable.h"
56#include "Galeri_grid_Generator.h"
57#include "Galeri_quadrature_Hex.h"
58#include "Galeri_problem_ScalarLaplacian.h"
59#include "Galeri_viz_MEDIT.h"
60
61#include "AztecOO.h"
62#include "Ifpack.h"
64
65using namespace Galeri;
66
68{
69 public:
70 static inline double
71 getElementLHS(const double& x,
72 const double& y,
73 const double& z,
74 const double& phi_i,
75 const double& phi_i_derx,
76 const double& phi_i_dery,
77 const double& phi_i_derz,
78 const double& phi_j,
79 const double& phi_j_derx,
80 const double& phi_j_dery,
81 const double& phi_j_derz)
82 {
83 return(phi_i_derx * phi_j_derx +
84 phi_i_dery * phi_j_dery +
85 phi_i_derz * phi_j_derz);
86 }
87
88 static inline double
89 getElementRHS(const double& x,
90 const double& y,
91 const double& z,
92 const double& phi_i)
93
94 {
95 return(-getExactSolution('f', x, y, z) * phi_i);
96 }
97
98 static inline double
99 getBoundaryValue(const double& x, const double& y, const double& z)
100 {
101 return(getExactSolution('f', x, y, z));
102 }
103
104 static inline char
105 getBoundaryType(const int ID, const double& x, const double& y, const double& z)
106 {
107 return('d');
108 }
109
110 static inline double
111 getExactSolution(const char& what, const double& x,
112 const double& y, const double& z)
113 {
114 if (what == 'f')
115 return(exp(x) + exp(y) + exp(z));
116 else if (what == 'x')
117 return(exp(x));
118 else if (what == 'y')
119 return(exp(y));
120 else if (what == 'z')
121 return(exp(z));
122 else
123 return(0.0);
124 }
125};
126
127// =========== //
128// main driver //
129// =========== //
130
131int main(int argc, char *argv[])
132{
133#ifdef HAVE_MPI
134 MPI_Init(&argc,&argv);
135 Epetra_MpiComm comm(MPI_COMM_WORLD);
136#else
138#endif
139
140 Galeri::core::Workspace::setNumDimensions(3);
141
142 Galeri::grid::Loadable domain, boundary;
143
144 int numGlobalElementsX = 2 * comm.NumProc();
145 int numGlobalElementsY = 2;
146 int numGlobalElementsZ = 2;
147
148 int mx = comm.NumProc();
149 int my = 1;
150 int mz = 1;
151
152 Galeri::grid::Generator::
153 getCubeWithHexs(comm, numGlobalElementsX, numGlobalElementsY, numGlobalElementsZ,
154 mx, my, mz, domain, boundary);
155
156 Epetra_Map matrixMap(domain.getNumGlobalVertices(), 0, comm);
157
158 Epetra_FECrsMatrix A(Copy, matrixMap, 0);
159 Epetra_FEVector LHS(matrixMap);
160 Epetra_FEVector RHS(matrixMap);
161
162 Galeri::problem::ScalarLaplacian<Laplacian> problem("Hex", 1, 8);
163
164 problem.integrate(domain, A, RHS);
165
166 LHS.PutScalar(0.0);
167
168 problem.imposeDirichletBoundaryConditions(boundary, A, RHS, LHS);
169
170 // ============================================================ //
171 // Solving the linear system is the next step, using the IFPACK //
172 // factory. This is done by using the IFPACK factory, then //
173 // asking for IC preconditioner, and setting few parameters //
174 // using a Teuchos::ParameterList. //
175 // ============================================================ //
176
177 Ifpack Factory;
178 Ifpack_Preconditioner* Prec = Factory.Create("IC", &A, 0);
179
180 Teuchos::ParameterList list;
181
182 list.set("fact: level-of-fill", 1);
183 IFPACK_CHK_ERR(Prec->SetParameters(list));
184 IFPACK_CHK_ERR(Prec->Initialize());
185 IFPACK_CHK_ERR(Prec->Compute());
186
187 Epetra_LinearProblem linearProblem(&A, &LHS, &RHS);
188
189 AztecOO solver(linearProblem);
190 solver.SetAztecOption(AZ_solver, AZ_cg);
191 solver.SetPrecOperator(Prec);
192 solver.Iterate(1550, 1e-9);
193
194 // visualization using MEDIT -- a VTK module is available as well
195 Galeri::viz::MEDIT::write(domain, "sol", LHS);
196
197 // now compute the norm of the solution
198 problem.computeNorms(domain, LHS);
199
200#ifdef HAVE_MPI
201 MPI_Finalize();
202#endif
203}
Copy
#define IFPACK_CHK_ERR(ifpack_err)
int main(int argc, char *argv[])
#define RHS(a)
Definition: MatGenFD.c:60
int NumProc() const
int PutScalar(double ScalarConstant)
Ifpack_Preconditioner: basic class for preconditioning in Ifpack.
virtual int Compute()=0
Computes all it is necessary to apply the preconditioner.
virtual int SetParameters(Teuchos::ParameterList &List)=0
Sets all parameters for the preconditioner.
virtual int Initialize()=0
Computes all it is necessary to initialize the preconditioner.
static Ifpack_Preconditioner * Create(EPrecType PrecType, Epetra_RowMatrix *Matrix, const int overlap=0, bool overrideSerialDefault=false)
Creates an instance of Ifpack_Preconditioner given the enum value of the preconditioner type (can not...
Definition: Ifpack.cpp:253
static double getElementRHS(const double &x, const double &y, const double &z, const double &phi_i)
static double getExactSolution(const char &what, const double &x, const double &y, const double &z)
static double getBoundaryValue(const double &x, const double &y, const double &z)
static char getBoundaryType(const int ID, const double &x, const double &y, const double &z)
static double getElementLHS(const double &x, const double &y, const double &z, const double &phi_i, const double &phi_i_derx, const double &phi_i_dery, const double &phi_i_derz, const double &phi_j, const double &phi_j_derx, const double &phi_j_dery, const double &phi_j_derz)