Actual source code: ex30.c
1: /*$Id: ex19.c,v 1.30 2001/08/07 21:31:17 bsmith Exp $*/
3: static char help[] = "Steady-state 2D subduction flow, pressure and temperature solver.\n\\n\
4: The flow is driven by the subducting slab.\n\
5: -ivisc <#> = rheology option.\n\
6: 0 --- constant viscosity.\n\
7: 1 --- olivine diffusion creep rheology (T-dependent, newtonian).\n\
8: 2 --- weak temperature dependent rheology (1200/T, newtonian).\n\
9: -ibound <#> = boundary condition \n\
10: 0 --- isoviscous analytic.\n\
11: 1 --- stress free. \n\
12: 2 --- stress is von neumann. \n\
13: -icorner <#> = i index of wedge corner point.\n\
14: -jcorner <#> = j index of wedge corner point.\n\
15: -slab_dip <#> = dip of the subducting slab in DEGREES.\n\
16: -back_arc <#> = distance from trench to back-arc in KM.(if unspecified then no back-arc). \n\
17: -u_back_arcocity <#> = full spreading rate of back arc as a factor of slab velocity. \n\
18: -width <#> = width of domain in KM.\n\
19: -depth <#> = depth of domain in KM.\n\
20: -lid_depth <#> = depth to the base of the lithosphere in KM.\n\
21: -slab_dip <#> = dip of the subducting slab in DEGREES.\n\
22: -slab_velocity <#> = velocity of slab in CM/YEAR.\n\
23: -slab_age <#> = age of slab in MILLIONS OF YEARS. \n\
24: -potentialT <#> = mantle potential temperature in degrees CENTIGRADE.\n\
25: -kappa <#> = thermal diffusivity in M^2/SEC. \n\
26: -peclet <#> = dimensionless Peclet number (default 111.691)\n\\n";
28: /*T
29: Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
30: Concepts: DA^using distributed arrays;
31: Concepts: multicomponent
32: Processors: n
33: T*/
35: /* ------------------------------------------------------------------------
37: This problem is modeled by the partial differential equation system
38:
39: -Grad(P) + Div[eta (Grad(v) + Grad(v)^T)] = 0
40: Div(U,W) = 0
42: which is uniformly discretized on a staggered mesh:
43: ------w_ij------
44: / /
45: u_i-1j P_ij u_ij
46: / /
47: ------w_ij-1----
49: ------------------------------------------------------------------------- */
51: /*
52: Include "petscda.h" so that we can use distributed arrays (DAs).
53: Include "petscsnes.h" so that we can use SNES solvers. Note that this
54: file automatically includes:
55: petsc.h - base PETSc routines petscvec.h - vectors
56: petscsys.h - system routines petscmat.h - matrices
57: petscis.h - index sets petscksp.h - Krylov subspace methods
58: petscviewer.h - viewers petscpc.h - preconditioners
59: petscksp.h - linear solvers
60: */
61: #include petscsnes.h
62: #include petscda.h
64: /*
65: User-defined routines and data structures
66: */
68: /*
69: The next two structures are essentially the same. The difference is that
70: the first does not contain derivative information (as used by ADIC) while the
71: second does. The first is used only to contain the solution at the previous time-step
72: which is a constant in the computation of the current time step and hence passive
73: (meaning not active in terms of derivatives).
74: */
76: typedef struct {
77: PetscScalar u,w,p,T;
78: } Field;
80: typedef struct {
81: int mglevels;
82: PetscReal width, depth, scaled_width, scaled_depth, peclet, potentialT;
83: PetscReal slab_dip, slab_age, slab_velocity, lid_depth, kappa;
84: PetscReal u_back_arc, x_back_arc;
85: int icorner, jcorner, ivisc, ifromm, ibound;
86: PetscTruth PreLoading, back_arc;
87: } Parameter;
89: typedef struct {
90: Vec Xold,func;
91: Parameter *param;
92: PetscReal fnorm_ini, fnorm;
93: } AppCtx;
95: PetscReal HALFPI = 3.14159265358979323846/2.0;
97: int SetParams(Parameter *param, int mx, int mz);
98: extern int FormInitialGuess(SNES,Vec,void*);
99: extern int FormFunctionLocal(DALocalInfo*,Field**,Field**,void*);
100: extern PassiveScalar HorizVelocity(PassiveScalar x, PassiveScalar z, PassiveScalar c, PassiveScalar d);
101: extern PassiveScalar VertVelocity(PassiveScalar x, PassiveScalar z, PassiveScalar c, PassiveScalar d);
102: extern PassiveScalar Pressure(PassiveScalar x, PassiveScalar z, PassiveScalar c, PassiveScalar d);
103: extern PassiveScalar LidVelocity(PassiveScalar x, PassiveScalar xBA, PetscTruth BA);
104: extern PetscScalar Viscosity(PetscScalar T, int iVisc);
105: extern PetscScalar UInterp(Field **x, int i, int j, PetscScalar fr);
106: extern PetscScalar WInterp(Field **x, int i, int j, PetscScalar fr);
107: extern PetscScalar PInterp(Field **x, int i, int j, PetscScalar fr);
108: extern PetscScalar TInterp(Field **x, int i, int j, PetscScalar fr);
109: extern void CalcJunk(PassiveScalar,PassiveScalar,PassiveScalar,PassiveScalar*,PassiveScalar*,PassiveScalar*,PassiveScalar*,PassiveScalar*);
110: int Initialize(DMMG*);
112: /*-----------------------------------------------------------------------*/
115: int main(int argc,char **argv)
116: /*-----------------------------------------------------------------------*/
117: {
118: DMMG *dmmg; /* multilevel grid structure */
119: AppCtx *user; /* user-defined work context */
120: Parameter param;
121: int mx,mz,its;
122: int i,ierr,tmpVisc,tmpBound;
123: MPI_Comm comm;
124: SNES snes;
125: DA da;
126: Vec res;
127: PetscReal SEC_PER_YR = 3600.00*24.00*356.2500;
129: PetscInitialize(&argc,&argv,(char *)0,help);
130: PetscOptionsSetValue("-preload","off");
131: PetscOptionsSetValue("-mat_type","seqaij");
132: PetscOptionsSetValue("-snes_monitor",PETSC_NULL);
133: PetscOptionsSetValue("-ksp_truemonitor",PETSC_NULL);
134: PetscOptionsSetValue("-pc_type","lu");
135: PetscOptionsInsert(&argc,&argv,PETSC_NULL);
137: comm = PETSC_COMM_WORLD;
138: PreLoadBegin(PETSC_TRUE,"SetUp");
140: param.PreLoading = PreLoading; //where is PreLoading defined?
141: DMMGCreate(comm,1,&user,&dmmg);
142: param.mglevels = DMMGGetLevels(dmmg);
144: /*
145: Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
146: for principal unknowns (x) and governing residuals (f)
147: */
148: DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,2,0,0,&da);
149: DMMGSetDM(dmmg,(DM)da);
150: DADestroy(da);
151: DAGetInfo(DMMGGetDA(dmmg),0,&mx,&mz,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
152: PETSC_IGNORE,PETSC_IGNORE);
154: /*
155: Problem parameters
156: */
157: SetParams(¶m,mx,mz);
159: DASetFieldName(DMMGGetDA(dmmg),0,"x-velocity");
160: DASetFieldName(DMMGGetDA(dmmg),1,"y-velocity");
161: DASetFieldName(DMMGGetDA(dmmg),2,"pressure");
162: DASetFieldName(DMMGGetDA(dmmg),3,"temperature");
164: /*======================================================================*/
165:
166: PetscMalloc(param.mglevels*sizeof(AppCtx),&user);
167: for (i=0; i<param.mglevels; i++) {
168: VecDuplicate(dmmg[i]->x, &(user[i].Xold));
169: VecDuplicate(dmmg[i]->x, &(user[i].func));
170: user[i].param = ¶m;
171: dmmg[i]->user = &user[i];
172: }
173: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: Create user context, set problem data, create vector data structures.
175: Also, compute the initial guess.
176: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177:
178: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: Create nonlinear solver context
180: Process NOT adiC(100): WInterp UInterp PInterp TInterp FormFunctionLocal Viscosity
181: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
182: /* DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);*/
183: DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,0,0);
184: DMMGSetInitialGuess(dmmg,FormInitialGuess);
185:
186: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187: Solve the nonlinear system
188: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189: PreLoadStage("Solve");
190: Initialize(dmmg);
192: if (param.ivisc>0) {
193: PetscPrintf(PETSC_COMM_WORLD,"Doing Constant Viscosity Solve\n", its);
194: tmpVisc = param.ivisc; param.ivisc=0; tmpBound = param.ibound; param.ibound=0;
195: snes = DMMGGetSNES(dmmg);
196: SNESSetTolerances(snes,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,2,PETSC_DEFAULT);
197: DMMGSolve(dmmg);
198: SNESSetTolerances(snes,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,50,PETSC_DEFAULT);
199: VecCopy(DMMGGetx(dmmg),user->Xold);
200: param.ivisc = tmpVisc; param.ibound=tmpBound;
201: PetscPrintf(PETSC_COMM_WORLD,"Doing Variable Viscosity Solve\n", its);
202: }
203: if (param.ifromm==1)
204: PetscPrintf(PETSC_COMM_WORLD,"Using Fromm advection scheme\n", its);
205: DMMGSolve(dmmg);
206: snes = DMMGGetSNES(dmmg);
207: SNESGetIterationNumber(snes,&its);
208: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %d\n", its);
209:
210: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
211: Output stuff to Matlab socket.
212: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
213: SNESGetFunction(snes, &res, PETSC_NULL, PETSC_NULL);
214: VecView(res, PETSC_VIEWER_SOCKET_WORLD);
215: VecView(DMMGGetx(dmmg),PETSC_VIEWER_SOCKET_WORLD);
216: PetscViewerSocketPutInt(PETSC_VIEWER_SOCKET_WORLD,1,&mx);
217: PetscViewerSocketPutInt(PETSC_VIEWER_SOCKET_WORLD,1,&mz);
218: PetscViewerSocketPutReal(PETSC_VIEWER_SOCKET_WORLD,1,1,&(param.slab_dip));
219: PetscViewerSocketPutReal(PETSC_VIEWER_SOCKET_WORLD,1,1,&(param.scaled_width));
220: PetscViewerSocketPutReal(PETSC_VIEWER_SOCKET_WORLD,1,1,&(param.scaled_depth));
221: PetscViewerSocketPutReal(PETSC_VIEWER_SOCKET_WORLD,1,1,&(param.lid_depth));
222: PetscViewerSocketPutReal(PETSC_VIEWER_SOCKET_WORLD,1,1,&(param.potentialT));
223: PetscViewerSocketPutReal(PETSC_VIEWER_SOCKET_WORLD,1,1,&(param.x_back_arc));
224: PetscViewerSocketPutInt(PETSC_VIEWER_SOCKET_WORLD,1,&(param.icorner));
225: PetscViewerSocketPutInt(PETSC_VIEWER_SOCKET_WORLD,1,&(param.jcorner));
226: PetscViewerSocketPutInt(PETSC_VIEWER_SOCKET_WORLD,1,&(param.ivisc));
227: PetscViewerSocketPutInt(PETSC_VIEWER_SOCKET_WORLD,1,&(param.ibound));
228: PetscViewerSocketPutInt(PETSC_VIEWER_SOCKET_WORLD,1,&(its));
230: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231: Free work space. All PETSc objects should be destroyed when they
232: are no longer needed.
233: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
234: for (i=0; i<param.mglevels; i++) {
235: VecDestroy(user[i].Xold);
236: VecDestroy(user[i].func);
237: }
238: PetscFree(user);
239: DMMGDestroy(dmmg);
240: PreLoadEnd();
241:
242: PetscFinalize();
243: return 0;
244: }
246: /* ------------------------------------------------------------------- */
249: /* ------------------------------------------------------------------- */
250: int Initialize(DMMG *dmmg)
251: {
252: AppCtx *user = (AppCtx*)dmmg[dmmg[0]->nlevels-1]->user;
253: Parameter *param = user->param;
254: DA da;
255: int i,j,mx,mz,ierr,xs,ys,xm,ym,iw,jw;
256: int mglevel,ic,jc;
257: PetscReal dx, dz, c, d, beta, itb, cb, sb, sbi, skt;
258: PetscReal xp, zp, r, st, ct, th, fr, xPrimeSlab;
259: Field **x;
261: da = (DA)(dmmg[param->mglevels-1]->dm); /* getting the fine grid */
263: beta = param->slab_dip;
264: CalcJunk(param->kappa,param->slab_age,beta,&skt,&cb,&sb,&itb,&sbi);
265: c = beta*sb/(beta*beta-sb*sb);
266: d = (beta*cb-sb)/(beta*beta-sb*sb);
268: DAGetInfo(da,0,&mx,&mz,0,0,0,0,0,0,0,0);
269: dx = param->scaled_width/((PetscReal)(mx-2));
270: dz = param->scaled_depth/((PetscReal)(mz-2));
272: ic = param->icorner; jc = param->jcorner;
273: xPrimeSlab = (ic-1)*dx;
275: fr = (1.0-dz/dx*itb)/2.0;
276: PetscPrintf(PETSC_COMM_WORLD,"interpolation fraction = %g\n", fr);
277: if (fr<0.0)
278: SETERRQ(1," USER ERROR: Grid shear exceeds limit! Decrease dz/dx!");
280: /*
281: Get local grid boundaries (for 2-dimensional DA):
282: xs, ys - starting grid indices (no ghost points)
283: xm, ym - scaled_widths of local grid (no ghost points)
284: */
285: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
287: /*
288: Get a pointer to vector data.
289: - For default PETSc vectors, VecGetArray() returns a pointer to
290: the data array. Otherwise, the routine is implementation dependent.
291: - You MUST call VecRestoreArray() when you no longer need access to
292: the array.
293: */
294: DAVecGetArray(da,((AppCtx*)dmmg[param->mglevels-1]->user)->Xold,(void**)&x);
296: /*
297: Compute initial guess (analytic soln) over the locally owned part of the grid
298: Initial condition is isoviscous corner flow and uniform temperature in interior
299: */
300: for (j=ys; j<ys+ym; j++) {
301: jw = j - jc + 1;
302: for (i=xs; i<xs+xm; i++) {
303: iw = i - ic + 1;
305: if (i<ic) { /* slab */
306: x[j][i].p = 0.0;
307: x[j][i].u = cb;
308: x[j][i].w = sb;
309: xp = (i-0.5)*dx; zp = (xPrimeSlab - xp)*tan(beta);
310: x[j][i].T = erf(zp*param->lid_depth/2.0/skt);
311: } else if (j<jc) { /* lid */
312: x[j][i].p = 0.0;
313: x[j][i].u = 0.0;
314: x[j][i].w = 0.0;
315: zp = (j-0.5)*dz;
316: x[j][i].T = zp;
317: } else { /* wedge */
318: /* horizontal velocity */
319: zp = (jw-0.5)*dz; xp = iw*dx+zp*itb;
320: x[j][i].u = HorizVelocity(xp,zp,c,d);
321: /* vertical velocity */
322: zp = jw*dz; xp = (iw-0.5)*dx+zp*itb;
323: x[j][i].w = VertVelocity(xp,zp,c,d);
324: /* pressure */
325: zp = (jw-0.5)*dz; xp = (iw-0.5)*dx+zp*itb;
326: x[j][i].p = Pressure(xp,zp,c,d);
327: /* temperature */
328: x[j][i].T = 1.0;
330: }
331: }
332: }
334: /* Trash initial guess */
335: /*
336: for (j=ys; j<ys+ym; j++) {
337: for (i=xs; i<xs+xm; i++) {
338: x[j][i].u = 0.0;
339: x[j][i].w = 0.0;
340: x[j][i].p = 0.0;
341: x[j][i].T = 0.0;
342: }
343: }
344: */
346: /* Restore x to Xold */
347: DAVecRestoreArray(da,((AppCtx*)dmmg[param->mglevels-1]->user)->Xold,(void**)&x);
348:
349: /* Restrict Xold to coarser levels */
350: for (mglevel=param->mglevels-1; mglevel>0; mglevel--) {
351: MatRestrict(dmmg[mglevel]->R, ((AppCtx*)dmmg[mglevel]->user)->Xold, ((AppCtx*)dmmg[mglevel-1]->user)->Xold);
352: VecPointwiseMult(dmmg[mglevel]->Rscale,((AppCtx*)dmmg[mglevel-1]->user)->Xold,((AppCtx*)dmmg[mglevel-1]->user)->Xold);
353: }
354:
355: return 0;
356: }
358: /*---------------------------------------------------------------------*/
361: int FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
362: /*---------------------------------------------------------------------*/
363: {
364: AppCtx *user = (AppCtx*)ptr;
365: Parameter *param = user->param;
366: PetscTruth back_arc;
367: int ierr,i,j,mx,mz,sh,ish,ic,jc,iw,jw,ilim,jlim;
368: int xints,xinte,yints,yinte,ivisc,ifromm,ibound;
369: PetscScalar dudxN,dudxS,dudxW,dudxE,dudzN,dudzS,dudzE,dudzW,dudxC,dudzC;
370: PetscScalar dwdzE,dwdzW,dwdxW,dwdxE,dwdzN,dwdzS,dwdxN,dwdxS,dwdxC,dwdzC;
371: PetscScalar pN,pS,pE,pW,uE,uW,uC,uN,uS,wN,wS,wE,wW,wC;
372: PetscScalar uNE,uNW,uSE,uSW,wNE,wNW,wSE,wSW;
373: PetscScalar vE, vN, vS, vW, TE, TN, TS, TW, TC, pC;
374: PetscScalar fN,fS,fE,fW,dTdxW,dTdxE,dTdzN,dTdzS,TNE,TNW,TSE,TSW,dTdxN,dTdxS;
375: PetscScalar etaN,etaS,etaE,etaW,etaC;
376: PassiveScalar dx, dz, dxp, dzp, c, d, beta, itb, sbi, pe;
377: PassiveScalar xp, zp, cb, sb, skt, z_scale, fr, x_back_arc, u_back_arc;
378: PassiveScalar eps=0.000000001, alpha_g_on_cp_units_inverse_km=4.0e-5*9.8;
383: /*
384: Define geometric and numeric parameters
385: */
386: back_arc = param->back_arc; x_back_arc = param->x_back_arc; u_back_arc=param->u_back_arc;
387: pe = param->peclet; beta = param->slab_dip;
388: ivisc = param->ivisc; ifromm = param->ifromm; ibound = param->ibound;
389: z_scale = param->lid_depth * alpha_g_on_cp_units_inverse_km;
390: ic = param->icorner; jc = param->jcorner;
391: CalcJunk(param->kappa,param->slab_age,beta,&skt,&cb,&sb,&itb,&sbi);
392: c = beta*sb/(beta*beta-sb*sb);
393: d = (beta*cb-sb)/(beta*beta-sb*sb);
394:
395: /*
396: Define global and local grid parameters
397: */
398: mx = info->mx; mz = info->my;
399: ilim = mx-1; jlim = mz-1;
400: dx = param->scaled_width/((double)(mx-2));
401: dz = param->scaled_depth/((double)(mz-2));
402: dxp = dx; dzp = dz*sbi;
403: xints = info->xs; xinte = info->xs+info->xm;
404: yints = info->ys; yinte = info->ys+info->ym;
405: fr = (1.0-dz/dx*itb)/2.0;
407: /*
408: Stokes equation, no buoyancy terms
409: Steady state advection-diffusion equation for temperature
410: */
412: for (j=yints; j<yinte; j++) {
413: jw = j-jc+1;
414: for (i=xints; i<xinte; i++) {
415: iw = i-ic+1;
417: /* X-MOMENTUM/VELOCITY */
418: if (i<ic) { /* within the slab */
419: f[j][i].u = x[j][i].u - cb;
420: } else if (j<jc) { /* within the lid */
421: f[j][i].u = x[j][i].u - u_back_arc*LidVelocity(xp,x_back_arc,back_arc);
422: } else if (i==ilim) { /* On the INFLOW boundary */
423: if (ibound==0) { /* isoviscous analytic soln */
424: zp = (jw-0.5)*dz; xp = iw*dx+zp*itb;
425: f[j][i].u = x[j][i].u - HorizVelocity(xp,zp,c,d);
426: } else { /* normal stress = 0 boundary condition */
427: uE = x[j][i].u; uW = x[j][i-1].u; pC = x[j][i].p;
428: TC = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
429: etaC = Viscosity(TC,ivisc);
430: f[j][i].u = 2.0*etaC*( uE - uW )/dxp - pC;
431: }
432: } else if (j==jlim) { /* On the OUTFLOW boundary */
433: if (ibound==0) { /* isoviscous analytic soln */
434: zp = (jw-0.5)*dz; xp = iw*dx+zp*itb;
435: f[j][i].u = x[j][i].u - HorizVelocity(xp,zp,c,d);
436: } else { /* shear stress = 0 boundary condition */
437: uN = x[j][i].u; wE = x[j-1][i+1].w;
438: uS = x[j-1][i].u; wW = x[j-1][i].w;
439: uW = UInterp(x,i-1,j-1,fr); uE = UInterp(x,i,j-1,fr);
440: f[j][i].u = sbi*( uN - uS )/dzp
441: - itb*( uE - uW )/dxp
442: + ( wE - wW )/dxp;
443: }
444: } else { /* Mantle wedge, horizontal velocity */
446: pW = x[j][i].p; pE = x[j][i+1].p;
448: TN = param->potentialT * TInterp(x,i,j,fr) * exp( j *dz*z_scale );
449: TS = param->potentialT * TInterp(x,i,j-1,fr) * exp( (j-1.0)*dz*z_scale );
450: TE = param->potentialT * x[j][i+1].T * exp( (j-0.5)*dz*z_scale );
451: TW = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
452: etaN = Viscosity(TN,ivisc); etaS = Viscosity(TS,ivisc);
453: etaW = Viscosity(TW,ivisc); etaE = Viscosity(TE,ivisc);
454: /*if (j==jc) etaS = 1.0;*/
456: /* ------ BEGIN VAR VISCOSITY USE ONLY ------- */
457: dwdxN = etaN * ( x[j][i+1].w - x[j][i].w )/dxp;
458: dwdxS = etaS * ( x[j-1][i+1].w - x[j-1][i].w )/dxp;
459: if (i<ilim-1) { wE = WInterp(x,i+1,j-1,fr); }
460: else { wE = ( x[j][i].w + x[j-1][i].w )/2.0; }
461: wC = WInterp(x,i,j-1,fr);
462: wW = WInterp(x,i-1,j-1,fr); if (i==ic) wW = sb;
463: dwdxE = etaE * ( wE - wC )/dxp;
464: dwdxW = etaW * ( wC - wW )/dxp;
465: /* ------ END VAR VISCOSITY USE ONLY ------- */
467: /* ------ BGN ISOVISC BETA != 0 USE ONLY ------- */
468: uNE = UInterp(x,i,j,fr); uNW = UInterp(x,i-1,j,fr);
469: uSE = UInterp(x,i,j-1,fr); uSW = UInterp(x,i-1,j-1,fr);
470: if (j==jc) {
471: xp = (iw+0.5)*dx+(j-1)*dz*itb; uSE = u_back_arc*LidVelocity(xp,x_back_arc,back_arc);
472: xp = (iw-0.5)*dx+(j-1)*dz*itb; uSW = u_back_arc*LidVelocity(xp,x_back_arc,back_arc);
473: }
474: dudxN = etaN * ( uNE - uNW )/dxp; dudxS = etaS * ( uSE - uSW )/dxp;
475: dudzE = etaE * ( uNE - uSE )/dzp; dudzW = etaW * ( uNW - uSW )/dzp;
476: /* ------ END ISOVISC BETA != 0 USE ONLY ------- */
478: dudzN = etaN * ( x[j+1][i].u - x[j][i].u )/dzp;
479: dudzS = etaS * ( x[j][i].u - x[j-1][i].u )/dzp;
480: dudxE = etaE * ( x[j][i+1].u - x[j][i].u )/dxp;
481: dudxW = etaW * ( x[j][i].u - x[j][i-1].u )/dxp;
482: if (j==jc) {
483: if (ibound==0) { /* apply isoviscous boundary condition */
484: xp = iw*dx;
485: dudzS = etaS * sb*(HorizVelocity(xp,eps,c,d)-HorizVelocity(xp,-eps,c,d))/eps/2.0;
486: } else /* force u=0 on the lid-wedge interface (off-grid point) */
487: dudzS = etaS * ( 2.0*x[j][i].u - 2.0*x[j-1][i].u )/dzp;
488: }
490: f[j][i].u = -( pE - pW )/dxp /* X-MOMENTUM EQUATION*/
491: +( dudxE - dudxW )/dxp * (1.0+itb*itb)
492: +( dudzN - dudzS )/dzp * sbi*sbi
493: -( ( dudxN - dudxS )/dzp
494: +( dudzE - dudzW )/dxp ) * itb*sbi;
496: if (ivisc>0) {
497: f[j][i].u = f[j][i].u + ( dudxE - dudxW )/dxp
498: + ( dwdxN - dwdxS )/dzp * sbi
499: - ( dwdxE - dwdxW )/dxp * itb;
500: }
501: }
503: /* Z-MOMENTUM/VELOCITY */
504: if (i<ic) { /* within the slab */
505: f[j][i].w = x[j][i].w - sb;
506: } else if (j<jc) { /* within the lid */
507: f[j][i].w = x[j][i].w - 0.0;
508: } else if (j==jlim) { /* On the OUTFLOW boundary */
509: if (ibound==0) { /* isoviscous analytic soln */
510: zp = jw*dz; xp = (iw-0.5)*dx+zp*itb;
511: f[j][i].w = x[j][i].w - VertVelocity(xp,zp,c,d);
512: } else { /* normal stress = 0 boundary condition */
513: wN = x[j][i].w; wS = x[j-1][i].w; pC = x[j][i].p;
514: wW = WInterp(x,i-1,j-1,fr); if (i==ic) wW = sb;
515: if (i==ilim) wE = ( x[j][i].w + x[j-1][i].w )/2.0;
516: else wE = WInterp(x,i,j-1,fr);
517: TC = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
518: etaC = Viscosity(TC,ivisc);
519: f[j][i].w = 2.0*etaC*( sbi*( wN - wS )/dzp
520: -itb*( wE - wW )/dxp ) - pC;
521: }
522: } else if (i==ilim) { /* On the INFLOW boundary */
523: if (ibound==0) { /* isoviscous analytic soln */
524: zp = jw*dz; xp = (iw-0.5)*dx+zp*itb;
525: f[j][i].w = x[j][i].w - VertVelocity(xp,zp,c,d);
526: } else { /* shear stress = 0 boundary condition */
527: uN = x[j+1][i-1].u; wE = x[j][i].w;
528: uS = x[j][i-1].u; wW = x[j][i-1].w;
529: uW = UInterp(x,i-2,j,fr); uE = UInterp(x,i-1,j,fr);
530: f[j][i].w = sbi*( uN - uS )/dzp
531: - itb*( uE - uW )/dxp
532: + ( wE - wW )/dxp;
533: if (j==jlim-1) {
534: f[j][i].w = x[j][i].w - x[j-1][i].w;
535: }
536: }
537: } else { /* Mantle wedge, vertical velocity */
538:
539: pE = PInterp(x,i,j,fr); pW = PInterp(x,i-1,j,fr); pS = x[j][i].p; pN = x[j+1][i].p;
540: if ( (i==ic) && (ibound==0) ) { zp = jw*dz; xp = zp*itb; pW = Pressure(xp,zp,c,d); }
541:
542: TN = param->potentialT * x[j+1][i].T * exp( (j+0.5)*dz*z_scale );
543: TS = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
544: TE = param->potentialT * TInterp(x,i,j,fr) * exp( j *dz*z_scale );
545: TW = param->potentialT * TInterp(x,i-1,j,fr) * exp( j *dz*z_scale );
546: etaN = Viscosity(TN,ivisc); etaS = Viscosity(TS,ivisc);
547: etaW = Viscosity(TW,ivisc); etaE = Viscosity(TE,ivisc);
549: /* ------ BGN VAR VISCOSITY USE ONLY ------- */
550: dudzE = etaE * ( x[j+1][i].u - x[j][i].u )/dzp;
551: dudzW = etaW * ( x[j+1][i-1].u - x[j][i-1].u )/dzp;
552: uE = UInterp(x,i,j,fr); uC = UInterp(x,i-1,j,fr);
553: uW = UInterp(x,i-2,j,fr); if (i==ic) uW = 2.0*cb - uC;
554: dudxE = etaE * ( uE - uC )/dxp;
555: dudxW = etaW * ( uC - uW )/dxp;
556: /* ------ END VAR VISCOSITY USE ONLY ------- */
558: /* ------ BGN ISOVISC BETA != 0 USE ONLY ------- */
559: wNE = WInterp(x,i,j,fr); wSE = WInterp(x,i,j-1,fr);
560: wNW = WInterp(x,i-1,j,fr); wSW = WInterp(x,i-1,j-1,fr);
561: if (i==ic) { wNW = wSW = sb; }
562: dwdzE = etaE * ( wNE - wSE )/dzp; dwdzW = etaW * ( wNW - wSW )/dzp;
563: dwdxN = etaN * ( wNE - wNW )/dxp; dwdxS = etaS * ( wSE - wSW )/dxp;
564: /* ------ END ISOVISC BETA != 0 USE ONLY ------- */
566: dwdzN = etaN * ( x[j+1][i].w - x[j][i].w )/dzp;
567: dwdzS = etaS * ( x[j][i].w - x[j-1][i].w )/dzp;
568: dwdxE = etaE * ( x[j][i+1].w - x[j][i].w )/dxp;
569: dwdxW = etaW * ( x[j][i].w - x[j][i-1].w )/dxp;
570: if (i==ic) {
571: if (ibound==0) { /* apply isoviscous boundary condition */
572: zp = jw*dz; xp = itb*zp;
573: dwdxW = etaW * ( VertVelocity(xp+eps,zp,c,d) - VertVelocity(xp,zp,c,d) )/eps;
574: } else /* force w=sin(beta) on the slab-wedge interface (off-grid point) */
575: dwdxW = etaW * ( 2.0*x[j][i].w - 2.0*sb )/dxp;
576: }
578: /* Z-MOMENTUM */
579: f[j][i].w = ( pE - pW )/dxp * itb /* constant viscosity terms */
580: -( pN - pS )/dzp * sbi
581: +( dwdzN - dwdzS )/dzp * sbi*sbi
582: +( dwdxE - dwdxW )/dxp * (itb*itb+1.0)
583: -( ( dwdzE - dwdzW )/dxp
584: +( dwdxN - dwdxS )/dzp ) * itb*sbi;
586: if (ivisc>0) {
587: f[j][i].w = f[j][i].w + ( dwdxE - dwdxW )/dxp * itb*itb
588: + ( dwdzN - dwdzS )/dzp * sbi*sbi
589: -( ( dwdzE - dwdzW )/dxp
590: +( dwdxN - dwdxS )/dzp ) * itb*sbi
591: - ( dudxE - dudxW )/dxp * itb
592: + ( dudzE - dudzW )/dxp * sbi;
593: }
594: }
595:
596: /* CONTINUITY/PRESSURE */
597: if ( (j<jc) || (i<ic-1) ) { /* within slab or lid */
598: f[j][i].p = x[j][i].p - 0.0;
600: } else if ( (j==jlim)&&(i==ic-1) ) {
601: f[j][i].p = x[j][i].p - 0.0;
602:
603: } else if ( (ibound==0) && ((i==ilim)||(j==jlim)) ) { /* isoviscous inflow/outflow BC */
604: zp = (jw-0.5)*dz; xp = (iw-0.5)*dx+zp*itb;
605: f[j][i].p = x[j][i].p - Pressure(xp,zp,c,d);
607: } else if (i==ic-1) { /* just west of the slab-wedge interface constrain
608: pressure using the x-momentum equation */
610: pE = x[j][i+1].p; pW = x[j][i].p; /* pW IS THE UNKNOWN */
611:
612: TN = param->potentialT * TInterp(x,i,j,fr) * exp( j *dz*z_scale );
613: TS = param->potentialT * TInterp(x,i,j-1,fr) * exp( (j-1.0)*dz*z_scale );
614: TE = param->potentialT * x[j][i+1].T * exp( (j-0.5)*dz*z_scale );
615: TW = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
616: etaN = Viscosity(TN,ivisc); etaS = Viscosity(TS,ivisc);
617: etaW = Viscosity(TW,ivisc); etaE = Viscosity(TE,ivisc);
619: /* ------ BGN VAR VISCOSITY USE ONLY ------- */
620: dwdxN = etaN * ( x[j][i+1].w - (2*sb - x[j][i+1].w) )/dxp;
621: dwdxS = etaS * ( x[j-1][i+1].w - (2*sb - x[j-1][i+1].w) )/dxp;
622: wE = WInterp(x,i+1,j-1,fr); wC = sb; wW = 2*sb - wE;
623: dwdxE = etaE * ( wE - wC )/dxp;
624: dwdxW = etaW * ( wC - wW )/dxp;
625: /* ------ END VAR VISCOSITY USE ONLY ------- */
627: /* ------ BGN BETA != 0 USE ONLY ------- */
628: uNE = UInterp(x,i,j,fr); uNW = 2*cb - uNE;
629: uSE = UInterp(x,i,j-1,fr); uSW = 2*cb - uSE;
630: if (j==jc) uSE = 0.0;
631: dudxN = etaN * ( uNE - uNW )/dxp; dudxS = etaS * ( uSE - uSW )/dxp;
632: dudzE = etaE * ( uNE - uSE )/dzp; dudzW = 0.0;
633: if (j==jc) dudxS = 0.0;
634: /* ------ BGN BETA != 0 USE ONLY ------- */
636: dudxE = etaE * ( x[j][i+1].u - x[j][i].u )/dxp;
637: dudxW = 0.0;
639: f[j][i].p = -( pE - pW )/dxp /* X-MOMENTUM */
640: +( dudxE - dudxW )/dxp * (1.0+itb*itb)
641: -( ( dudxN - dudxS )/dzp
642: +( dudzE - dudzW )/dxp ) * itb*sbi;
644: if (ivisc>0) {
645: f[j][i].p = f[j][i].p + ( dudxE - dudxW )/dxp
646: + ( dwdxN - dwdxS )/dzp * sbi
647: - ( dwdxE - dwdxW )/dxp * itb;
648: }
649: } else { /* interior of the domain */
650: uW = x[j][i-1].u; uE = x[j][i].u;
651: wS = x[j-1][i].w; wN = x[j][i].w;
652: wW = WInterp(x,i-1,j-1,fr); if (i==ic) wW = sb;
653: if (i==ilim) wE = ( x[j][i].w + x[j-1][i].w )/2.0;
654: else wE = WInterp(x,i,j-1,fr);
656: f[j][i].p = ( uE - uW )/dxp /* CONTINUITY ON PRESSURE POINTS */
657: -( wE - wW )/dxp * itb
658: +( wN - wS )/dzp * sbi;
660: if ( (ivisc==10)&&(i<ilim)&&(j<jlim) ) {
661: wW = x[j][i].w; wN = WInterp(x,i,j,fr);
662: wE = x[j][i+1].w; wS = WInterp(x,i,j-1,fr);
663: uW = UInterp(x,i-1,j,fr); uE = UInterp(x,i,j,fr);
664:
665: f[j][i].p = f[j][i].p +
666: ( uE - uW )/dxp - /* CONTINUITY ON VACANT POINTS */
667: ( wE - wW )/dxp * itb +
668: ( wN - wS )/dzp * sbi;
669: }
670:
671: }
672:
673: /* TEMPERATURE EQUATION */
674: zp = (j-0.5)*dz; xp = (iw-0.5)*dx+zp*itb;
675: if (i<=1) { /* dirichlet on boundary along slab side */
676: f[j][i].T = x[j][i].T - 1.0;
677: } else if ( (j<=2) && (i<ic) ) { /* slab inflow dirichlet */
678: f[j][i].T = x[j][i].T - erf(-xp*sb*param->lid_depth/2.0/skt);
679: } else if (j==0) { /* force T=0 on surface */
680: f[j][i].T = x[j][i].T + x[j+1][i].T;
681: } else if (j>=jlim-1) { /* neumann on outflow boundary */
682: if (x[j][i].w<0.0)
683: f[j][i].T = x[j][i].T - 1.0;
684: else
685: f[j][i].T = x[j][i].T - x[j-1][i].T;
686: } else if (i>=ilim-1) /* (xp>=20.0) */ {
687: if (back_arc && (x[j][i].u>0))
688: f[j][i].T = x[j][i].T - x[j][i-1].T;
689: else
690: f[j][i].T = x[j][i].T - PetscMin(zp,1.0);
691: } else { /* interior of the domain */
693: uW = x[j][i-1].u; uE = x[j][i].u;
694: wS = x[j-1][i].w; wN = x[j][i].w;
695: wE = WInterp(x,i,j-1,fr); wW = WInterp(x,i-1,j-1,fr);
697: if (i==ic) { /* Just east of slab-wedge interface */
698: if (j<jc) {
699: wW = wE = wN = wS = 0.0;
700: uW = uE = 0.0;
701: } else {
702: wW = sb;
703: }
704: }
706: if (i==ic-1) wE = sb; /* Just west of slab-wedge interface */
708: TNE = TInterp(x,i,j,fr); TNW = TInterp(x,i-1,j,fr);
709: TSE = TInterp(x,i,j-1,fr); TSW = TInterp(x,i-1,j-1,fr);
710: dTdxN = ( TNE - TNW )/dxp;
711: dTdxS = ( TSE - TSW )/dxp;
712: dTdzN = ( x[j+1][i].T - x[j][i].T )/dzp;
713: dTdzS = ( x[j][i].T - x[j-1][i].T )/dzp;
714: dTdxE = ( x[j][i+1].T - x[j][i].T )/dxp;
715: dTdxW = ( x[j][i].T - x[j][i-1].T )/dxp;
716:
717: f[j][i].T = ( ( dTdzN - dTdzS )/dzp * sbi*sbi + /* diffusion term */
718: ( dTdxE - dTdxW )/dxp * (1.0+itb*itb) -
719: ( dTdxN - dTdxS )/dzp * 2.0*itb*sbi )*dx*dz/pe;
721: if ( (j<jc) && (i>=ic) ) { /* don't advect in lid */
722: fE = fW = fN = fS = 0.0;
724: } else if ( (ifromm==0) ||(i>=ilim-2)||(j>=jlim-2)||(i<=2) ) { /* finite volume advection */
725: TN = ( x[j][i].T + x[j+1][i].T )/2.0; TS = ( x[j][i].T + x[j-1][i].T )/2.0;
726: TE = ( x[j][i].T + x[j][i+1].T )/2.0; TW = ( x[j][i].T + x[j][i-1].T )/2.0;
727: fN = wN*TN*dxp; fS = wS*TS*dxp;
728: fE = ( uE*sb - wE*cb )*TE*dzp;
729: fW = ( uW*sb - wW*cb )*TW*dzp;
730:
731: } else { /* Fromm advection scheme */
732: vN = wN; vS = wS; vE = uE*sb - wE*cb; vW = uW*sb - wW*cb;
733: fE = ( vE *(-x[j][i+2].T + 5.0*(x[j][i+1].T+x[j][i].T)-x[j][i-1].T)/8.0
734: - fabs(vE)*(-x[j][i+2].T + 3.0*(x[j][i+1].T-x[j][i].T)+x[j][i-1].T)/8.0 )*dzp;
735: fW = ( vW *(-x[j][i+1].T + 5.0*(x[j][i].T+x[j][i-1].T)-x[j][i-2].T)/8.0
736: - fabs(vW)*(-x[j][i+1].T + 3.0*(x[j][i].T-x[j][i-1].T)+x[j][i-2].T)/8.0 )*dzp;
737: fN = ( vN *(-x[j+2][i].T + 5.0*(x[j+1][i].T+x[j][i].T)-x[j-1][i].T)/8.0
738: - fabs(vN)*(-x[j+2][i].T + 3.0*(x[j+1][i].T-x[j][i].T)+x[j-1][i].T)/8.0 )*dxp;
739: fS = ( vS *(-x[j+1][i].T + 5.0*(x[j][i].T+x[j-1][i].T)-x[j-2][i].T)/8.0
740: - fabs(vS)*(-x[j+1][i].T + 3.0*(x[j][i].T-x[j-1][i].T)+x[j-2][i].T)/8.0 )*dxp;
741: }
742:
743: f[j][i].T = f[j][i].T -
744: ( fE - fW + fN - fS );
745: }
746: }
747: }
748:
749: return(0);
750: }
752: /* ------------------------------------------------------------------- */
755: /*
756: FormInitialGuess - Forms initial approximation.
758: Input Parameters:
759: user - user-defined application context
760: X - vector
762: Output Parameter:
763: X - vector
764: */
765: /* ------------------------------------------------------------------- */
766: int FormInitialGuess(SNES snes,Vec X,void *ptr)
767: {
768: DMMG dmmg = (DMMG)ptr;
769: AppCtx *user = (AppCtx*)dmmg->user;
770: int ierr;
772: VecCopy(user->Xold, X);
774: /* calculate the residual on fine mesh, but only the first time this is called */
775: if (user->fnorm_ini == 0.0) {
776: SNESComputeFunction(snes,user->Xold,user->func);
777: VecNorm(user->func,NORM_2,&user->fnorm_ini);
778: }
779:
780: return 0;
781: }
783: /*--------------------------UTILITY FUNCTION BELOW THIS LINE-----------------------------*/
785: /*---------------------------------------------------------------------*/
788: int SetParams(Parameter *param, int mx, int mz)
789: /*---------------------------------------------------------------------*/
790: {
792: PetscReal SEC_PER_YR = 3600.00*24.00*356.2500;
794: /* domain geometry */
795: param->icorner = (int)(mx/6+1); /* gridpoints */
796: param->jcorner = (int)((mz-2)/12+1);/* gridpoints */
797: param->width = 1200.0; /* km */
798: param->depth = 600.0; /* km */
799: param->slab_dip = HALFPI; /* 90 degrees */
800: PetscOptionsGetInt(PETSC_NULL, "-icorner",&(param->icorner),PETSC_NULL);
801: PetscOptionsGetInt(PETSC_NULL, "-jcorner",&(param->jcorner),PETSC_NULL);
802: PetscOptionsGetReal(PETSC_NULL,"-width",&(param->width),PETSC_NULL);
803: PetscOptionsGetReal(PETSC_NULL,"-depth",&(param->depth),PETSC_NULL);
804: PetscOptionsGetReal(PETSC_NULL,"-slab_dip",&(param->slab_dip),PETSC_NULL);
806: /* back-arc */
807: param->back_arc = PETSC_FALSE; /* no back arc spreading */
808: param->x_back_arc = 0.0; /* km */
809: param->u_back_arc = 1.0; /* full spreading at velocity of slab */
810: PetscOptionsHasName(PETSC_NULL,"-back_arc",&(param->back_arc));
811: PetscOptionsGetReal(PETSC_NULL,"-back_arc",&(param->x_back_arc),PETSC_NULL);
812: PetscOptionsGetReal(PETSC_NULL,"-back_arc_velocity",&(param->u_back_arc),PETSC_NULL);
813: if (param->back_arc) {
814: PetscPrintf(PETSC_COMM_WORLD,"Dist to back arc = %g km, ",param->x_back_arc);
815: PetscPrintf(PETSC_COMM_WORLD,"Full spreading rate of back arc (scaled) = %g \n",param->u_back_arc);
816: }
818: /* physics parameters */
819: param->slab_velocity = 5.0; /* cm/yr */
820: param->slab_age = 50.0; /* Ma */
821: param->kappa = 0.7272e-6; /* m^2/sec */
822: param->potentialT = 1300.0; /* degrees C */
823: param->ivisc = 0; /* 0=constant, 1=diffusion creep, 2=simple T dependent */
824: PetscOptionsGetReal(PETSC_NULL,"-slab_velocity",&(param->slab_velocity),PETSC_NULL);
825: PetscOptionsGetReal(PETSC_NULL,"-slab_age",&(param->slab_age),PETSC_NULL);
826: PetscOptionsGetReal(PETSC_NULL,"-kappa",&(param->kappa),PETSC_NULL);
827: PetscOptionsGetReal(PETSC_NULL,"-potentialT",&(param->potentialT),PETSC_NULL);
828: PetscOptionsGetInt(PETSC_NULL, "-ivisc",&(param->ivisc),PETSC_NULL);
830: /* boundaries */
831: param->ibound = param->ivisc; /* 0=isovisc analytic, 1,2,...= stress free */
832: PetscOptionsGetInt(PETSC_NULL,"-ibound",&(param->ibound),PETSC_NULL);
834: /* misc */
835: param->ifromm = 1; /* advection scheme: 0=finite vol, 1=Fromm */
836: PetscOptionsGetInt(PETSC_NULL,"-ifromm",&(param->ifromm),PETSC_NULL);
838: /* unit conversions and derived parameters */
839: param->lid_depth = (param->jcorner - 1) * param->depth/((double)(mz-2));
840: PetscOptionsGetReal(PETSC_NULL,"-lid_depth",&(param->lid_depth),PETSC_NULL);
841: param->slab_dip = param->slab_dip*HALFPI/90.0;
842: param->scaled_width = param->width/param->lid_depth;
843: param->scaled_depth = param->depth/param->lid_depth;
844: param->x_back_arc = param->x_back_arc/param->lid_depth;
845: param->peclet = param->slab_velocity/100.0/SEC_PER_YR /* m/sec */
846: * param->lid_depth*1000.0 /* m */
847: / param->kappa; /* m^2/sec */
848: PetscOptionsGetReal(PETSC_NULL,"-peclet",&(param->peclet),PETSC_NULL);
849: PetscPrintf(PETSC_COMM_WORLD,"Lid depth = %g km, ",param->lid_depth);
850: PetscPrintf(PETSC_COMM_WORLD,"Peclet number = %g\n",param->peclet);
851: if ( (param->ibound==0) && (param->ivisc>0) )
852: PetscPrintf(PETSC_COMM_WORLD,"Warning: isoviscous BC may be inconsistent w/ var viscosity!!\n");
854: return 0;
855: }
858: /*---------------------------------------------------------------------*/
861: PetscScalar Viscosity(PetscScalar T, int iVisc)
862: /*---------------------------------------------------------------------*/
863: {
864: PetscScalar result;
865: double p1 = 1.32047792e-12, p2 = 335.0e3/8.314510;
866: double t1 = 1200.0;
867: /*
868: p1 = exp( -p2/(1473 K) ) so the cutoff for high viscosity
869: occurs at 1200 C. Below this cutoff, all eta( T<1200 ) = 1.0;
870: p2=Ea/R is from van Keken's subroutine
871: */
873: if (iVisc==0) { /* constant viscosity */
874: result = 1.0;
875: } else if (iVisc==1) { /* diffusion creep rheology */
876: result = p1*PetscExpScalar(p2/(T+273.0));
877: } else if (iVisc==2) { /* ad hoc T-dependent rheology */
878: result = t1/T;
879: } else if (iVisc==3) {
880: result = 1.0;
881: }
883: if (result<1.0)
884: return result;
885: else
886: return 1.0;
887: }
889: /*---------------------------------------------------------------------*/
892: PassiveScalar LidVelocity(PassiveScalar x, PassiveScalar xBA, PetscTruth BA)
893: /*---------------------------------------------------------------------*/
894: {
895: PassiveScalar localize = 10.0;
897: if (BA)
898: return ( 1.0 + tanh( (x - xBA)*localize ) )/2.0;
899: else
900: return 0.0;
901: }
903: /*---------------------------------------------------------------------*/
906: PetscScalar UInterp(Field **x, int i, int j, PetscScalar fr)
907: /*---------------------------------------------------------------------*/
908: {
909: PetscScalar p,m;
910: p = (1.0-fr)*x[j+1][i].u + fr*x[j+1][i+1].u;
911: m = (1.0-fr)*x[j][i+1].u + fr*x[j][i].u;
912: return (p + m)/2.0;
913: }
915: /*---------------------------------------------------------------------*/
918: PetscScalar WInterp(Field **x, int i, int j, PetscScalar fr)
919: /*---------------------------------------------------------------------*/
920: {
921: PetscScalar p,m;
922: p = (1.0-fr)*x[j+1][i].w + fr*x[j+1][i+1].w;
923: m = (1.0-fr)*x[j][i+1].w + fr*x[j][i].w;
924: return (p + m)/2.0;
925: }
927: /*---------------------------------------------------------------------*/
930: PetscScalar PInterp(Field **x, int i, int j, PetscScalar fr)
931: /*---------------------------------------------------------------------*/
932: {
933: PetscScalar p,m;
934: p = (1.0-fr)*x[j+1][i].p + fr*x[j+1][i+1].p;
935: m = (1.0-fr)*x[j][i+1].p + fr*x[j][i].p;
936: return (p + m)/2.0;
937: }
939: /*---------------------------------------------------------------------*/
942: PetscScalar TInterp(Field **x, int i, int j, PetscScalar fr)
943: /*---------------------------------------------------------------------*/
944: {
945: PetscScalar p,m;
946: p = (1.0-fr)*x[j+1][i].T + fr*x[j+1][i+1].T;
947: m = (1.0-fr)*x[j][i+1].T + fr*x[j][i].T;
948: return (p + m)/2.0;
949: }
951: /*---------------------------------------------------------------------*/
954: PetscScalar HorizVelocity(PetscScalar x, PetscScalar z, PetscScalar c, PetscScalar d)
955: /*---------------------------------------------------------------------*/
956: {
957: PetscScalar r, st, ct, th;
958: r = sqrt(x*x+z*z);
959: st = z/r; ct = x/r; th = atan(z/x);
960: return ct*(c*th*st+d*(st+th*ct)) + st*(c*(st-th*ct)+d*th*st);
961: }
963: /*---------------------------------------------------------------------*/
966: PetscScalar VertVelocity(PetscScalar x, PetscScalar z, PetscScalar c, PetscScalar d)
967: /*---------------------------------------------------------------------*/
968: {
969: PetscScalar r, st, ct, th;
970: r = sqrt(x*x+z*z);
971: st = z/r; ct = x/r; th = atan(z/x);
972: return st*(c*th*st+d*(st+th*ct)) - ct*(c*(st-th*ct)+d*th*st);
973: }
975: /*---------------------------------------------------------------------*/
978: PetscScalar Pressure(PetscScalar x, PetscScalar z, PetscScalar c, PetscScalar d)
979: /*---------------------------------------------------------------------*/
980: {
981: PetscScalar r, st, ct;
982: r = sqrt(x*x+z*z);
983: st = z/r; ct = x/r;
984: return (-2.0*(c*ct-d*st)/r);
985: }
987: /*---------------------------------------------------------------------*/
990: void CalcJunk(double kappa,double slab_age,double beta,double *skt,double *cb,double *sb, double *itb,
991: double *sbi)
992: /*---------------------------------------------------------------------*/
993: {
994: PetscReal SEC_PER_YR = 3600.00*24.00*356.2500;
996: *skt = sqrt(kappa*slab_age*SEC_PER_YR);
997: *cb = cos(beta); *sb = sin(beta);
998: *itb = 1.0/tan(beta); *sbi = 1.0/(*sb);
999: }