[MITgcm-support] inverting a large penta-diagonal matrix to solve for magnetic field
David Trossman
david.s.trossman at gmail.com
Thu Sep 30 16:22:42 EDT 2021
Using seaice_map2vec in pkg/seaice/seaice_fgmres.F as a guide, I was able
to resolve the problem. The remaining issue now is the "forcing" field in
the equations I'm solving, which I should be able to figure out by myself.
Thanks!
-David
On Mon, Sep 27, 2021 at 8:55 AM David Trossman <david.s.trossman at gmail.com>
wrote:
> 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/20210930/2c286c8c/attachment-0001.html>
More information about the MITgcm-support
mailing list