If it is not a BAIJ matrix, perhaps a figure of the matrix (for some small b and N) might help us understand the structure. We understand vector valued (non-scalar) PDEs and have worked with cases of 1000s of entries at each grid point but don't understand the index notation you are using below.
Also did you previously use an AIJ that resulted in poor performance? Barry > On Mar 11, 2021, at 7:48 AM, Mathieu Dutour <[email protected]> wrote: > > On Thu, 11 Mar 2021 at 13:52, Mark Adams <[email protected] > <mailto:[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. > > 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
