[petsc-users] Block unstructured grid

[Matthew Knepley]

On Thu, Mar 11, 2021 at 8:49 AM Mathieu Dutour <mathieu.dutour at gmail.com>
wrote:

> On Thu, 11 Mar 2021 at 13:52, Mark Adams <mfadams at lbl.gov> wrote:
>
>> Mathieu,
>> We have "FieldSplit" support for fields, but I don't know if it has ever
>> been pushed to 1000's of fields so it might fall down. It might work.
>> FieldSplit lets you manipulate the ordering, say field major (j) or node
>> major (i).
>>
> I just looked at it and FieldSplit appears to be used in preconditioner so
> not exactly relevant.
>
> What was unsatisfactory?
>> It sounds like you made a rectangular matrix A(1000,3e5) . Is that
>> correct?
>>
> That is incorrect. The matrix is of size (N, N) with N = 1000 * 3e^5. It
> is a square
> matrix coming from an implicit scheme.
>
> Since the other answer appears to have the same misunderstanding, let me
> try
> to re-explain my point:
> --- In many contexts we need a partial differential equation that is not
> scalar.
> For example, the shallow water equation has b = 3 fields: H, HU, HV. There
> are other
> examples like wave modelling where we have something like b = 1000 fields
> (in a
> discretization).
> --- So, if we work with say an unstructured grid with N nodes then the
> total number
> of variables of the system will be N_tot = 3N or N_tot = 1000N.
>
> The linear system has N_tot unknowns and N_tot equations. The entries
> can be written as idx = (i , j) with 1 <= i <= b   and   1 <= j <= N.
>

This notation does not make sense to me. You are labeling an unknown with
an index (i, j) where i in [1, b] and j in [1, N].
Usually, we label an unknown in a matrix using (i, j) where i, j in [1, N].
Perhaps what you want is a multiindex

  alpha = (i, k) where i in [1, N] and k in [1, b]

Then a matrix entry is indexed using

  (alpha, beta) = ((i, k), (j, l))

As Pierre points out, this is exactly what we do in MATBAIJ, where we call
"i" a "block index". If you only have block structure,
then BAIJ is the best you can do. If you have product structure, then KAIJ
is better.

  Thanks,

     Matt


> Thus the non-zero entries in the matrix will be of two kinds:
> --- (idx1, idx2)   with idx1 = (i , j) and idx2 = (i' , j) , 1 <= i, i' <=
> b and 1 <= j <= N.
> Together those define a block in the matrix.
>
> --- (idx1, idx2) with idx1 = (i , j) and idx2 = (i, j'), 1<= i <= b and
> 1<= j, j' <= N.
> For each unknown idx1, there will be about 6 unknowns idx2 of this form.
>
> Otherwise, the block matrices do not have the same coefficients, so a
> tensor
> product approach does not appear to be workable.
>
>   Mathieu
>
>>

-- 
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.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210311/15185e7d/attachment.html>