Epetra Package Browser (Single Doxygen Collection) Development
cc_main.cc
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright 2001 Sandia Corporation
00007 // 
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00039 // 
00040 // ************************************************************************
00041 //@HEADER
00042 */
00043 
00044 #include <stdio.h>
00045 #include <stdlib.h>
00046 #include <ctype.h>
00047 #include <assert.h>
00048 #include <string.h>
00049 #include <math.h>
00050 #include "Petra_Comm.h"
00051 #include "Petra_Map.h"
00052 #include "Petra_Time.h"
00053 #include "Petra_RDP_MultiVector.h"
00054 #include "Petra_RDP_Vector.h"
00055 #include "Petra_RDP_CRS_Matrix.h"
00056 #include "Trilinos_LinearProblem.h"
00057 #ifdef PETRA_MPI
00058 #include "mpi.h"
00059 #endif
00060 #ifndef __cplusplus
00061 #define __cplusplus
00062 #endif
00063 
00064 // prototype
00065 #include"basis.h"
00066 
00067 int main(int argc, char *argv[])
00068 {
00069   int ierr = 0, i, j;
00070   bool debug = false;
00071 
00072 #ifdef PETRA_MPI
00073 
00074   // Initialize MPI
00075 
00076   MPI_Init(&argc,&argv);
00077   int size, rank; // Number of MPI processes, My process ID
00078 
00079   MPI_Comm_size(MPI_COMM_WORLD, &size);
00080   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00081 
00082 #else
00083 
00084   int size = 1; // Serial case (not using MPI)
00085   int rank = 0;
00086 
00087 #endif
00088 
00089 
00090 #ifdef PETRA_MPI
00091   Petra_Comm & Comm = *new Petra_Comm( MPI_COMM_WORLD );
00092 #else
00093   Petra_Comm & Comm = *new Petra_Comm();
00094 #endif
00095 
00096   int MyPID = Comm.MyPID();
00097   int NumProc = Comm.NumProc();
00098 
00099   // Get the number of local equations from the command line
00100   if (argc!=2) { 
00101     cout << "Usage: " << argv[0] << " number_of_elements" << endl;
00102     exit(1);
00103   }
00104   int NumGlobalElements = atoi(argv[1])+1;
00105   int IndexBase = 0;
00106 
00107   if (NumGlobalElements < NumProc) {
00108     cout << "numGlobalBlocks = " << NumGlobalElements 
00109    << " cannot be < number of processors = " << NumProc << endl;
00110     exit(1);
00111   }
00112 
00113   // Construct a Source Map that puts approximately the same 
00114   // Number of equations on each processor in uniform global ordering
00115 
00116   Petra_Map& StandardMap = *new Petra_Map(NumGlobalElements, 0, Comm);
00117   int NumMyElements = StandardMap.NumMyElements();
00118   int * StandardMyGlobalElements = new int[NumMyElements];
00119   StandardMap.MyGlobalElements(StandardMyGlobalElements);
00120 
00121   // Create a standard Petra_CRS_Graph
00122   //Petra_CRS_Graph& StandardGraph = *new Petra_CRS_Graph(Copy, StandardMap, 3);
00123   //assert(!StandardGraph.IndicesAreGlobal());
00124   //assert(!StandardGraph.IndicesAreLocal());
00125   
00126   // Construct an Overlapped Map of StandardMap that include 
00127   // the endpoints from two neighboring processors.
00128 
00129   int OverlapNumMyElements;
00130   int OverlapMinMyGID;
00131 
00132   OverlapNumMyElements = NumMyElements + 2;
00133   if ((MyPID==0)||(MyPID==NumProc-1)) OverlapNumMyElements--;
00134 
00135   if (MyPID==0) OverlapMinMyGID = StandardMap.MinMyGID();
00136   else OverlapMinMyGID = StandardMap.MinMyGID()-1;
00137 
00138   int * OverlapMyGlobalElements = new int[OverlapNumMyElements];
00139 
00140   for (i=0; i< OverlapNumMyElements; i++) 
00141                        OverlapMyGlobalElements[i] = OverlapMinMyGID + i;
00142 
00143   Petra_Map& OverlapMap = *new Petra_Map(-1, OverlapNumMyElements, 
00144            OverlapMyGlobalElements, 0, Comm);
00145   
00146   int pS=3;
00147   int pO=3;
00148   // Create Linear Objects for Solve
00149   Petra_RDP_Vector& du = *new Petra_RDP_Vector(StandardMap);
00150   Petra_RDP_Vector& rhs = *new Petra_RDP_Vector(StandardMap);
00151   Petra_RDP_Vector& soln = *new Petra_RDP_Vector(StandardMap);
00152   Petra_RDP_CRS_Matrix& A = *new Petra_RDP_CRS_Matrix(Copy, StandardMap, pS);
00153 
00154   // Create Linear Objects for Fill (Solution Vector)
00155   Petra_RDP_Vector& rhs1 = *new Petra_RDP_Vector(OverlapMap);
00156   Petra_RDP_Vector& u = *new Petra_RDP_Vector(OverlapMap);
00157   Petra_RDP_Vector& x = *new Petra_RDP_Vector(OverlapMap);
00158   Petra_RDP_CRS_Matrix& A1 = *new Petra_RDP_CRS_Matrix(Copy, OverlapMap, pO);
00159 
00160   // Initialize Solution
00161   i=u.PutScalar(1.0);
00162   i=soln.PutScalar(1.0);
00163 
00164   int row;
00165   double eta;
00166   double *xx = new double[2];
00167   double *uu = new double[2];
00168   double jac;
00169   int column;
00170   int *indicies = new int[3];
00171   double *RowValues = new double[OverlapMap.NumGlobalElements()];
00172   Basis basis;
00173   double residual, difference;
00174   double relTol=1.0e-4;
00175   double absTol=1.0e-9;
00176 
00177   // Create the nodal position variables
00178   double Length=1.0;
00179   double dx=Length/((double) NumGlobalElements-1);
00180   for (i=0; i < OverlapNumMyElements; i++) {
00181     x[i]=dx*((double) OverlapMinMyGID+i);
00182   }
00183   
00184   // Begin Nonlinear solver LOOP ************************************
00185   
00186   for (int NLS=0; NLS<2; NLS++) {
00187 
00188 
00189   i=A1.PutScalar(0.0);
00190   i=A.PutScalar(0.0);
00191   i=rhs1.PutScalar(0.0);
00192   i=rhs.PutScalar(0.0);
00193 
00194   // Loop Over # of Finite Elements on Processor
00195   for (int ne=0; ne < OverlapNumMyElements-1; ne++) {
00196     
00197     // Loop Over Gauss Points (5th order GQ)
00198     for(int gp=0; gp < 3; gp++) {
00199       xx[0]=x[ne];
00200       xx[1]=x[ne+1];
00201       uu[0]=u[ne];
00202       uu[1]=u[ne+1];
00203       basis.getBasis(gp, xx, uu);
00204               
00205       // Loop over Nodes in Element
00206       for (i=0; i< 2; i++) {
00207     rhs1[ne+i]+=basis.wt*basis.dx
00208       *((1.0/(basis.dx*basis.dx))*basis.duu*
00209          basis.dphide[i]+basis.uu*basis.uu*basis.phi[i]);
00210     //  printf("Proc=%d, GID=%d, rhs=%e owned%d\n",MyPID ,
00211     //OverlapMap.GID(i),rhs[i],StandardMap.MyGID(OverlapMap.GID(i)));
00212 
00213   // Loop over Trial Functions
00214   for(j=0;j < 2; j++) {
00215     jac=basis.wt*basis.dx*((1.0/(basis.dx*basis.dx))*basis.dphide[j]*
00216          basis.dphide[i]+2.0*basis.uu*basis.phi[j]*basis.phi[i]);
00217     row=OverlapMap.GID(ne+i);
00218     column=OverlapMap.GID(ne+j);
00219     ierr=A1.SumIntoGlobalValues(row, 1, &jac, &column);
00220     if (ierr!=0) {      
00221       // printf("SumInto failed at (%d,%d)!!\n",row,column);
00222       ierr=A1.InsertGlobalValues(row, 1, &jac, &column);
00223       // if (ierr==0) printf("Insert SUCCEEDED at (%d,%d)!!\n",row,column);
00224     } //else if (ierr==0) 
00225       // printf("SumInto SUCCEEDED at (%d,%d)!!\n",row,column);
00226     
00227   }
00228       }
00229     }
00230   }
00231 
00232   Comm.Barrier();
00233 
00234   Petra_Import & Importer = *new Petra_Import(StandardMap, OverlapMap);
00235   assert(rhs.Import(rhs1, Importer, Insert)==0);
00236   assert(A.Import(A1, Importer, Insert)==0);
00237   delete &Importer;
00238 
00239   // Insert Boundary Conditions
00240   // U(0)=1.0
00241   if (MyPID==0) {
00242     u[0]=1.0;
00243     rhs[0]=0.0;
00244     column=0;
00245     jac=1.0;
00246     A.ReplaceGlobalValues(0, 1, &jac, &column);
00247     column=1;
00248     jac=0.0;
00249     A.ReplaceGlobalValues(0, 1, &jac, &column);
00250   }
00251 
00252   Comm.Barrier();
00253 
00254   assert(A.FillComplete()==0);
00255 
00256   /*
00257   // Print Matrix
00258   int StandardNumMyRows = A.NumMyRows();
00259   int * StandardIndices; 
00260   int StandardNumEntries; 
00261   double * StandardValues;
00262   for (i=0; i< StandardNumMyRows; i++) {
00263     A.ExtractMyRowView(i, StandardNumEntries, 
00264             StandardValues, StandardIndices);
00265     for (j=0; j < StandardNumEntries; j++) {
00266       printf("MyPID=%d, J[%d,%d]=%e\n",MyPID,i,j,StandardValues[j]);
00267     }
00268   }
00269   */
00270 
00271   // check if Converged   
00272   ierr=rhs.Norm2(&residual);
00273   ierr=du.Norm2(&difference);
00274   if (MyPID==0) printf("\n***********************************************\n");
00275   if (MyPID==0) printf("Iteration %d  Residual L2=%e   Update L2=%e\n"
00276        ,NLS,residual,difference);
00277   if (MyPID==0) printf("***********************************************\n");  
00278   if ((residual < absTol)&&(difference < relTol)) {
00279     if (MyPID==0) printf("\n\nConvergence Achieved!!!!\n");
00280     return 0;
00281   }    
00282   
00283   Trilinos_LinearProblem *Problem = new  Trilinos_LinearProblem(&A,&du,&rhs);
00284   Problem->SetPDL(hard);
00285   Problem->Iterate(400, 1.0e-8);
00286   delete Problem;
00287 
00288   // Update Solution    
00289   for (i=0;i<NumMyElements;i++) soln[i] -= du[i];
00290    
00291   Petra_Import & Importer2 = *new Petra_Import(OverlapMap, StandardMap);
00292   assert(u.Import(soln, Importer2, Insert)==0);
00293   delete &Importer2;
00294   
00295   for (i=0;i<OverlapNumMyElements;i++) 
00296     printf("Proc=%d GID=%d u=%e soln=%e\n",MyPID,
00297      OverlapMap.GID(i),u[i],soln[i]);
00298 
00299   } // End NLS Loop *****************************************************
00300  
00301 
00302 
00303   delete &OverlapMap;
00304   delete [] OverlapMyGlobalElements;
00305 
00306   //delete &StandardGraph;
00307   delete &StandardMap;
00308   delete [] StandardMyGlobalElements;
00309 
00310   delete &Comm;
00311 
00312 #ifdef PETRA_MPI
00313   MPI_Finalize() ;
00314 #endif
00315 
00316 /* end main
00317 */
00318 return ierr ;
00319 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines