/[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.8 - (hide annotations) (download)
Tue May 23 16:24:46 2017 UTC (8 years, 1 month ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, HEAD
Changes since 1.7: +12 -4 lines
update description how to generate tangent linear code for exact
jacobian times vector after adding new subroutines

1 mlosch 1.8 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_jacvec.F,v 1.7 2016/04/22 08:50:34 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 mlosch 1.5 #ifdef SEAICE_ALLOW_JFNK
60 mlosch 1.1 C Local variables:
61     _RL utp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
62     _RL vtp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
63     C u/vIceResP :: residual computed with u/vtp
64     _RL uIceResP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
65     _RL vIceResP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
66    
67     C i,j,bi,bj :: loop indices
68     INTEGER i,j,bi,bj
69     _RL epsilon, reps
70     CEOP
71 mlosch 1.2 C Instructions for using TAF or TAMC to generate exact Jacobian times
72 mlosch 1.8 C vector operations (if SEAICE_ALLOW_MOM_ADVECTION is defined, the
73     C file list also needs to include seaice_mom_advection.f,
74     C mom_calc_hfacz.f, mom_calc_ke.f, mom_calc_relvort3.f,
75     C mom_vi_u_coriolis.f, mom_vi_u_coriolis_c4.f, mom_vi_u_grad_ke.f,
76     C mom_vi_v_coriolis.f, mom_vi_v_coriolis_c4.f, mom_vi_v_grad_ke.f
77     C plus flow information for diagnostics_fill.f:
78     CCCCCCCADJ SUBROUTINE DIAGNOSTICS_FILL INPUT = 1,2,3,4,5,6,7,8
79     CCCCCCCADJ SUBROUTINE DIAGNOSTICS_FILL OUTPUT =
80     C )
81 mlosch 1.2 C
82     C 1. make small_f
83 mlosch 1.6 C 2. cat seaice_calc_residual.f seaice_oceandrag_coeffs.f \
84 mlosch 1.8 C seaice_bottomdrag_coeffs.f seaice_calc_stressdiv.f \
85 mlosch 1.6 C seaice_calc_strainrates.f seaice_calc_viscosities.f \
86     C seaice_calc_rhs.f seaice_calc_lhs.f > taf_input.f
87     C 3. staf -v1 -forward -toplevel seaice_calc_residual \
88     C -input uIceLoc,viceLoc -output uIceRes,vIceRes taf_input.f
89 mlosch 1.2 C 4. insert content of taf_input_ftl.f at the end of this file
90     C 5. add the following code and comment out the finite difference code
91 mlosch 1.4 C
92     C Instruction for using TAF 2.4 and higher (or staf with default -v2
93     C starting with version 2.0):
94     C
95     C 1. make small_f
96 mlosch 1.6 C 2. files="seaice_calc_residual.f seaice_oceandrag_coeffs.f \
97 mlosch 1.8 C seaice_bottomdrag_coeffs.f seaice_calc_stressdiv.f \
98 mlosch 1.6 C seaice_calc_strainrates.f seaice_calc_viscosities.f \
99     C seaice_calc_rhs.f seaice_calc_lhs.f"
100     C 3. staf -forward -toplevel seaice_calc_residual \
101     C -input uIceLoc,viceLoc -output uIceRes,vIceRes $files
102     C 4. copy files seaice_*_tl.f to the corresponding seaice_*.f files,
103 mlosch 1.4 C e.g. with this bash script:
104 mlosch 1.6 C for file in $files; do
105     C nfile=`echo $file | awk -F. '{printf "%s_tl.f", $1}'`;
106     C \cp -f $nfile $file
107 mlosch 1.4 C done
108 mlosch 1.6 C 5. add the following code, change "call g_seaice_calc_residual"
109     C to "call seaice_calc_residual_tl", and comment out the finite
110     C difference code
111 mlosch 1.2 CML _RL g_duIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
112     CML _RL g_dvIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
113     CML _RL g_uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
114     CML _RL g_vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
115     CML
116     CMLC Initialise
117     CML DO bj=myByLo(myThid),myByHi(myThid)
118     CML DO bi=myBxLo(myThid),myBxHi(myThid)
119     CML DO J=1-Oly,sNy+Oly
120     CML DO I=1-Olx,sNx+Olx
121     CML g_duIce(I,J,bi,bj) = duice(I,J,bi,bj)
122     CML g_dvIce(I,J,bi,bj) = dvice(I,J,bi,bj)
123     CML g_uIceRes(I,J,bi,bj) = 0. _d 0
124     CML g_vIceRes(I,J,bi,bj) = 0. _d 0
125     CML uIceResP(I,J,bi,bj) = 0. _d 0
126     CML vIceResP(I,J,bi,bj) = 0. _d 0
127     CML ENDDO
128     CML ENDDO
129     CML ENDDO
130     CML ENDDO
131     CML
132     CML CALL G_SEAICE_CALC_RESIDUAL( uIce, g_duice, vIce,
133     CML $g_dvice, uiceresp, g_uiceres, viceresp, g_viceres, newtoniter,
134     CML $kryloviter, mytime, myiter, mythid )
135 mlosch 1.4 CMLCML For staf -v2 replace the above with the below call
136     CMLCML CALL SEAICE_CALC_RESIDUAL_TL( uIce, g_duice, vIce,
137     CMLCML $g_dvice, uiceresp, g_uiceres, viceresp, g_viceres, newtoniter,
138     CMLCML $kryloviter, mytime, myiter, mythid )
139 mlosch 1.2 CML
140     CML DO bj=myByLo(myThid),myByHi(myThid)
141     CML DO bi=myBxLo(myThid),myBxHi(myThid)
142     CML DO J=1-Oly,sNy+Oly
143     CML DO I=1-Olx,sNx+Olx
144     CML duice(I,J,bi,bj)=g_uiceres(I,J,bi,bj)
145     CML dvice(I,J,bi,bj)=g_viceres(I,J,bi,bj)
146     CML ENDDO
147     CML ENDDO
148     CML ENDDO
149     CML ENDDO
150    
151 mlosch 1.1 C Initialise
152 mlosch 1.4 epsilon = SEAICE_JFNKepsilon
153 mlosch 1.1 reps = 1. _d 0/epsilon
154    
155     DO bj=myByLo(myThid),myByHi(myThid)
156     DO bi=myBxLo(myThid),myBxHi(myThid)
157     DO J=1-Oly,sNy+Oly
158     DO I=1-Olx,sNx+Olx
159     utp(I,J,bi,bj) = uIce(I,J,bi,bj) + epsilon * duIce(I,J,bi,bj)
160     vtp(I,J,bi,bj) = vIce(I,J,bi,bj) + epsilon * dvIce(I,J,bi,bj)
161     ENDDO
162     ENDDO
163     ENDDO
164     ENDDO
165    
166     C Compute new residual F(u)
167     CALL SEAICE_CALC_RESIDUAL(
168     I utp, vtp,
169     O uIceResP, vIceResP,
170     I newtonIter, krylovIter, myTime, myIter, myThid )
171    
172     C approximate Jacobian times vector by one-sided finite differences
173     C and store in du/vIce
174     DO bj = myByLo(myThid),myByHi(myThid)
175     DO bi = myBxLo(myThid),myBxHi(myThid)
176     DO I = 1, sNx
177     DO J = 1, sNy
178     duIce(I,J,bi,bj) =
179     & (uIceResP(I,J,bi,bj)-uIceRes(I,J,bi,bj))*reps
180     dvIce(I,J,bi,bj) =
181     & (vIceResP(I,J,bi,bj)-vIceRes(I,J,bi,bj))*reps
182     ENDDO
183     ENDDO
184     ENDDO
185     ENDDO
186    
187 mlosch 1.5 #endif /* SEAICE_ALLOW_JFNK */
188 mlosch 1.1
189     RETURN
190     END

  ViewVC Help
Powered by ViewVC 1.1.22