[MITgcm-support] inverting a large penta-diagonal matrix to solve for magnetic field

David Trossman david.s.trossman at gmail.com
Wed Sep 22 16:33:29 EDT 2021


I solved the penta-diagonal matrix inversion algorithm issue, but now the
solution appears to have its tiles out of order:
[image: 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: 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: 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
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20210922/3d0bf281/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 1137691 bytes
Desc: not available
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20210922/3d0bf281/attachment-0003.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 748546 bytes
Desc: not available
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20210922/3d0bf281/attachment-0004.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 1017881 bytes
Desc: not available
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20210922/3d0bf281/attachment-0005.png>


More information about the MITgcm-support mailing list