[Trilinos-Users] FECrs Vector/Matrix

[Williams, Alan B]

David,
I'm glad to hear you've got it working.

Certainly the BC enforcement via extraction and replacement of rows will
only work on the processor that owns those rows.

Alan
 

> -----Original Message-----
> From: David Neckels [mailto:dneckels at ucar.edu] 
> Sent: Thursday, November 29, 2007 9:52 AM
> To: Williams, Alan B
> Cc: trilinos-users at software.sandia.gov
> Subject: Re: [Trilinos-Users] FECrs Vector/Matrix
> 
> Hi Alan,
> 
> Okay. Using your previous explanation as a guide, I reworked 
> things and figured out what I was doing wrong.
> 
> After creating the matrix and vectors with the unique map, 
> and calling InsertValues, I was calling a version of 
> SumIntoGlobalValues that exists only in the base class, 
> Epetra_CrsMatrix:
> 
> int SumIntoGlobalValues(int GlobalRow, int NumEntries, double* Values,
> int* Indices);
> 
> Obviously, since the function is in the base class it didn't 
> know about the extra (shared rows), and so it would return an error.
> I now call:
> int SumIntoGlobalValues(int numRows, const int* rows, int 
> numCols, const int* cols, const double* values, int 
> format=Epetra_FECrsMatrix::COLUMN_MAJOR);
> 
> from Epetra_FECrsMatrix, which works correctly.
> 
> Also, after assembly I apply dirichlet boundary conditions by 
> replacing the row for a dirichlet dof. I call
> 
> int ret = matrix->ExtractGlobalRowView(row_id, nentries, 
> old_row, old_indices);
> 
> followed by
> 
> matrix->ReplaceGlobalValues(1, &row_id, nentries, old_indices, 
> matrix->&new_row[0]);
> 
> At first this wasn't working correctly for shared dofs, since 
> the ExtractView didn't know about the shared dofs on 
> non-owning side (it is a function in the base class CrsMatrix).
> What seems to work is to call:
> matrix->GlobalAssemble(false /* don't fillcomplete */ );
> and then extract the row and replace, on the owning processor.....
> 
> I don't know how I was getting away with my previous 
> unique/shared map approach. It was giving the right answer, 
> but I notice that my matrices are much better conditioned 
> now.... A fluke, I guess.
> 
> Thanks a lot for the help.
> 
> -David N.
> 
> Williams, Alan B wrote:
> > If you can give me some example code that reproduces the -1 
> > return-errors you're seeing, I'll figure out what's going on.
> >
> > Perhaps there is a bug in FECrsMatrix.
> >
> > Alan
> >
> >
> >   
> >> -----Original Message-----
> >> From: David Neckels [mailto:dneckels at ucar.edu]
> >> Sent: Wednesday, November 28, 2007 6:07 PM
> >> To: Williams, Alan B
> >> Cc: trilinos-users at software.sandia.gov
> >> Subject: Re: [Trilinos-Users] FECrs Vector/Matrix
> >>
> >> Hi Alan,
> >>
> >> Thanks for responding, good to hear from you, hope all is well!
> >>
> >> I appreciate the basic explanation, I think it helps 
> understand the 
> >> theory on using the classes better.
> >> However, I still don't seem to have luck with things. I have now 
> >> tried:
> >>
> >> matrix = new Epetra_FECrsMatrix(Copy, *unique_map, max_df);
> >>
> >> x = new Epetra_FEVector(*unique_map); //( solving Ax=b) b = new 
> >> Epetra_FEVector(*unique_map);
> >>
> >> Call InsertGlobalValues for the shared set of dofs to build the 
> >> matrix stencil.
> >>
> >> Then,
> >> during assembly call SumIntoGlobalValues for the same set of dofs 
> >> (the shared and unique dofs).
> >> However, for shared dofs this function now returns -1 !!
> >>
> >> I then close the assembly:
> >> matrix->GlobalAssembly();
> >> b->GlobalAssembly()
> >>
> >> This gives the wrong result...
> >> I notice that calling the InsertGlobalValues alone does appear to 
> >> correctly sum the entries, but it doesn't seem to play well with 
> >> SumIntoGlobalValues for some reason.
> >>
> >> For this example I am not using ReplaceGlobalValues; I have no 
> >> constraints...
> >>
> >> If you like I can try to boil this down to a simple example....??
> >>
> >> Thanks,
> >> David N.
> >>
> >>
> >> // insert rows with all (shared) dofs // sum into global with all 
> >> (shared) dofs --here is I used unique_map to create the matrix, 
> >> above, I would get bad return values??
> >> // Call ExtractGlobalView and replace constrained rows of 
> the matrix 
> >> for shared dofs (unique_map in the matrix create caused 
> this to fail) 
> >> // Also sum into shared dofs of b
> >>
> >> // After assembly:
> >> matrix->GlobalAssemble(*umap, *umap);
> >> b->GlobalAssemble();
> >>
> >>
> >>
> >>
> >> Williams, Alan B wrote:
> >>     
> >>> Hi David,
> >>>
> >>> The "FE" epetra objects are intended to handle the situation you 
> >>> describe, where shared nodes exist on multiple processors,
> >>>       
> >> so you're
> >>     
> >>> on the right track.
> >>>
> >>> As you know, the non-FE epetra objects tend to be used for the 
> >>> uniquely-owned case, where the Epetra_Map attribute has 
> >>> non-overlapping entries.
> >>>
> >>> The general idea is that the FE objects allow a
> >>>       
> >> finite-element caller
> >>     
> >>> to assemble data on the local processor (i.e. overlapping 
> data for 
> >>> shared nodes), which may end up on another processor in the 
> >>> non-overlapping distribution.
> >>>
> >>> The FE objects inherit their non-FE counterparts
> >>>       
> >> (FECrsMatrix inherits
> >>     
> >>> CrsMatrix, etc), and the Epetra_Map objects that are passed as 
> >>> constructor arguments should describe the uniquely-owned, 
> >>> non-overlapping distribution (which is what you want the 
> underlying 
> >>> CrsMatrix to end up with).
> >>>
> >>> When you assemble data into an FE object, it checks whether
> >>>       
> >> the data
> >>     
> >>> is local according to the Epetra_Map object it was
> >>>       
> >> constructed with. 
> >>     
> >>> If so, that data is passed straight through to the non-FE
> >>>       
> >> base-class
> >>     
> >>> object. If the data belongs on another processor, then it is held 
> >>> locally in "temporary" data structures until the
> >>>       
> >> GlobalAssemble method is called.
> >>     
> >>> GlobalAssemble uses Import/Export operations to send that
> >>>       
> >> data to the
> >>     
> >>> appropriate owning processors, where it is placed into the non-FE 
> >>> base-class object.
> >>>
> >>> So you should create your FECrsMatrix with the
> >>>       
> >> non-overlapping row-map.
> >>     
> >>> FECrsMatrix shouldn't need your shared-map at all. Note
> >>>       
> >> that the map
> >>     
> >>> arguments which are passed to the GlobalAssemble method are the 
> >>> domain-map and range-map for the underlying CrsMatrix, 
> which is the 
> >>> target of the global assembly. In general, a CrsMatrix has 4 maps:
> >>> row-map, column-map, domain-map and range-map. In some cases the 
> >>> domain-map is different than the column-map, and the range-map is 
> >>> different than the row-map. See the doxygen documentation for 
> >>> Epetra_CrsGraph for a description of these 4 map attributes. (At 
> >>> trilinos.sandia.gov, go to packages => epetra => documentation.)
> >>>
> >>> The GlobalAssemble method should work without any map 
> arguments in 
> >>> most cases. (It internally obtains the domain-map and
> >>>       
> >> range-map from
> >>     
> >>> the underlying CrsMatrix.)
> >>>
> >>> Let me know if this doesn't clear things up for you, or 
> if you are 
> >>> still seeing errors.
> >>>
> >>> Note that non-zero return-values for SumIntoGlobalValues don't 
> >>> necessarily indicate an error. If the return-value is a positive 
> >>> number, it simply means that the matrix did an internal
> >>>       
> >> allocation in
> >>     
> >>> the process of storing the input values, because it wasn't
> >>>       
> >> constructed
> >>     
> >>> with enough information to pre-allocate all needed nonzero
> >>>       
> >> positions.
> >>     
> >>> If you see a negative return-value, then that's an error.
> >>>
> >>> Alan
> >>>
> >>>
> >>>   
> >>>       
> >>>> -----Original Message-----
> >>>> From: trilinos-users-bounces at software.sandia.gov
> >>>> [mailto:trilinos-users-bounces at software.sandia.gov] On Behalf Of 
> >>>> David Neckels
> >>>> Sent: Wednesday, November 28, 2007 3:12 PM
> >>>> To: trilinos-users at software.sandia.gov
> >>>> Subject: [Trilinos-Users] FECrs Vector/Matrix
> >>>>
> >>>> Hi,
> >>>>
> >>>> I installed the Epetra_FECrsMatrix and Epetra_FEVector
> >>>>         
> >> classes in a
> >>     
> >>>> parallel finite element code, but had some trouble figured
> >>>>         
> >> out how to
> >>     
> >>>> use the classes and their Epetra_Map's.  I settled on
> >>>>         
> >> something that
> >>     
> >>>> seems to work, but wasn't sure if it was correct?
> >>>>
> >>>> My code partitions on element boundary's, so faces, edges,
> >>>>         
> >> and nodes
> >>     
> >>>> on a shared element side exist on more than one processor.
> >>>> Hence although each degree of freedom is owned by only one
> >>>>         
> >> processor,
> >>     
> >>>> it may reside on multiple processors.
> >>>>
> >>>> I have two corresponding Epetra_Map's:
> >>>> unique_map: contains only the owned dofs on a processor
> >>>>         
> >> (partitions
> >>     
> >>>> the problem by processors).
> >>>> shared_map: has overlap, i.e. a dof will appear on each 
> processor 
> >>>> where the dof resides.
> >>>>
> >>>> This was the sequence that seemed to work:
> >>>>
> >>>> matrix = new Epetra_FECrsMatrix(Copy, *shared_map,
> >>>>         
> >> max_df); x = new
> >>     
> >>>> Epetra_FEVector(*shared_map); //( solving Ax=b)  b = new 
> >>>> Epetra_FEVector(*unique_map);
> >>>>
> >>>> // insert rows with all (shared) dofs // sum into global with all
> >>>> (shared) dofs --here is I used unique_map to create the matrix, 
> >>>> above, I would get bad return values??
> >>>> // Call ExtractGlobalView and replace constrained rows of
> >>>>         
> >> the matrix
> >>     
> >>>> for shared dofs (unique_map in the matrix create caused
> >>>>         
> >> this to fail)
> >>     
> >>>> // Also sum into shared dofs of b
> >>>>
> >>>> // After assembly:
> >>>> matrix->GlobalAssemble(*umap, *umap);
> >>>> b->GlobalAssemble();
> >>>>
> >>>> // Call Aztec solver....
> >>>>
> >>>> This dance between the unique and shared maps seemed strange. 
> >>>>  I just wondered if a more thorough example for using the
> >>>>         
> >> FE classes
> >>     
> >>>> existed, or if someone knew what the correct approach was? 
> >>>>         
> >>  For now
> >>     
> >>>> things seem to work....
> >>>>
> >>>> I have attached the routine where I set these maps up.
> >>>>
> >>>> Thanks for any help.
> >>>>
> >>>> -David N.
> >>>>
> >>>>
> >>>>
> >>>>     
> >>>>         
> >>>   
> >>>       
> >>
> >>     
> >
> >   
> 
> 
>