[MITgcm-support] inverting a large penta-diagonal matrix to solve for magnetic field
Martin Losch
Martin.Losch at awi.de
Tue Sep 21 06:00:48 EDT 2021
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
-------------- 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/20210921/8c1143d4/attachment.p7s>
More information about the MITgcm-support
mailing list