/[MITgcm]/MITgcm/pkg/seaice/seaice_jacvec.F
ViewVC logotype

Annotation of /MITgcm/pkg/seaice/seaice_jacvec.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (hide annotations) (download)
Tue Nov 6 13:09:30 2012 UTC (12 years, 8 months ago) by mlosch
Branch: MAIN
Changes since 1.1: +46 -2 lines
add instructions to use TAF to generate jacobian-times-vector code and
replace finite-difference code with it

1 mlosch 1.2 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_jacvec.F,v 1.1 2012/10/16 07:00:21 mlosch Exp $
2 mlosch 1.1 C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: SEAICE_JACVEC
8     C !INTERFACE:
9     SUBROUTINE SEAICE_JACVEC(
10     I uIceLoc, vIceLoc, uIceRes, vIceRes,
11     U duIce, dvIce,
12     I newtonIter, krylovIter, myTime, myIter, myThid )
13    
14     C !DESCRIPTION: \bv
15     C *==========================================================*
16     C | SUBROUTINE SEAICE_JACVEC
17     C | o For Jacobian-free Newton-Krylov solver compute
18     C | Jacobian times vector by finite difference approximation
19     C *==========================================================*
20     C | written by Martin Losch, Oct 2012
21     C *==========================================================*
22     C \ev
23    
24     C !USES:
25     IMPLICIT NONE
26    
27     C === Global variables ===
28     #include "SIZE.h"
29     #include "EEPARAMS.h"
30     #include "PARAMS.h"
31     #include "DYNVARS.h"
32     #include "GRID.h"
33     #include "SEAICE_SIZE.h"
34     #include "SEAICE_PARAMS.h"
35     #include "SEAICE.h"
36    
37     C !INPUT/OUTPUT PARAMETERS:
38     C === Routine arguments ===
39     C myTime :: Simulation time
40     C myIter :: Simulation timestep number
41     C myThid :: my Thread Id. number
42     C newtonIter :: current iterate of Newton iteration
43     C krylovIter :: current iterate of Krylov iteration
44     _RL myTime
45     INTEGER myIter
46     INTEGER myThid
47     INTEGER newtonIter
48     INTEGER krylovIter
49     C u/vIceLoc :: local copies of the current ice velocity
50     _RL uIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
51     _RL vIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
52     C u/vIceRes :: initial residual of this Newton iterate
53     _RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54     _RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
55     C du/vIce :: correction of ice velocities
56     _RL duIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
57     _RL dvIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
58    
59     #if ( defined (SEAICE_CGRID) && \
60     defined (SEAICE_ALLOW_JFNK) && \
61     defined (SEAICE_ALLOW_DYNAMICS) )
62     C Local variables:
63     _RL utp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
64     _RL vtp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
65     C u/vIceResP :: residual computed with u/vtp
66     _RL uIceResP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
67     _RL vIceResP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
68    
69     C i,j,bi,bj :: loop indices
70     INTEGER i,j,bi,bj
71     _RL epsilon, reps
72     CEOP
73 mlosch 1.2 C Instructions for using TAF or TAMC to generate exact Jacobian times
74     C vector operations:
75     C
76     C 1. make small_f
77     C 2. cat seaice_calc_residual.f seaice_oceandrag_coeffs.f seaice_calc_strainrates.f seaice_calc_viscosities.f seaice_calc_rhs.f seaice_calc_lhs.f > taf_input.f
78     C 3. staf -v1 -forward -toplevel seaice_calc_residual -input uIceLoc,viceLoc -output uIceRes,vIceRes taf_input.f
79     C 4. insert content of taf_input_ftl.f at the end of this file
80     C 5. add the following code and comment out the finite difference code
81     CML _RL g_duIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
82     CML _RL g_dvIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
83     CML _RL g_uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
84     CML _RL g_vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
85     CML
86     CML IF (.false.) then
87     CMLC Initialise
88     CML DO bj=myByLo(myThid),myByHi(myThid)
89     CML DO bi=myBxLo(myThid),myBxHi(myThid)
90     CML DO J=1-Oly,sNy+Oly
91     CML DO I=1-Olx,sNx+Olx
92     CML g_duIce(I,J,bi,bj) = duice(I,J,bi,bj)
93     CML g_dvIce(I,J,bi,bj) = dvice(I,J,bi,bj)
94     CML g_uIceRes(I,J,bi,bj) = 0. _d 0
95     CML g_vIceRes(I,J,bi,bj) = 0. _d 0
96     CML uIceResP(I,J,bi,bj) = 0. _d 0
97     CML vIceResP(I,J,bi,bj) = 0. _d 0
98     CML ENDDO
99     CML ENDDO
100     CML ENDDO
101     CML ENDDO
102     CML
103     CML CALL G_SEAICE_CALC_RESIDUAL( uIce, g_duice, vIce,
104     CML $g_dvice, uiceresp, g_uiceres, viceresp, g_viceres, newtoniter,
105     CML $kryloviter, mytime, myiter, mythid )
106     CML
107     CML DO bj=myByLo(myThid),myByHi(myThid)
108     CML DO bi=myBxLo(myThid),myBxHi(myThid)
109     CML DO J=1-Oly,sNy+Oly
110     CML DO I=1-Olx,sNx+Olx
111     CML duice(I,J,bi,bj)=g_uiceres(I,J,bi,bj)
112     CML dvice(I,J,bi,bj)=g_viceres(I,J,bi,bj)
113     CML ENDDO
114     CML ENDDO
115     CML ENDDO
116     CML ENDDO
117    
118 mlosch 1.1 C Initialise
119     epsilon = 1. _d -06
120     reps = 1. _d 0/epsilon
121    
122     DO bj=myByLo(myThid),myByHi(myThid)
123     DO bi=myBxLo(myThid),myBxHi(myThid)
124     DO J=1-Oly,sNy+Oly
125     DO I=1-Olx,sNx+Olx
126     utp(I,J,bi,bj) = uIce(I,J,bi,bj) + epsilon * duIce(I,J,bi,bj)
127     vtp(I,J,bi,bj) = vIce(I,J,bi,bj) + epsilon * dvIce(I,J,bi,bj)
128     ENDDO
129     ENDDO
130     ENDDO
131     ENDDO
132    
133     C Compute new residual F(u)
134     CALL SEAICE_CALC_RESIDUAL(
135     I utp, vtp,
136     O uIceResP, vIceResP,
137     I newtonIter, krylovIter, myTime, myIter, myThid )
138    
139     C approximate Jacobian times vector by one-sided finite differences
140     C and store in du/vIce
141     DO bj = myByLo(myThid),myByHi(myThid)
142     DO bi = myBxLo(myThid),myBxHi(myThid)
143     DO I = 1, sNx
144     DO J = 1, sNy
145     duIce(I,J,bi,bj) =
146     & (uIceResP(I,J,bi,bj)-uIceRes(I,J,bi,bj))*reps
147     dvIce(I,J,bi,bj) =
148     & (vIceResP(I,J,bi,bj)-vIceRes(I,J,bi,bj))*reps
149     ENDDO
150     ENDDO
151     ENDDO
152     ENDDO
153    
154     #endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID and SEAICE_ALLOW_JFNK */
155    
156     RETURN
157     END

  ViewVC Help
Powered by ViewVC 1.1.22