[MITgcm-support] inverting a large penta-diagonal matrix to solve for magnetic field
David Trossman
david.s.trossman at gmail.com
Tue Sep 21 12:14:27 EDT 2021
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/20210921/6945ef32/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/20210921/6945ef32/attachment-0002.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/20210921/6945ef32/attachment-0003.png>
More information about the MITgcm-support
mailing list