On Fri, Jun 1, 2018 at 12:07 PM, Samuel Lanthaler <[email protected]> wrote:
> On 06/01/2018 03:42 PM, Matthew Knepley wrote: > > On Fri, Jun 1, 2018 at 9:21 AM, Samuel Lanthaler <[email protected]> > wrote: > >> Hi, >> >> I was wondering what the most efficient way to use MatPtAP would be in >> the following situation: I am discretizing a PDE system. The discretization >> yields a matrix A that has a band structure (with k upper and lower bands, >> say). In order to implement the boundary conditions, I use a transformation >> matrix P which is essentially the unit matrix, except for the entries >> P_{ij} where i,j<k and n-i,n-j<k, so >> >> P = [ B, 0, 0, 0, ..., 0, 0 ] >> [ 0, 1, 0, 0, ..., 0, 0 ] >> [ ] >> [ ] >> [ ..., 1, 0 ] >> [ 0, 0, 0, 0, ..., 0, C ] >> >> with B,C are (k-by-k) matrices. >> Right now, I'm simply constructing A, P and calling >> >> CALL MatPtAP(petsc_matA,petsc_matP,MAT_INITIAL_MATRIX,PETSC_DEFAU >> LT_REAL,petsc_matPtAP,ierr) >> >> where I haven't done anything to pestc_matPtAP, prior to this call. Is >> this the way to do it? >> >> I'm asking because, currently, setting up the matrices A and P takes very >> little time, whereas the operation MatPtAP is taking quite long, which >> seems very odd... The matrices are of type MPIAIJ. In my problem, the total >> matrix dimension is around 10'000 and the matrix blocks (B,C) are of size >> ~100. >> > > Are you sure this is what you want to do? Usually BC are local, since by > definition PDE are local, and > are applied pointwise. What kind of BC do you have here? > > > The boundary conditions are a mixture of Dirichlet and Neumann; in my > case, the PDE is a system involving 8 variables on a disk, where the > periodic direction is discretized using a Fourier series expansion, the > radial direction uses B-splines. > > In reality, I have two matrices A,B, and want to solve the eigenvalue > problem \lambda*B*x = A*x. > I found it quite convenient to use a transformation P to a different set > of variables y, such that x=P*y and x satisfies the BC iff certain > components of y are 0. The latter is enforced by inserting spurious > eigenvalues at the relevant components of y in the transformed eigenvalue > problem \lambda*Pt*B*P*y=Pt*A*P*y. After solving the EVP in terms of y, I > get back x=P*y. > Is this an inherently bad/inefficient way of enforcing BC's? Thanks. > I cannot quite see what is happening. However it appears that you are coupling all the unknowns on a boundary to enforce the boundary condition. If this is true, its pretty expensive. For example, in 2D you have N^2 unknowns for N boundary unknowns. However, dense coupling on the boundary produces N^2 storage and N^3 computation, which is not great. Am I missing something? Its fine to do this for some sizes (2D BEM guys do it all the time), but for bigger stuff it gets untenable. Thanks, Matt > > > Thanks, > > Matt > > >> Thanks in advance for any ideas. >> >> Cheers, >> Samuel >> > > > > -- > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > -- Norbert Wiener > > https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/%7Emk51/> > > > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
