cc_main.cc

Go to the documentation of this file.
00001 //@HEADER
00002 /*
00003 ************************************************************************
00004 
00005               Epetra: Linear Algebra Services Package 
00006                 Copyright (2001) Sandia Corporation
00007 
00008 Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 license for use of this work by or on behalf of the U.S. Government.
00010 
00011 This library is free software; you can redistribute it and/or modify
00012 it under the terms of the GNU Lesser General Public License as
00013 published by the Free Software Foundation; either version 2.1 of the
00014 License, or (at your option) any later version.
00015  
00016 This library is distributed in the hope that it will be useful, but
00017 WITHOUT ANY WARRANTY; without even the implied warranty of
00018 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 Lesser General Public License for more details.
00020  
00021 You should have received a copy of the GNU Lesser General Public
00022 License along with this library; if not, write to the Free Software
00023 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 USA
00025 Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00026 
00027 ************************************************************************
00028 */
00029 //@HEADER
00030 
00031 #include <stdio.h>
00032 #include <stdlib.h>
00033 #include <ctype.h>
00034 #include <assert.h>
00035 #include <string.h>
00036 #include <math.h>
00037 #include "Petra_Comm.h"
00038 #include "Petra_Map.h"
00039 #include "Petra_Time.h"
00040 #include "Petra_RDP_MultiVector.h"
00041 #include "Petra_RDP_Vector.h"
00042 #include "Petra_RDP_CRS_Matrix.h"
00043 #include "Trilinos_LinearProblem.h"
00044 #ifdef PETRA_MPI
00045 #include "mpi.h"
00046 #endif
00047 #ifndef __cplusplus
00048 #define __cplusplus
00049 #endif
00050 
00051 // prototype
00052 #include"basis.h"
00053 
00054 int main(int argc, char *argv[])
00055 {
00056   int ierr = 0, i, j;
00057   bool debug = false;
00058 
00059 #ifdef PETRA_MPI
00060 
00061   // Initialize MPI
00062 
00063   MPI_Init(&argc,&argv);
00064   int size, rank; // Number of MPI processes, My process ID
00065 
00066   MPI_Comm_size(MPI_COMM_WORLD, &size);
00067   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00068 
00069 #else
00070 
00071   int size = 1; // Serial case (not using MPI)
00072   int rank = 0;
00073 
00074 #endif
00075 
00076 
00077 #ifdef PETRA_MPI
00078   Petra_Comm & Comm = *new Petra_Comm( MPI_COMM_WORLD );
00079 #else
00080   Petra_Comm & Comm = *new Petra_Comm();
00081 #endif
00082 
00083   int MyPID = Comm.MyPID();
00084   int NumProc = Comm.NumProc();
00085 
00086   // Get the number of local equations from the command line
00087   if (argc!=2) { 
00088     cout << "Usage: " << argv[0] << " number_of_elements" << endl;
00089     exit(1);
00090   }
00091   int NumGlobalElements = atoi(argv[1])+1;
00092   int IndexBase = 0;
00093 
00094   if (NumGlobalElements < NumProc) {
00095     cout << "numGlobalBlocks = " << NumGlobalElements 
00096    << " cannot be < number of processors = " << NumProc << endl;
00097     exit(1);
00098   }
00099 
00100   // Construct a Source Map that puts approximately the same 
00101   // Number of equations on each processor in uniform global ordering
00102 
00103   Petra_Map& StandardMap = *new Petra_Map(NumGlobalElements, 0, Comm);
00104   int NumMyElements = StandardMap.NumMyElements();
00105   int * StandardMyGlobalElements = new int[NumMyElements];
00106   StandardMap.MyGlobalElements(StandardMyGlobalElements);
00107 
00108   // Create a standard Petra_CRS_Graph
00109   //Petra_CRS_Graph& StandardGraph = *new Petra_CRS_Graph(Copy, StandardMap, 3);
00110   //assert(!StandardGraph.IndicesAreGlobal());
00111   //assert(!StandardGraph.IndicesAreLocal());
00112   
00113   // Construct an Overlapped Map of StandardMap that include 
00114   // the endpoints from two neighboring processors.
00115 
00116   int OverlapNumMyElements;
00117   int OverlapMinMyGID;
00118 
00119   OverlapNumMyElements = NumMyElements + 2;
00120   if ((MyPID==0)||(MyPID==NumProc-1)) OverlapNumMyElements--;
00121 
00122   if (MyPID==0) OverlapMinMyGID = StandardMap.MinMyGID();
00123   else OverlapMinMyGID = StandardMap.MinMyGID()-1;
00124 
00125   int * OverlapMyGlobalElements = new int[OverlapNumMyElements];
00126 
00127   for (i=0; i< OverlapNumMyElements; i++) 
00128                        OverlapMyGlobalElements[i] = OverlapMinMyGID + i;
00129 
00130   Petra_Map& OverlapMap = *new Petra_Map(-1, OverlapNumMyElements, 
00131            OverlapMyGlobalElements, 0, Comm);
00132   
00133   int pS=3;
00134   int pO=3;
00135   // Create Linear Objects for Solve
00136   Petra_RDP_Vector& du = *new Petra_RDP_Vector(StandardMap);
00137   Petra_RDP_Vector& rhs = *new Petra_RDP_Vector(StandardMap);
00138   Petra_RDP_Vector& soln = *new Petra_RDP_Vector(StandardMap);
00139   Petra_RDP_CRS_Matrix& A = *new Petra_RDP_CRS_Matrix(Copy, StandardMap, pS);
00140 
00141   // Create Linear Objects for Fill (Solution Vector)
00142   Petra_RDP_Vector& rhs1 = *new Petra_RDP_Vector(OverlapMap);
00143   Petra_RDP_Vector& u = *new Petra_RDP_Vector(OverlapMap);
00144   Petra_RDP_Vector& x = *new Petra_RDP_Vector(OverlapMap);
00145   Petra_RDP_CRS_Matrix& A1 = *new Petra_RDP_CRS_Matrix(Copy, OverlapMap, pO);
00146 
00147   // Initialize Solution
00148   i=u.PutScalar(1.0);
00149   i=soln.PutScalar(1.0);
00150 
00151   int row;
00152   double eta;
00153   double *xx = new double[2];
00154   double *uu = new double[2];
00155   double jac;
00156   int column;
00157   int *indicies = new int[3];
00158   double *RowValues = new double[OverlapMap.NumGlobalElements()];
00159   Basis basis;
00160   double residual, difference;
00161   double relTol=1.0e-4;
00162   double absTol=1.0e-9;
00163 
00164   // Create the nodal position variables
00165   double Length=1.0;
00166   double dx=Length/((double) NumGlobalElements-1);
00167   for (i=0; i < OverlapNumMyElements; i++) {
00168     x[i]=dx*((double) OverlapMinMyGID+i);
00169   }
00170   
00171   // Begin Nonlinear solver LOOP ************************************
00172   
00173   for (int NLS=0; NLS<2; NLS++) {
00174 
00175 
00176   i=A1.PutScalar(0.0);
00177   i=A.PutScalar(0.0);
00178   i=rhs1.PutScalar(0.0);
00179   i=rhs.PutScalar(0.0);
00180 
00181   // Loop Over # of Finite Elements on Processor
00182   for (int ne=0; ne < OverlapNumMyElements-1; ne++) {
00183     
00184     // Loop Over Gauss Points (5th order GQ)
00185     for(int gp=0; gp < 3; gp++) {
00186       xx[0]=x[ne];
00187       xx[1]=x[ne+1];
00188       uu[0]=u[ne];
00189       uu[1]=u[ne+1];
00190       basis.getBasis(gp, xx, uu);
00191               
00192       // Loop over Nodes in Element
00193       for (i=0; i< 2; i++) {
00194     rhs1[ne+i]+=basis.wt*basis.dx
00195       *((1.0/(basis.dx*basis.dx))*basis.duu*
00196          basis.dphide[i]+basis.uu*basis.uu*basis.phi[i]);
00197     //  printf("Proc=%d, GID=%d, rhs=%e owned%d\n",MyPID ,
00198     //OverlapMap.GID(i),rhs[i],StandardMap.MyGID(OverlapMap.GID(i)));
00199 
00200   // Loop over Trial Functions
00201   for(j=0;j < 2; j++) {
00202     jac=basis.wt*basis.dx*((1.0/(basis.dx*basis.dx))*basis.dphide[j]*
00203          basis.dphide[i]+2.0*basis.uu*basis.phi[j]*basis.phi[i]);
00204     row=OverlapMap.GID(ne+i);
00205     column=OverlapMap.GID(ne+j);
00206     ierr=A1.SumIntoGlobalValues(row, 1, &jac, &column);
00207     if (ierr!=0) {      
00208       // printf("SumInto failed at (%d,%d)!!\n",row,column);
00209       ierr=A1.InsertGlobalValues(row, 1, &jac, &column);
00210       // if (ierr==0) printf("Insert SUCCEEDED at (%d,%d)!!\n",row,column);
00211     } //else if (ierr==0) 
00212       // printf("SumInto SUCCEEDED at (%d,%d)!!\n",row,column);
00213     
00214   }
00215       }
00216     }
00217   }
00218 
00219   Comm.Barrier();
00220 
00221   Petra_Import & Importer = *new Petra_Import(StandardMap, OverlapMap);
00222   assert(rhs.Import(rhs1, Importer, Insert)==0);
00223   assert(A.Import(A1, Importer, Insert)==0);
00224   delete &Importer;
00225 
00226   // Insert Boundary Conditions
00227   // U(0)=1.0
00228   if (MyPID==0) {
00229     u[0]=1.0;
00230     rhs[0]=0.0;
00231     column=0;
00232     jac=1.0;
00233     A.ReplaceGlobalValues(0, 1, &jac, &column);
00234     column=1;
00235     jac=0.0;
00236     A.ReplaceGlobalValues(0, 1, &jac, &column);
00237   }
00238 
00239   Comm.Barrier();
00240 
00241   assert(A.FillComplete()==0);
00242 
00243   /*
00244   // Print Matrix
00245   int StandardNumMyRows = A.NumMyRows();
00246   int * StandardIndices; 
00247   int StandardNumEntries; 
00248   double * StandardValues;
00249   for (i=0; i< StandardNumMyRows; i++) {
00250     A.ExtractMyRowView(i, StandardNumEntries, 
00251             StandardValues, StandardIndices);
00252     for (j=0; j < StandardNumEntries; j++) {
00253       printf("MyPID=%d, J[%d,%d]=%e\n",MyPID,i,j,StandardValues[j]);
00254     }
00255   }
00256   */
00257 
00258   // check if Converged   
00259   ierr=rhs.Norm2(&residual);
00260   ierr=du.Norm2(&difference);
00261   if (MyPID==0) printf("\n***********************************************\n");
00262   if (MyPID==0) printf("Iteration %d  Residual L2=%e   Update L2=%e\n"
00263        ,NLS,residual,difference);
00264   if (MyPID==0) printf("***********************************************\n");  
00265   if ((residual < absTol)&&(difference < relTol)) {
00266     if (MyPID==0) printf("\n\nConvergence Achieved!!!!\n");
00267     return 0;
00268   }    
00269   
00270   Trilinos_LinearProblem *Problem = new  Trilinos_LinearProblem(&A,&du,&rhs);
00271   Problem->SetPDL(hard);
00272   Problem->Iterate(400, 1.0e-8);
00273   delete Problem;
00274 
00275   // Update Solution    
00276   for (i=0;i<NumMyElements;i++) soln[i] -= du[i];
00277    
00278   Petra_Import & Importer2 = *new Petra_Import(OverlapMap, StandardMap);
00279   assert(u.Import(soln, Importer2, Insert)==0);
00280   delete &Importer2;
00281   
00282   for (i=0;i<OverlapNumMyElements;i++) 
00283     printf("Proc=%d GID=%d u=%e soln=%e\n",MyPID,
00284      OverlapMap.GID(i),u[i],soln[i]);
00285 
00286   } // End NLS Loop *****************************************************
00287  
00288 
00289 
00290   delete &OverlapMap;
00291   delete [] OverlapMyGlobalElements;
00292 
00293   //delete &StandardGraph;
00294   delete &StandardMap;
00295   delete [] StandardMyGlobalElements;
00296 
00297   delete &Comm;
00298 
00299 #ifdef PETRA_MPI
00300   MPI_Finalize() ;
00301 #endif
00302 
00303 /* end main
00304 */
00305 return ierr ;
00306 }

Generated on Thu Sep 18 12:37:56 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1