On Fri, Jun 1, 2018 at 12:07 PM, Samuel Lanthaler <s.lanthaler at gmail.com>
wrote:
> On 06/01/2018 03:42 PM, Matthew Knepley wrote:
>
> On Fri, Jun 1, 2018 at 9:21 AM, Samuel Lanthaler <s.lanthaler at gmail.com>
> 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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180601/22c7736f/attachment.html>