[MITgcm-support] inverting a large penta-diagonal matrix to solve for magnetic field
David Trossman
david.s.trossman at gmail.com
Mon Sep 27 09:55:43 EDT 2021
This is helpful, Martin. I will take a look at how the mapping is done for
sea ice and use that as a guide. I will let you know if I have any
additional questions and will let you know if I solve the issue.
-David
On Mon, Sep 27, 2021 at 5:00 AM Martin Losch <Martin.Losch at awi.de> wrote:
> 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
>
> _______________________________________________
> 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/20210927/56e68e3e/attachment-0001.html>
More information about the MITgcm-support
mailing list