[MITgcm-support] inverting a large penta-diagonal matrix to solve for magnetic field
Martin Losch
Martin.Losch at awi.de
Mon Sep 27 05:59:25 EDT 2021
HI David,
I am not entirely sure how your solution algorithm works, but we have something similar, e.g. in the sea ice model. where we use a global Krylov method (pkg/seaice/seaice_krylov.F), or a Jacobian free Newton Krylov method (pkg/seaice/seaice_jfnk.F). In all cases, the 2d fields are copied to 1D vectors (s/r seaice_map2vec in pkg/seaice/seaice_fgmres.F). The vector is then still split between the tiles, i.e. the vector field has the dimensions (sNx*sNy, nSx, nSy), and all operations are local to the tile and processor, except for the scalar product (also defined in seaice_fmgres.F), where there is a global sum.
In your algorithm, you seem to need some overlap between tiles, so your vector should maybe have the dimensions ((sNx+2)*(sNy+2), nSx, nSy) with the appropriate mapping (i=1-1,sNx+1, etc), but I can’t be sure about that.
Does that make sense?
Martin
> On 22. Sep 2021, at 22:33, David Trossman <david.s.trossman at gmail.com> wrote:
>
> I solved the penta-diagonal matrix inversion algorithm issue, but now the solution appears to have its tiles out of order:
> <image.png>
> If I want to place an array in i,j,bi,bj dimensions into an array with a single dimension (e.g., the diagonal of a matrix), perform my (matrix) operations, and put the result back in i,j,bi,bj dimensions, how should I order the array with a single dimension? I think I need to order it by cycling over the longitudes for a given latitude and then repeat for each latitude. Do I need to manually code this or is there an equivalent procedure in the code already for a different application? I hope I'm being clear (my forever-problem)... Any advice would be appreciated. Thanks,
> David
>
> On Tue, Sep 21, 2021 at 11:14 AM David Trossman <david.s.trossman at gmail.com> wrote:
> Thanks for the response, Martin. Using myBxLo(myThid),myBxHi(myThid) and myByLo(myThid),myByHi(myThid) instead of 1,nSx and 1,nSy doesn't make a difference. The code I sent you happened to use the latter because I suspected that it would make a difference too, but it doesn't. Ngrid definitely has the correct value. Also, I figured that I should show you an example input to the penta-diagonal solver subroutine.
>
> Here is the forcing term (which is variable B in the subroutine code I sent you):
> <image.png>
> I've n-ple-checked that this is correct. Here is the diagonal entry of the penta-diagonal matrix (which is variable D in the subroutine code I sent you):
> <image.png>
> You can see that this is also correct, as it should correlate highly with the seafloor depths.
>
> You mentioned doing a run without blank tiles. I didn't know this was an option. It looked to me like I just needed to comment out blankList in data.exch2. I tried this and I still see the "hole" over Asia. What else do I need to do to try this? I suspect that the "hole" isn't the main problem, though. My main problem is my penta-diagonal matrix solver subroutine. I can't just invert the matrix using some LU decomposition subroutine (too much memory). So that leaves me stuck because I don't expect anyone to know how to correctly program the very specific type of matrix solver I need in Fortran...
> Cheers,
> David
>
> On Tue, Sep 21, 2021 at 5:17 AM Martin Losch <Martin.Losch at awi.de> wrote:
> Hi David,
>
> have you tried to do this in a run without blank tiles? That way the “hole” over Asia has some reasonable coordinates (and values=0 instead of NaN).
> Also you may want to use myBxLo(myThid),myBxHi(myThid) instead of 1,nSx and similar for y?
>
> I am assuming that Ngrid has the correct value. Do you need the overlap here? You I,J loop ranges include the overlaps, but are these required for the solution (as they are just duplicates from the neighboring tiles). If you do need overlaps, then the code implies that you only need one and not OLx/OLy?
>
> Martin
>
> > On 14. Sep 2021, at 18:01, David Trossman <david.s.trossman at gmail.com> wrote:
> >
> > Hi all,
> > I've been banging my head against the wall with this problem for so long I may have too much brain damage to come up with a solution.
> >
> > The problem is this: I want to solve some equations to solve for the magnetic field associated with the ocean circulation, and these equations involve the inverse of a penta-diagonal matrix. I translated what would be the following in Matlab (I know, I know...):
> >
> > difFluxw = spdiags([-Cw(:),[Cw(2:end);nan]], [0, -1], Ngrid,Ngrid);
> > difFluxe = spdiags([-Ce(:),[nan;Ce(1:end-1)]], [0, 1], Ngrid,Ngrid);
> > difFluxs = spdiags([-Cs(:),[Cs(1+Nlon:end);nan*ones(Nlon,1)]],[0, -Nlon], Ngrid,Ngrid);
> > difFluxn = spdiags([-Cn(:),[nan*ones(Nlon,1);Cn(1:end-Nlon)]],[0, Nlon], Ngrid,Ngrid);
> > advFluxw = spdiags([-aCw(:),-[aCw(2:end);nan]], [0, -1], Ngrid,Ngrid);
> > advFluxe = spdiags([ aCe(:), [nan;aCe(1:end-1)]], [0, 1], Ngrid,Ngrid);
> > advFluxs = spdiags([-aCs(:),-[aCs(1+Nlon:end);nan*ones(Nlon,1)]],[0, -Nlon], Ngrid,Ngrid);
> > advFluxn = spdiags([ aCn(:), [nan*ones(Nlon,1);aCn(1:end-Nlon)]],[0, Nlon], Ngrid,Ngrid);
> > Fluxw = difFluxw + advFluxw;
> > Fluxe = difFluxe + advFluxe;
> > Fluxs = difFluxs + advFluxs;
> > Fluxn = difFluxn + advFluxn;
> > bjunk=find(lon==min(lon(:)) & lat>min(lat(:)) & lat<max(lat(:)));
> > Fluxw(bjunk,:) = 0;
> > Fluxw = Fluxw + sparse(bjunk,bjunk, -Cw(bjunk), Ngrid,Ngrid);
> > Fluxw = Fluxw + sparse(bjunk,bjunk-1+Nlon, Cw(bjunk), Ngrid,Ngrid);
> > % flux through east side of cell:
> > bjunk=find(lon==max(lon(:)) & lat > min(lat(:)) & lat<max(lat(:)));
> > Fluxe(bjunk,:) = 0;
> > Fluxe = Fluxe + sparse(bjunk,bjunk, -Ce(bjunk), Ngrid,Ngrid);
> > Fluxe = Fluxe + sparse(bjunk,bjunk+1-Nlon, Ce(bjunk),Ngrid,Ngrid);
> > integHcoef = spdiags(Hcoef(:).*vol(:), 0, Ngrid,Ngrid);
> > metrik.Lfluxes = Fluxw + Fluxe + Fluxs + Fluxn ;
> > metrik.L = metrik.Lfluxes + integHcoef;
> >
> > into Fortran. Nevermind what the variables mean; all that's important is that the variable metrik.L is a pentadiagonal matrix that I save in Fortran as two off-diagonal arrays that are Ngrid-Nlon in length, two off-diagonal arrays that are Ngrid-1 in length, and a diagonal array of length Ngrid. I've written a pent-diagonal matrix solver for this type of situation as follows:
> >
> > SUBROUTINE PENTA(myThid,Ngrid,Nlon,E0,A0,D0,C0,F0,B0,X0)
> >
> > c RESULTS: matrix has 5 bands, EADCF, with D being the main diagonal,
> > c E and A are the lower diagonals, and C and F are the upper diagonals.
> >
> > c E is defined for rows i = n1:Ngrid, but is defined as E(1) to E(Ngrid-Nlon)
> > c A is defined for rows i = 2:Ngrid, but is defined as A(1) to A(Ngrid-1)
> > c D is defined for rows i = 1:Ngrid
> > c C is defined for rows i = 1:Ngrid-1, but the last element isn't used
> > c F is defined for rows i = 1:(Ngrid-Nlon), but the last Nlon-1 elements aren't used
> >
> > c B is the right-hand side
> > c X is the solution vector
> >
> > IMPLICIT NONE
> > #include "SIZE.h"
> > #include "EEPARAMS.h"
> > #include "GRID.h"
> > #include "SURFACE.h"
> >
> > INTEGER, intent(in) :: myThid
> > INTEGER, intent(in) :: Ngrid
> > INTEGER, intent(in) :: Nlon
> > _RL, dimension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nSx,1:nSy),
> > & intent(in) :: E0
> > _RL, dimension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nSx,1:nSy),
> > & intent(in) :: A0
> > _RL, dimension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nSx,1:nSy),
> > & intent(in) :: D0
> > _RL, dimension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nSx,1:nSy),
> > & intent(in) :: C0
> > _RL, dimension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nSx,1:nSy),
> > & intent(in) :: F0
> > _RL, dimension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nSx,1:nSy),
> > & intent(out) :: B0
> > _RL, dimension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nSx,1:nSy),
> > & intent(out) :: X0
> >
> > INTEGER bi,bj
> > INTEGER I,J,IJ
> > INTEGER Nn
> > _RL XMULT
> > _RL, dimension(1:Ngrid*nSx*nSy) :: E
> > _RL, dimension(1:Ngrid*nSx*nSy) :: A
> > _RL, dimension(1:Ngrid*nSx*nSy) :: D
> > _RL, dimension(1:Ngrid*nSx*nSy) :: C
> > _RL, dimension(1:Ngrid*nSx*nSy) :: F
> > _RL, dimension(1:Ngrid*nSx*nSy) :: B
> > _RL, dimension(1:Ngrid*nSx*nSy) :: X
> >
> > Nn=Ngrid*nSx*nSy
> > IJ=0
> > DO bj=1,nSy
> > DO bi=1,nSx
> > DO I = 1-OLx,sNx+OLx
> > DO J = 1-OLy,sNy+OLy
> > IJ=IJ+1
> > E(IJ)=E0(I,J,bi,bj)
> > A(IJ)=A0(I,J,bi,bj)
> > D(IJ)=D0(I,J,bi,bj)
> > C(IJ)=C0(I,J,bi,bj)
> > F(IJ)=F0(I,J,bi,bj)
> > ENDDO
> > ENDDO
> > ENDDO
> > ENDDO
> > DO IJ = 2,Nn-1
> > XMULT = A(IJ-1)/D(IJ-1)
> > D(IJ) = D(IJ) - XMULT*C(IJ-1)
> > C(IJ) = C(IJ) - XMULT*F(IJ-1)
> > B(IJ) = B(IJ) - XMULT*B(IJ-1)
> > XMULT = E(IJ-1)/D(IJ-1)
> > A(IJ) = A(IJ) - XMULT*C(IJ-1)
> > D(IJ+1) = D(IJ+1) - XMULT*F(IJ-1)
> > B(IJ+1) = B(IJ+1) - XMULT*B(IJ-1)
> > ENDDO
> > XMULT = A(Nn-1)/D(Nn-1)
> > D(Nn) = D(Nn) - XMULT*C(Nn-1)
> > X(Nn) = (B(Nn) - XMULT*B(Nn-1))/D(Nn)
> > DO IJ = Nn-1,Nn-Nlon+1,-1
> > X(IJ) = (B(IJ) - C(IJ)*X(IJ+1))/D(IJ)
> > ENDDO
> > DO IJ = Nn-Nlon,1,-1
> > X(IJ) = (B(IJ) - F(IJ)*X(IJ+n1) -
> > & C(IJ)*X(IJ+1))/D(IJ)
> > ENDDO
> > IJ=0
> > DO bj=1,nSy
> > DO bi=1,nSx
> > DO I = 1-OLx,sNx+OLx
> > DO J = 1-OLy,sNy+OLy
> > IJ=IJ+1
> > B0(I,J,bi,bj)=B(IJ)
> > X0(I,J,bi,bj)=X(IJ)
> > ENDDO
> > ENDDO
> > ENDDO
> > ENDDO
> > RETURN
> > END
> >
> > X0 is what I want, but the above code isn't correct. I get something like the following for the magnetic field:
> > <image.png>
> > It seems like I'm not considering an ordering for the tiling correctly, or the pent-diagonal matrix solver isn't quite correct in other ways (although, I'm adapting it from an existing code with two off-diagonal arrays of length Ngrid-2, two off-diagonal arrays of length Ngrid-1, and a diagonal array of length Ngrid: http://www.math.uakron.edu/~kreider/anpde/penta.f). It's also possible that I need to interpolate over the "hole" (of NaNs) over Asia in the MITgcm. I recognize that you don't know whether my inputs to the subroutine are correct, but just assume they are.
> >
> > One possible approach to figuring out what's wrong is to interpolate from the LLC grid to a regularly spaced lat-lon grid, but I'm not sure if this capability exists inline as the model runs. No, I don't want to do this offline using averaged output because I'm trying to make this calculation adjoint-able (plus, a function of the average of each field is not the same as the average of a function of each field). I figured that someone on this list would be able to point out what I'm doing incorrectly in the above Fortran code or knows how to do the interpolation to a regular lat-lon grid, but maybe I'm asking for too much... In any case, thanks for any help.
> > -David
> > _______________________________________________
> > MITgcm-support mailing list
> > MITgcm-support at mitgcm.org
> > http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support
-------------- next part --------------
A non-text attachment was scrubbed...
Name: smime.p7s
Type: application/pkcs7-signature
Size: 1665 bytes
Desc: not available
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20210927/b96a3764/attachment.p7s>
More information about the MITgcm-support
mailing list