[MITgcm-support] inverting a large penta-diagonal matrix to solve for magnetic field
David Trossman
david.s.trossman at gmail.com
Tue Sep 14 12:01:22 EDT 2021
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: 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20210914/f5872506/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 998486 bytes
Desc: not available
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20210914/f5872506/attachment-0001.png>
More information about the MITgcm-support
mailing list