Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
cc_main.cc
Go to the documentation of this file.
1/*
2//@HEADER
3// ************************************************************************
4//
5// Epetra: Linear Algebra Services Package
6// Copyright 2011 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ************************************************************************
41//@HEADER
42*/
43
44#include <stdio.h>
45#include <stdlib.h>
46#include <ctype.h>
47#include <assert.h>
48#include <string.h>
49#include <math.h>
50#include "Petra_Comm.h"
51#include "Petra_Map.h"
52#include "Petra_Time.h"
53#include "Petra_RDP_MultiVector.h"
54#include "Petra_RDP_Vector.h"
55#include "Petra_RDP_CRS_Matrix.h"
56#include "Trilinos_LinearProblem.h"
57#ifdef PETRA_MPI
58#include "mpi.h"
59#endif
60#ifndef __cplusplus
61#define __cplusplus
62#endif
63
64// prototype
65#include"basis.h"
66
67int main(int argc, char *argv[])
68{
69 int ierr = 0, i, j;
70 bool debug = false;
71
72#ifdef PETRA_MPI
73
74 // Initialize MPI
75
76 MPI_Init(&argc,&argv);
77 int size, rank; // Number of MPI processes, My process ID
78
79 MPI_Comm_size(MPI_COMM_WORLD, &size);
80 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
81
82#else
83
84 int size = 1; // Serial case (not using MPI)
85 int rank = 0;
86
87#endif
88
89
90#ifdef PETRA_MPI
91 Petra_Comm & Comm = *new Petra_Comm( MPI_COMM_WORLD );
92#else
93 Petra_Comm & Comm = *new Petra_Comm();
94#endif
95
96 int MyPID = Comm.MyPID();
97 int NumProc = Comm.NumProc();
98
99 // Get the number of local equations from the command line
100 if (argc!=2) {
101 cout << "Usage: " << argv[0] << " number_of_elements" << endl;
102 exit(1);
103 }
104 int NumGlobalElements = atoi(argv[1])+1;
105 int IndexBase = 0;
106
107 if (NumGlobalElements < NumProc) {
108 cout << "numGlobalBlocks = " << NumGlobalElements
109 << " cannot be < number of processors = " << NumProc << endl;
110 exit(1);
111 }
112
113 // Construct a Source Map that puts approximately the same
114 // Number of equations on each processor in uniform global ordering
115
116 Petra_Map& StandardMap = *new Petra_Map(NumGlobalElements, 0, Comm);
117 int NumMyElements = StandardMap.NumMyElements();
118 int * StandardMyGlobalElements = new int[NumMyElements];
119 StandardMap.MyGlobalElements(StandardMyGlobalElements);
120
121 // Create a standard Petra_CRS_Graph
122 //Petra_CRS_Graph& StandardGraph = *new Petra_CRS_Graph(Copy, StandardMap, 3);
123 //assert(!StandardGraph.IndicesAreGlobal());
124 //assert(!StandardGraph.IndicesAreLocal());
125
126 // Construct an Overlapped Map of StandardMap that include
127 // the endpoints from two neighboring processors.
128
129 int OverlapNumMyElements;
130 int OverlapMinMyGID;
131
132 OverlapNumMyElements = NumMyElements + 2;
133 if ((MyPID==0)||(MyPID==NumProc-1)) OverlapNumMyElements--;
134
135 if (MyPID==0) OverlapMinMyGID = StandardMap.MinMyGID();
136 else OverlapMinMyGID = StandardMap.MinMyGID()-1;
137
138 int * OverlapMyGlobalElements = new int[OverlapNumMyElements];
139
140 for (i=0; i< OverlapNumMyElements; i++)
141 OverlapMyGlobalElements[i] = OverlapMinMyGID + i;
142
143 Petra_Map& OverlapMap = *new Petra_Map(-1, OverlapNumMyElements,
144 OverlapMyGlobalElements, 0, Comm);
145
146 int pS=3;
147 int pO=3;
148 // Create Linear Objects for Solve
149 Petra_RDP_Vector& du = *new Petra_RDP_Vector(StandardMap);
150 Petra_RDP_Vector& rhs = *new Petra_RDP_Vector(StandardMap);
151 Petra_RDP_Vector& soln = *new Petra_RDP_Vector(StandardMap);
152 Petra_RDP_CRS_Matrix& A = *new Petra_RDP_CRS_Matrix(Copy, StandardMap, pS);
153
154 // Create Linear Objects for Fill (Solution Vector)
155 Petra_RDP_Vector& rhs1 = *new Petra_RDP_Vector(OverlapMap);
156 Petra_RDP_Vector& u = *new Petra_RDP_Vector(OverlapMap);
157 Petra_RDP_Vector& x = *new Petra_RDP_Vector(OverlapMap);
158 Petra_RDP_CRS_Matrix& A1 = *new Petra_RDP_CRS_Matrix(Copy, OverlapMap, pO);
159
160 // Initialize Solution
161 i=u.PutScalar(1.0);
162 i=soln.PutScalar(1.0);
163
164 int row;
165 double eta;
166 double *xx = new double[2];
167 double *uu = new double[2];
168 double jac;
169 int column;
170 int *indicies = new int[3];
171 double *RowValues = new double[OverlapMap.NumGlobalElements()];
172 Basis basis;
173 double residual, difference;
174 double relTol=1.0e-4;
175 double absTol=1.0e-9;
176
177 // Create the nodal position variables
178 double Length=1.0;
179 double dx=Length/((double) NumGlobalElements-1);
180 for (i=0; i < OverlapNumMyElements; i++) {
181 x[i]=dx*((double) OverlapMinMyGID+i);
182 }
183
184 // Begin Nonlinear solver LOOP ************************************
185
186 for (int NLS=0; NLS<2; NLS++) {
187
188
189 i=A1.PutScalar(0.0);
190 i=A.PutScalar(0.0);
191 i=rhs1.PutScalar(0.0);
192 i=rhs.PutScalar(0.0);
193
194 // Loop Over # of Finite Elements on Processor
195 for (int ne=0; ne < OverlapNumMyElements-1; ne++) {
196
197 // Loop Over Gauss Points (5th order GQ)
198 for(int gp=0; gp < 3; gp++) {
199 xx[0]=x[ne];
200 xx[1]=x[ne+1];
201 uu[0]=u[ne];
202 uu[1]=u[ne+1];
203 basis.getBasis(gp, xx, uu);
204
205 // Loop over Nodes in Element
206 for (i=0; i< 2; i++) {
207 rhs1[ne+i]+=basis.wt*basis.dx
208 *((1.0/(basis.dx*basis.dx))*basis.duu*
209 basis.dphide[i]+basis.uu*basis.uu*basis.phi[i]);
210 // printf("Proc=%d, GID=%d, rhs=%e owned%d\n",MyPID ,
211 //OverlapMap.GID(i),rhs[i],StandardMap.MyGID(OverlapMap.GID(i)));
212
213 // Loop over Trial Functions
214 for(j=0;j < 2; j++) {
215 jac=basis.wt*basis.dx*((1.0/(basis.dx*basis.dx))*basis.dphide[j]*
216 basis.dphide[i]+2.0*basis.uu*basis.phi[j]*basis.phi[i]);
217 row=OverlapMap.GID(ne+i);
218 column=OverlapMap.GID(ne+j);
219 ierr=A1.SumIntoGlobalValues(row, 1, &jac, &column);
220 if (ierr!=0) {
221 // printf("SumInto failed at (%d,%d)!!\n",row,column);
222 ierr=A1.InsertGlobalValues(row, 1, &jac, &column);
223 // if (ierr==0) printf("Insert SUCCEEDED at (%d,%d)!!\n",row,column);
224 } //else if (ierr==0)
225 // printf("SumInto SUCCEEDED at (%d,%d)!!\n",row,column);
226
227 }
228 }
229 }
230 }
231
232 Comm.Barrier();
233
234 Petra_Import & Importer = *new Petra_Import(StandardMap, OverlapMap);
235 assert(rhs.Import(rhs1, Importer, Insert)==0);
236 assert(A.Import(A1, Importer, Insert)==0);
237 delete &Importer;
238
239 // Insert Boundary Conditions
240 // U(0)=1.0
241 if (MyPID==0) {
242 u[0]=1.0;
243 rhs[0]=0.0;
244 column=0;
245 jac=1.0;
246 A.ReplaceGlobalValues(0, 1, &jac, &column);
247 column=1;
248 jac=0.0;
249 A.ReplaceGlobalValues(0, 1, &jac, &column);
250 }
251
252 Comm.Barrier();
253
254 assert(A.FillComplete()==0);
255
256 /*
257 // Print Matrix
258 int StandardNumMyRows = A.NumMyRows();
259 int * StandardIndices;
260 int StandardNumEntries;
261 double * StandardValues;
262 for (i=0; i< StandardNumMyRows; i++) {
263 A.ExtractMyRowView(i, StandardNumEntries,
264 StandardValues, StandardIndices);
265 for (j=0; j < StandardNumEntries; j++) {
266 printf("MyPID=%d, J[%d,%d]=%e\n",MyPID,i,j,StandardValues[j]);
267 }
268 }
269 */
270
271 // check if Converged
272 ierr=rhs.Norm2(&residual);
273 ierr=du.Norm2(&difference);
274 if (MyPID==0) printf("\n***********************************************\n");
275 if (MyPID==0) printf("Iteration %d Residual L2=%e Update L2=%e\n"
276 ,NLS,residual,difference);
277 if (MyPID==0) printf("***********************************************\n");
278 if ((residual < absTol)&&(difference < relTol)) {
279 if (MyPID==0) printf("\n\nConvergence Achieved!!!!\n");
280 return 0;
281 }
282
283 Trilinos_LinearProblem *Problem = new Trilinos_LinearProblem(&A,&du,&rhs);
284 Problem->SetPDL(hard);
285 Problem->Iterate(400, 1.0e-8);
286 delete Problem;
287
288 // Update Solution
289 for (i=0;i<NumMyElements;i++) soln[i] -= du[i];
290
291 Petra_Import & Importer2 = *new Petra_Import(OverlapMap, StandardMap);
292 assert(u.Import(soln, Importer2, Insert)==0);
293 delete &Importer2;
294
295 for (i=0;i<OverlapNumMyElements;i++)
296 printf("Proc=%d GID=%d u=%e soln=%e\n",MyPID,
297 OverlapMap.GID(i),u[i],soln[i]);
298
299 } // End NLS Loop *****************************************************
300
301
302
303 delete &OverlapMap;
304 delete [] OverlapMyGlobalElements;
305
306 //delete &StandardGraph;
307 delete &StandardMap;
308 delete [] StandardMyGlobalElements;
309
310 delete &Comm;
311
312#ifdef PETRA_MPI
313 MPI_Finalize() ;
314#endif
315
316/* end main
317*/
318return ierr ;
319}
int main(int argc, char *argv[])
Definition cc_main.cc:67
Definition basis.h:50
double wt
Definition basis.h:54
double * phi
Definition basis.h:53
void getBasis(int gp, double *x, double *u)
Definition basis.cc:61
double duu
Definition basis.h:54
double uu
Definition basis.h:54
double dx
Definition basis.h:58
double * dphide
Definition basis.h:53