On Thu, Mar 11, 2021 at 8:49 AM Mathieu Dutour <[email protected]> wrote:
> On Thu, 11 Mar 2021 at 13:52, Mark Adams <[email protected]> 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/>
