[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