← Back to branch summary

~ubuntu-branches/ubuntu/warty/petsc/warty

« back to all changes in this revision

Viewing changes to src/snes/examples/tutorials/ex18.c.html

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2004-06-07 13:41:43 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20040607134143-92p586zrauvie0le
Tags: 2.2.0-2
http://bugs.debian.org/249617
* Upstream patch level 2.
* New PETSC_BOPT_EXTRA option for different BOPT and lib names, with _c++
  symlinks only for plain and single (closes: #249617).
* New DEBIAN_DIST=contrib option to link with hypre, parmetis (closes:
  #249619).
* Combined petsc-c and petsc-fortran substvars into petsc-compilers.
* Extra quote in -dev prerm eliminates "too many arguments" problem.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/snes/examples/tutorials/ex18.c.html
3
3
<html>
4
4
<head>
5
5
<title></title>
6
 
<meta name="generator" content="c2html 0.9.1">
7
 
<meta name="date" content="2002-05-31T16:32:10+00:00">
 
6
<meta name="generator" content="c2html 0.9.4">
 
7
<meta name="date" content="2004-02-27T20:03:07+00:00">
8
8
</head>
9
9
 
10
10
<body bgcolor="#FFFFFF">
11
11
<pre width="80"><a name="line1">  1: </a><font color="#B22222">/* $Id: ex18.c,v 1.23 2001/08/07 21:31:17 bsmith Exp $ */</font>
12
12
 
13
13
 
14
 
<a name="line4">  4: </a>static char help[] =<font color="#666666">"Nonlinear Radiative Transport PDE with multigrid in 2d.n</font>
15
 
<a name="line5">  5: </a><font color="#666666">Uses 2-dimensional distributed arrays.n</font>
16
 
<a name="line6">  6: </a><font color="#666666">A 2-dim simplified Radiative Transport test problem is used, with analytic Jacobian. n</font>
17
 
<a name="line7">  7: </a><font color="#666666">n</font>
18
 
<a name="line8">  8: </a><font color="#666666">  Solves the linear systems via multilevel methods n</font>
19
 
<a name="line9">  9: </a><font color="#666666">n</font>
20
 
<a name="line10"> 10: </a><font color="#666666">The command linen</font>
21
 
<a name="line11"> 11: </a><font color="#666666">options are:n</font>
22
 
<a name="line12"> 12: </a><font color="#666666">  -tleft &lt;tl&gt;, where &lt;tl&gt; indicates the left Diriclet BC n</font>
23
 
<a name="line13"> 13: </a><font color="#666666">  -tright &lt;tr&gt;, where &lt;tr&gt; indicates the right Diriclet BC n</font>
24
 
<a name="line14"> 14: </a><font color="#666666">  -beta &lt;beta&gt;, where &lt;beta&gt; indicates the exponent in T nn"</font>;
 
14
<a name="line4">  4: </a>static char help[] =<font color="#666666">"Nonlinear Radiative Transport PDE with multigrid in 2d.\n\</font>
 
15
<a name="line5">  5: </a><font color="#666666">Uses 2-dimensional distributed arrays.\n\</font>
 
16
<a name="line6">  6: </a><font color="#666666">A 2-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \n\</font>
 
17
<a name="line7">  7: </a><font color="#666666">\n\</font>
 
18
<a name="line8">  8: </a><font color="#666666">  Solves the linear systems via multilevel methods \n\</font>
 
19
<a name="line9">  9: </a><font color="#666666">\n\</font>
 
20
<a name="line10"> 10: </a><font color="#666666">The command line\n\</font>
 
21
<a name="line11"> 11: </a><font color="#666666">options are:\n\</font>
 
22
<a name="line12"> 12: </a><font color="#666666">  -tleft &lt;tl&gt;, where &lt;tl&gt; indicates the left Diriclet BC \n\</font>
 
23
<a name="line13"> 13: </a><font color="#666666">  -tright &lt;tr&gt;, where &lt;tr&gt; indicates the right Diriclet BC \n\</font>
 
24
<a name="line14"> 14: </a><font color="#666666">  -beta &lt;beta&gt;, where &lt;beta&gt; indicates the exponent in T \n\n"</font>;
25
25
 
26
26
<a name="line16"> 16: </a><font color="#B22222">/*T</font>
27
27
<a name="line17"> 17: </a><font color="#B22222">   Concepts: <A href="../../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A>^solving a system of nonlinear equations</font>
58
58
<a name="line48"> 48: </a><font color="#B22222">/* User-defined application context */</font>
59
59
 
60
60
<a name="line50"> 50: </a><font color="#4169E1">typedef</font> <font color="#4169E1">struct</font> {
61
 
<a name="line51"> 51: </a>   PetscReal  tleft,tright;  <font color="#B22222">/* Dirichlet boundary conditions */</font>
62
 
<a name="line52"> 52: </a>   PetscReal  beta,bm1,coef; <font color="#B22222">/* nonlinear diffusivity parameterizations */</font>
 
61
<a name="line51"> 51: </a>   <A href="../../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A>  tleft,tright;  <font color="#B22222">/* Dirichlet boundary conditions */</font>
 
62
<a name="line52"> 52: </a>   <A href="../../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A>  beta,bm1,coef; <font color="#B22222">/* nonlinear diffusivity parameterizations */</font>
63
63
<a name="line53"> 53: </a>} AppCtx;
64
64
 
65
65
<a name="line55"> 55: </a><strong><font color="#228B22">#define POWFLOP 5 </font><font color="#B22222">/* assume a pow() takes five flops */</font><font color="#228B22"></font></strong>
68
68
<a name="line58"> 58: </a><strong><font color="#4169E1">extern int FormFunction(<A href="../../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A>,<A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A>,<A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A>,void*)</font></strong>;
69
69
<a name="line59"> 59: </a><strong><font color="#4169E1">extern int FormJacobian(<A href="../../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A>,<A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A>,<A href="../../../../docs/manualpages/Mat/Mat.html#Mat">Mat</A>*,<A href="../../../../docs/manualpages/Mat/Mat.html#Mat">Mat</A>*,<A href="../../../../docs/manualpages/Mat/MatStructure.html#MatStructure">MatStructure</A>*,void*)</font></strong>;
70
70
 
71
 
<a name="line61"> 61: </a><strong><font color="#4169E1"><a name="main"></a>int main(int argc,char **argv)</font></strong>
72
 
<a name="line62"> 62: </a>{
73
 
<a name="line63"> 63: </a>  <A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A>      *dmmg;
74
 
<a name="line64"> 64: </a>  <A href="../../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A>      snes;
75
 
<a name="line65"> 65: </a>  AppCtx    user;
76
 
<a name="line66"> 66: </a>  int       ierr,its,lits;
77
 
<a name="line67"> 67: </a>  PetscReal litspit;
78
 
<a name="line68"> 68: </a>  <A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>        da;
79
 
 
80
 
<a name="line70"> 70: </a>  <A href="../../../../docs/manualpages/Sys/PetscInitialize.html#PetscInitialize">PetscInitialize</A>(&amp;argc,&amp;argv,PETSC_NULL,help);
81
 
 
82
 
<a name="line72"> 72: </a>  <font color="#B22222">/* set problem parameters */</font>
83
 
<a name="line73"> 73: </a>  user.tleft  = 1.0;
84
 
<a name="line74"> 74: </a>  user.tright = 0.1;
85
 
<a name="line75"> 75: </a>  user.beta   = 2.5;
86
 
<a name="line76"> 76: </a>  <A href="../../../../docs/manualpages/Sys/PetscOptionsGetReal.html#PetscOptionsGetReal">PetscOptionsGetReal</A>(PETSC_NULL,<font color="#666666">"-tleft"</font>,&amp;user.tleft,PETSC_NULL);
87
 
<a name="line77"> 77: </a>  <A href="../../../../docs/manualpages/Sys/PetscOptionsGetReal.html#PetscOptionsGetReal">PetscOptionsGetReal</A>(PETSC_NULL,<font color="#666666">"-tright"</font>,&amp;user.tright,PETSC_NULL);
88
 
<a name="line78"> 78: </a>  <A href="../../../../docs/manualpages/Sys/PetscOptionsGetReal.html#PetscOptionsGetReal">PetscOptionsGetReal</A>(PETSC_NULL,<font color="#666666">"-beta"</font>,&amp;user.beta,PETSC_NULL);
89
 
<a name="line79"> 79: </a>  user.bm1  = user.beta - 1.0;
90
 
<a name="line80"> 80: </a>  user.coef = user.beta/2.0;
91
 
 
92
 
 
93
 
<a name="line83"> 83: </a>  <font color="#B22222">/*</font>
94
 
<a name="line84"> 84: </a><font color="#B22222">      Create the multilevel <A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A> data structure </font>
95
 
<a name="line85"> 85: </a><font color="#B22222">  */</font>
96
 
<a name="line86"> 86: </a>  <A href="../../../../docs/manualpages/SLES/DMMGCreate.html#DMMGCreate">DMMGCreate</A>(PETSC_COMM_WORLD,3,&amp;user,&amp;dmmg);
97
 
 
98
 
<a name="line88"> 88: </a>  <font color="#B22222">/*</font>
99
 
<a name="line89"> 89: </a><font color="#B22222">      Set the <A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A> (grid structure) for the grids.</font>
100
 
<a name="line90"> 90: </a><font color="#B22222">  */</font>
101
 
<a name="line91"> 91: </a>  <A href="../../../../docs/manualpages/DA/DACreate2d.html#DACreate2d">DACreate2d</A>(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,5,5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&amp;da);
102
 
<a name="line92"> 92: </a>  <A href="../../../../docs/manualpages/SLES/DMMGSetDM.html#DMMGSetDM">DMMGSetDM</A>(dmmg,(<A href="../../../../docs/manualpages/DA/DM.html#DM">DM</A>)da);
103
 
<a name="line93"> 93: </a>  <A href="../../../../docs/manualpages/DA/DADestroy.html#DADestroy">DADestroy</A>(da);
104
 
 
105
 
<a name="line95"> 95: </a>  <font color="#B22222">/*</font>
106
 
<a name="line96"> 96: </a><font color="#B22222">     Create the nonlinear solver, and tell the <A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A> structure to use it</font>
107
 
<a name="line97"> 97: </a><font color="#B22222">  */</font>
108
 
<a name="line98"> 98: </a>  <A href="../../../../docs/manualpages/SNES/DMMGSetSNES.html#DMMGSetSNES">DMMGSetSNES</A>(dmmg,FormFunction,FormJacobian);
109
 
 
110
 
<a name="line100">100: </a>  <font color="#B22222">/*</font>
111
 
<a name="line101">101: </a><font color="#B22222">      <A href="../../../../docs/manualpages/Profiling/PreLoadBegin.html#PreLoadBegin">PreLoadBegin</A>() means that the following section of code is run twice. The first time</font>
112
 
<a name="line102">102: </a><font color="#B22222">     through the flag PreLoading is on this the nonlinear solver is only run for a single step.</font>
113
 
<a name="line103">103: </a><font color="#B22222">     The second time through (the actually timed code) the maximum iterations is set to 10</font>
114
 
<a name="line104">104: </a><font color="#B22222">     Preload of the executable is done to eliminate from the timing the time spent bring the </font>
115
 
<a name="line105">105: </a><font color="#B22222">     executable into memory from disk (paging in).</font>
116
 
<a name="line106">106: </a><font color="#B22222">  */</font>
117
 
<a name="line107">107: </a>  <A href="../../../../docs/manualpages/Profiling/PreLoadBegin.html#PreLoadBegin">PreLoadBegin</A>(PETSC_TRUE,<font color="#666666">"Solve"</font>);
118
 
<a name="line108">108: </a>    <A href="../../../../docs/manualpages/SNES/DMMGSetInitialGuess.html#DMMGSetInitialGuess">DMMGSetInitialGuess</A>(dmmg,FormInitialGuess);
119
 
<a name="line109">109: </a>    <A href="../../../../docs/manualpages/SLES/DMMGSolve.html#DMMGSolve">DMMGSolve</A>(dmmg);
120
 
<a name="line110">110: </a>  <A href="../../../../docs/manualpages/Profiling/PreLoadEnd.html#PreLoadEnd">PreLoadEnd</A>();
121
 
<a name="line111">111: </a>  snes = DMMGGetSNES(dmmg);
122
 
<a name="line112">112: </a>  <A href="../../../../docs/manualpages/SNES/SNESGetIterationNumber.html#SNESGetIterationNumber">SNESGetIterationNumber</A>(snes,&amp;its);
123
 
<a name="line113">113: </a>  <A href="../../../../docs/manualpages/SNES/SNESGetNumberLinearIterations.html#SNESGetNumberLinearIterations">SNESGetNumberLinearIterations</A>(snes,&amp;lits);
124
 
<a name="line114">114: </a>  litspit = ((PetscReal)lits)/((PetscReal)its);
125
 
<a name="line115">115: </a>  <A href="../../../../docs/manualpages/Sys/PetscPrintf.html#PetscPrintf">PetscPrintf</A>(PETSC_COMM_WORLD,<font color="#666666">"Number of Newton iterations = %dn"</font>,its);
126
 
<a name="line116">116: </a>  <A href="../../../../docs/manualpages/Sys/PetscPrintf.html#PetscPrintf">PetscPrintf</A>(PETSC_COMM_WORLD,<font color="#666666">"Number of Linear iterations = %dn"</font>,lits);
127
 
<a name="line117">117: </a>  <A href="../../../../docs/manualpages/Sys/PetscPrintf.html#PetscPrintf">PetscPrintf</A>(PETSC_COMM_WORLD,<font color="#666666">"Average Linear its / Newton = %en"</font>,litspit);
128
 
 
129
 
<a name="line119">119: </a>  <A href="../../../../docs/manualpages/SLES/DMMGDestroy.html#DMMGDestroy">DMMGDestroy</A>(dmmg);
130
 
<a name="line120">120: </a>  <A href="../../../../docs/manualpages/Sys/PetscFinalize.html#PetscFinalize">PetscFinalize</A>();
131
 
 
132
 
<a name="line122">122: </a>  <font color="#4169E1">return</font> 0;
133
 
<a name="line123">123: </a>}
134
 
<a name="line124">124: </a><font color="#B22222">/* --------------------  Form initial approximation ----------------- */</font>
135
 
<a name="line125">125: </a><strong><font color="#4169E1"><a name="FormInitialGuess"></a>int FormInitialGuess(<A href="../../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A> snes,<A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A> X,void *ptr)</font></strong>
136
 
<a name="line126">126: </a>{
137
 
<a name="line127">127: </a>  <A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A>        dmmg = (<A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A>)ptr;
138
 
<a name="line128">128: </a>  AppCtx      *user = (AppCtx*)dmmg-&gt;user;
139
 
<a name="line129">129: </a>  int         i,j,ierr,xs,ys,xm,ym;
140
 
<a name="line130">130: </a>  PetscReal   tleft = user-&gt;tleft;
141
 
<a name="line131">131: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A> **x;
142
 
 
143
 
 
144
 
<a name="line135">135: </a>  <font color="#B22222">/* Get ghost points */</font>
145
 
<a name="line136">136: </a>  <A href="../../../../docs/manualpages/DA/DAGetCorners.html#DAGetCorners">DAGetCorners</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,&amp;xs,&amp;ys,0,&amp;xm,&amp;ym,0);
146
 
<a name="line137">137: </a>  <A href="../../../../docs/manualpages/DA/DAVecGetArray.html#DAVecGetArray">DAVecGetArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,X,(void**)&amp;x);
147
 
 
148
 
<a name="line139">139: </a>  <font color="#B22222">/* Compute initial guess */</font>
149
 
<a name="line140">140: </a>  <font color="#4169E1">for</font> (j=ys; j&lt;ys+ym; j++) {
150
 
<a name="line141">141: </a>    <font color="#4169E1">for</font> (i=xs; i&lt;xs+xm; i++) {
151
 
<a name="line142">142: </a>      x[j][i] = tleft;
152
 
<a name="line143">143: </a>    }
153
 
<a name="line144">144: </a>  }
154
 
<a name="line145">145: </a>  <A href="../../../../docs/manualpages/DA/DAVecRestoreArray.html#DAVecRestoreArray">DAVecRestoreArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,X,(void**)&amp;x);
155
 
<a name="line146">146: </a>  <font color="#4169E1">return</font>(0);
156
 
<a name="line147">147: </a>}
157
 
<a name="line148">148: </a><font color="#B22222">/* --------------------  Evaluate Function F(x) --------------------- */</font>
158
 
<a name="line149">149: </a><strong><font color="#4169E1"><a name="FormFunction"></a>int FormFunction(<A href="../../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A> snes,<A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A> X,<A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A> F,void* ptr)</font></strong>
159
 
<a name="line150">150: </a>{
160
 
<a name="line151">151: </a>  <A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A>    dmmg = (<A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A>)ptr;
161
 
<a name="line152">152: </a>  AppCtx  *user = (AppCtx*)dmmg-&gt;user;
162
 
<a name="line153">153: </a>  int     ierr,i,j,mx,my,xs,ys,xm,ym;
163
 
<a name="line154">154: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  zero = 0.0,one = 1.0;
164
 
<a name="line155">155: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  hx,hy,hxdhy,hydhx;
165
 
<a name="line156">156: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw,fn = 0.0,fs = 0.0,fe =0.0,fw = 0.0;
166
 
<a name="line157">157: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  tleft,tright,beta;
167
 
<a name="line158">158: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  **x,**f;
168
 
<a name="line159">159: </a>  <A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A>     localX;
169
 
 
170
 
<a name="line162">162: </a>  <A href="../../../../docs/manualpages/DA/DAGetLocalVector.html#DAGetLocalVector">DAGetLocalVector</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,&amp;localX);
171
 
<a name="line163">163: </a>  <A href="../../../../docs/manualpages/DA/DAGetInfo.html#DAGetInfo">DAGetInfo</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,PETSC_NULL,&amp;mx,&amp;my,0,0,0,0,0,0,0,0);
172
 
<a name="line164">164: </a>  hx    = one/(PetscReal)(mx-1);  hy    = one/(PetscReal)(my-1);
173
 
<a name="line165">165: </a>  hxdhy = hx/hy;               hydhx = hy/hx;
174
 
<a name="line166">166: </a>  tleft = user-&gt;tleft;         tright = user-&gt;tright;
175
 
<a name="line167">167: </a>  beta  = user-&gt;beta;
176
 
<a name="line168">168: </a>
177
 
<a name="line169">169: </a>  <font color="#B22222">/* Get ghost points */</font>
178
 
<a name="line170">170: </a>  <A href="../../../../docs/manualpages/DA/DAGlobalToLocalBegin.html#DAGlobalToLocalBegin">DAGlobalToLocalBegin</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,X,INSERT_VALUES,localX);
179
 
<a name="line171">171: </a>  <A href="../../../../docs/manualpages/DA/DAGlobalToLocalEnd.html#DAGlobalToLocalEnd">DAGlobalToLocalEnd</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,X,INSERT_VALUES,localX);
180
 
<a name="line172">172: </a>  <A href="../../../../docs/manualpages/DA/DAGetCorners.html#DAGetCorners">DAGetCorners</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,&amp;xs,&amp;ys,0,&amp;xm,&amp;ym,0);
181
 
<a name="line173">173: </a>  <A href="../../../../docs/manualpages/DA/DAVecGetArray.html#DAVecGetArray">DAVecGetArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,localX,(void**)&amp;x);
182
 
<a name="line174">174: </a>  <A href="../../../../docs/manualpages/DA/DAVecGetArray.html#DAVecGetArray">DAVecGetArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,F,(void**)&amp;f);
183
 
 
184
 
<a name="line176">176: </a>  <font color="#B22222">/* Evaluate function */</font>
185
 
<a name="line177">177: </a>  <font color="#4169E1">for</font> (j=ys; j&lt;ys+ym; j++) {
186
 
<a name="line178">178: </a>    <font color="#4169E1">for</font> (i=xs; i&lt;xs+xm; i++) {
187
 
<a name="line179">179: </a>      t0 = x[j][i];
188
 
 
189
 
<a name="line181">181: </a>      <font color="#4169E1">if</font> (i &gt; 0 &amp;&amp; i &lt; mx-1 &amp;&amp; j &gt; 0 &amp;&amp; j &lt; my-1) {
190
 
 
191
 
<a name="line183">183: </a>        <font color="#B22222">/* general interior volume */</font>
192
 
 
193
 
<a name="line185">185: </a>        tw = x[j][i-1];
194
 
<a name="line186">186: </a>        aw = 0.5*(t0 + tw);
195
 
<a name="line187">187: </a>        dw = PetscPowScalar(aw,beta);
196
 
<a name="line188">188: </a>        fw = dw*(t0 - tw);
197
 
 
198
 
<a name="line190">190: </a>        te = x[j][i+1];
199
 
<a name="line191">191: </a>        ae = 0.5*(t0 + te);
200
 
<a name="line192">192: </a>        de = PetscPowScalar(ae,beta);
201
 
<a name="line193">193: </a>        fe = de*(te - t0);
202
 
 
203
 
<a name="line195">195: </a>        ts = x[j-1][i];
204
 
<a name="line196">196: </a>        as = 0.5*(t0 + ts);
205
 
<a name="line197">197: </a>        ds = PetscPowScalar(as,beta);
206
 
<a name="line198">198: </a>        fs = ds*(t0 - ts);
207
 
<a name="line199">199: </a>
208
 
<a name="line200">200: </a>        tn = x[j+1][i];
209
 
<a name="line201">201: </a>        an = 0.5*(t0 + tn);
210
 
<a name="line202">202: </a>        dn = PetscPowScalar(an,beta);
211
 
<a name="line203">203: </a>        fn = dn*(tn - t0);
212
 
 
213
 
<a name="line205">205: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (i == 0) {
214
 
 
215
 
<a name="line207">207: </a>        <font color="#B22222">/* left-hand boundary */</font>
216
 
<a name="line208">208: </a>        tw = tleft;
217
 
<a name="line209">209: </a>        aw = 0.5*(t0 + tw);
218
 
<a name="line210">210: </a>        dw = PetscPowScalar(aw,beta);
219
 
<a name="line211">211: </a>        fw = dw*(t0 - tw);
220
 
 
221
 
<a name="line213">213: </a>        te = x[j][i+1];
222
 
<a name="line214">214: </a>        ae = 0.5*(t0 + te);
223
 
<a name="line215">215: </a>        de = PetscPowScalar(ae,beta);
224
 
<a name="line216">216: </a>        fe = de*(te - t0);
225
 
 
226
 
<a name="line218">218: </a>        <font color="#4169E1">if</font> (j &gt; 0) {
227
 
<a name="line219">219: </a>          ts = x[j-1][i];
228
 
<a name="line220">220: </a>          as = 0.5*(t0 + ts);
229
 
<a name="line221">221: </a>          ds = PetscPowScalar(as,beta);
230
 
<a name="line222">222: </a>          fs = ds*(t0 - ts);
231
 
<a name="line223">223: </a>        } <font color="#4169E1">else</font> {
232
 
<a name="line224">224: </a>           fs = zero;
233
 
<a name="line225">225: </a>        }
234
 
 
235
 
<a name="line227">227: </a>        <font color="#4169E1">if</font> (j &lt; my-1) {
236
 
<a name="line228">228: </a>          tn = x[j+1][i];
237
 
<a name="line229">229: </a>          an = 0.5*(t0 + tn);
238
 
<a name="line230">230: </a>          dn = PetscPowScalar(an,beta);
239
 
<a name="line231">231: </a>          fn = dn*(tn - t0);
240
 
<a name="line232">232: </a>        } <font color="#4169E1">else</font> {
241
 
<a name="line233">233: </a>          fn = zero;
242
 
<a name="line234">234: </a>        }
243
 
 
244
 
<a name="line236">236: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (i == mx-1) {
245
 
 
246
 
<a name="line238">238: </a>        <font color="#B22222">/* right-hand boundary */</font>
247
 
<a name="line239">239: </a>        tw = x[j][i-1];
248
 
<a name="line240">240: </a>        aw = 0.5*(t0 + tw);
249
 
<a name="line241">241: </a>        dw = PetscPowScalar(aw,beta);
250
 
<a name="line242">242: </a>        fw = dw*(t0 - tw);
251
 
<a name="line243">243: </a>
252
 
<a name="line244">244: </a>        te = tright;
253
 
<a name="line245">245: </a>        ae = 0.5*(t0 + te);
254
 
<a name="line246">246: </a>        de = PetscPowScalar(ae,beta);
255
 
<a name="line247">247: </a>        fe = de*(te - t0);
256
 
<a name="line248">248: </a>
257
 
<a name="line249">249: </a>        <font color="#4169E1">if</font> (j &gt; 0) {
258
 
<a name="line250">250: </a>          ts = x[j-1][i];
259
 
<a name="line251">251: </a>          as = 0.5*(t0 + ts);
260
 
<a name="line252">252: </a>          ds = PetscPowScalar(as,beta);
261
 
<a name="line253">253: </a>          fs = ds*(t0 - ts);
262
 
<a name="line254">254: </a>        } <font color="#4169E1">else</font> {
263
 
<a name="line255">255: </a>          fs = zero;
264
 
<a name="line256">256: </a>        }
265
 
<a name="line257">257: </a>
266
 
<a name="line258">258: </a>        <font color="#4169E1">if</font> (j &lt; my-1) {
267
 
<a name="line259">259: </a>          tn = x[j+1][i];
268
 
<a name="line260">260: </a>          an = 0.5*(t0 + tn);
269
 
<a name="line261">261: </a>          dn = PetscPowScalar(an,beta);
270
 
<a name="line262">262: </a>          fn = dn*(tn - t0);
271
 
<a name="line263">263: </a>        } <font color="#4169E1">else</font> {
272
 
<a name="line264">264: </a>          fn = zero;
273
 
<a name="line265">265: </a>        }
274
 
 
275
 
<a name="line267">267: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (j == 0) {
276
 
 
277
 
<a name="line269">269: </a>        <font color="#B22222">/* bottom boundary,and i &lt;&gt; 0 or mx-1 */</font>
278
 
<a name="line270">270: </a>        tw = x[j][i-1];
279
 
<a name="line271">271: </a>        aw = 0.5*(t0 + tw);
280
 
<a name="line272">272: </a>        dw = PetscPowScalar(aw,beta);
281
 
<a name="line273">273: </a>        fw = dw*(t0 - tw);
282
 
 
283
 
<a name="line275">275: </a>        te = x[j][i+1];
284
 
<a name="line276">276: </a>        ae = 0.5*(t0 + te);
285
 
<a name="line277">277: </a>        de = PetscPowScalar(ae,beta);
286
 
<a name="line278">278: </a>        fe = de*(te - t0);
287
 
 
288
 
<a name="line280">280: </a>        fs = zero;
289
 
 
290
 
<a name="line282">282: </a>        tn = x[j+1][i];
291
 
<a name="line283">283: </a>        an = 0.5*(t0 + tn);
292
 
<a name="line284">284: </a>        dn = PetscPowScalar(an,beta);
293
 
<a name="line285">285: </a>        fn = dn*(tn - t0);
294
 
 
295
 
<a name="line287">287: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (j == my-1) {
296
 
 
297
 
<a name="line289">289: </a>        <font color="#B22222">/* top boundary,and i &lt;&gt; 0 or mx-1 */</font>
298
 
<a name="line290">290: </a>        tw = x[j][i-1];
299
 
<a name="line291">291: </a>        aw = 0.5*(t0 + tw);
300
 
<a name="line292">292: </a>        dw = PetscPowScalar(aw,beta);
301
 
<a name="line293">293: </a>        fw = dw*(t0 - tw);
302
 
 
303
 
<a name="line295">295: </a>        te = x[j][i+1];
304
 
<a name="line296">296: </a>        ae = 0.5*(t0 + te);
305
 
<a name="line297">297: </a>        de = PetscPowScalar(ae,beta);
306
 
<a name="line298">298: </a>        fe = de*(te - t0);
307
 
 
308
 
<a name="line300">300: </a>        ts = x[j-1][i];
309
 
<a name="line301">301: </a>        as = 0.5*(t0 + ts);
310
 
<a name="line302">302: </a>        ds = PetscPowScalar(as,beta);
311
 
<a name="line303">303: </a>        fs = ds*(t0 - ts);
312
 
 
313
 
<a name="line305">305: </a>        fn = zero;
314
 
 
315
 
<a name="line307">307: </a>      }
316
 
 
317
 
<a name="line309">309: </a>      f[j][i] = - hydhx*(fe-fw) - hxdhy*(fn-fs);
318
 
 
319
 
<a name="line311">311: </a>    }
320
 
<a name="line312">312: </a>  }
321
 
<a name="line313">313: </a>  <A href="../../../../docs/manualpages/DA/DAVecRestoreArray.html#DAVecRestoreArray">DAVecRestoreArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,localX,(void**)&amp;x);
322
 
<a name="line314">314: </a>  <A href="../../../../docs/manualpages/DA/DAVecRestoreArray.html#DAVecRestoreArray">DAVecRestoreArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,F,(void**)&amp;f);
323
 
<a name="line315">315: </a>  <A href="../../../../docs/manualpages/DA/DARestoreLocalVector.html#DARestoreLocalVector">DARestoreLocalVector</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,&amp;localX);
324
 
<a name="line316">316: </a>  <A href="../../../../docs/manualpages/Profiling/PetscLogFlops.html#PetscLogFlops">PetscLogFlops</A>((22 + 4*POWFLOP)*ym*xm);
325
 
<a name="line317">317: </a>  <font color="#4169E1">return</font>(0);
326
 
<a name="line318">318: </a>}
327
 
<a name="line319">319: </a><font color="#B22222">/* --------------------  Evaluate Jacobian F(x) --------------------- */</font>
328
 
<a name="line320">320: </a><strong><font color="#4169E1"><a name="FormJacobian"></a>int FormJacobian(<A href="../../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A> snes,<A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A> X,<A href="../../../../docs/manualpages/Mat/Mat.html#Mat">Mat</A> *J,<A href="../../../../docs/manualpages/Mat/Mat.html#Mat">Mat</A> *B,<A href="../../../../docs/manualpages/Mat/MatStructure.html#MatStructure">MatStructure</A> *flg,void *ptr)</font></strong>
329
 
<a name="line321">321: </a>{
330
 
<a name="line322">322: </a>  <A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A>         dmmg = (<A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A>)ptr;
331
 
<a name="line323">323: </a>  AppCtx       *user = (AppCtx*)dmmg-&gt;user;
332
 
<a name="line324">324: </a>  <A href="../../../../docs/manualpages/Mat/Mat.html#Mat">Mat</A>          jac = *J;
333
 
<a name="line325">325: </a>  int          ierr,i,j,mx,my,xs,ys,xm,ym;
334
 
<a name="line326">326: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
335
 
<a name="line327">327: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
336
 
<a name="line328">328: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  tleft,tright,beta,bm1,coef;
337
 
<a name="line329">329: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  v[5],**x;
338
 
<a name="line330">330: </a>  <A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A>          localX;
339
 
<a name="line331">331: </a>  <A href="../../../../docs/manualpages/Mat/MatStencil.html#MatStencil">MatStencil</A>   col[5],row;
340
 
 
341
 
<a name="line334">334: </a>  <A href="../../../../docs/manualpages/DA/DAGetLocalVector.html#DAGetLocalVector">DAGetLocalVector</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,&amp;localX);
342
 
<a name="line335">335: </a>  *flg = SAME_NONZERO_PATTERN;
343
 
<a name="line336">336: </a>  <A href="../../../../docs/manualpages/DA/DAGetInfo.html#DAGetInfo">DAGetInfo</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,PETSC_NULL,&amp;mx,&amp;my,0,0,0,0,0,0,0,0);
344
 
<a name="line337">337: </a>  hx    = one/(PetscReal)(mx-1);  hy     = one/(PetscReal)(my-1);
345
 
<a name="line338">338: </a>  hxdhy = hx/hy;               hydhx  = hy/hx;
346
 
<a name="line339">339: </a>  tleft = user-&gt;tleft;         tright = user-&gt;tright;
347
 
<a name="line340">340: </a>  beta  = user-&gt;beta;               bm1    = user-&gt;bm1;                coef = user-&gt;coef;
348
 
 
349
 
<a name="line342">342: </a>  <font color="#B22222">/* Get ghost points */</font>
350
 
<a name="line343">343: </a>  <A href="../../../../docs/manualpages/DA/DAGlobalToLocalBegin.html#DAGlobalToLocalBegin">DAGlobalToLocalBegin</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,X,INSERT_VALUES,localX);
351
 
<a name="line344">344: </a>  <A href="../../../../docs/manualpages/DA/DAGlobalToLocalEnd.html#DAGlobalToLocalEnd">DAGlobalToLocalEnd</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,X,INSERT_VALUES,localX);
352
 
<a name="line345">345: </a>  <A href="../../../../docs/manualpages/DA/DAGetCorners.html#DAGetCorners">DAGetCorners</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,&amp;xs,&amp;ys,0,&amp;xm,&amp;ym,0);
353
 
<a name="line346">346: </a>  <A href="../../../../docs/manualpages/DA/DAVecGetArray.html#DAVecGetArray">DAVecGetArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,localX,(void**)&amp;x);
354
 
 
355
 
<a name="line348">348: </a>  <font color="#B22222">/* Evaluate Jacobian of function */</font>
356
 
<a name="line349">349: </a>  <font color="#4169E1">for</font> (j=ys; j&lt;ys+ym; j++) {
357
 
<a name="line350">350: </a>    <font color="#4169E1">for</font> (i=xs; i&lt;xs+xm; i++) {
358
 
<a name="line351">351: </a>      t0 = x[j][i];
359
 
 
360
 
<a name="line353">353: </a>      <font color="#4169E1">if</font> (i &gt; 0 &amp;&amp; i &lt; mx-1 &amp;&amp; j &gt; 0 &amp;&amp; j &lt; my-1) {
361
 
 
362
 
<a name="line355">355: </a>        <font color="#B22222">/* general interior volume */</font>
363
 
 
364
 
<a name="line357">357: </a>        tw = x[j][i-1];
365
 
<a name="line358">358: </a>        aw = 0.5*(t0 + tw);
366
 
<a name="line359">359: </a>        bw = PetscPowScalar(aw,bm1);
367
 
<a name="line360">360: </a>        <font color="#B22222">/* dw = bw * aw */</font>
368
 
<a name="line361">361: </a>        dw = PetscPowScalar(aw,beta);
369
 
<a name="line362">362: </a>        gw = coef*bw*(t0 - tw);
370
 
 
371
 
<a name="line364">364: </a>        te = x[j][i+1];
372
 
<a name="line365">365: </a>        ae = 0.5*(t0 + te);
373
 
<a name="line366">366: </a>        be = PetscPowScalar(ae,bm1);
374
 
<a name="line367">367: </a>        <font color="#B22222">/* de = be * ae; */</font>
375
 
<a name="line368">368: </a>        de = PetscPowScalar(ae,beta);
376
 
<a name="line369">369: </a>        ge = coef*be*(te - t0);
377
 
 
378
 
<a name="line371">371: </a>        ts = x[j-1][i];
379
 
<a name="line372">372: </a>        as = 0.5*(t0 + ts);
380
 
<a name="line373">373: </a>        bs = PetscPowScalar(as,bm1);
381
 
<a name="line374">374: </a>        <font color="#B22222">/* ds = bs * as; */</font>
382
 
<a name="line375">375: </a>        ds = PetscPowScalar(as,beta);
383
 
<a name="line376">376: </a>        gs = coef*bs*(t0 - ts);
384
 
<a name="line377">377: </a>
385
 
<a name="line378">378: </a>        tn = x[j+1][i];
386
 
<a name="line379">379: </a>        an = 0.5*(t0 + tn);
387
 
<a name="line380">380: </a>        bn = PetscPowScalar(an,bm1);
388
 
<a name="line381">381: </a>        <font color="#B22222">/* dn = bn * an; */</font>
389
 
<a name="line382">382: </a>        dn = PetscPowScalar(an,beta);
390
 
<a name="line383">383: </a>        gn = coef*bn*(tn - t0);
391
 
 
392
 
<a name="line385">385: </a>        v[0] = - hxdhy*(ds - gs);                                      col[0].j = j-1;       col[0].i = i;
393
 
<a name="line386">386: </a>        v[1] = - hydhx*(dw - gw);                                      col[1].j = j;         col[1].i = i-1;
394
 
<a name="line387">387: </a>        v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
395
 
<a name="line388">388: </a>        v[3] = - hydhx*(de + ge);                                      col[3].j = j;         col[3].i = i+1;
396
 
<a name="line389">389: </a>        v[4] = - hxdhy*(dn + gn);                                      col[4].j = j+1;       col[4].i = i;
397
 
<a name="line390">390: </a>        <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,5,col,v,INSERT_VALUES);
398
 
 
399
 
<a name="line392">392: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (i == 0) {
400
 
 
401
 
<a name="line394">394: </a>        <font color="#B22222">/* left-hand boundary */</font>
402
 
<a name="line395">395: </a>        tw = tleft;
403
 
<a name="line396">396: </a>        aw = 0.5*(t0 + tw);
404
 
<a name="line397">397: </a>        bw = PetscPowScalar(aw,bm1);
405
 
<a name="line398">398: </a>        <font color="#B22222">/* dw = bw * aw */</font>
406
 
<a name="line399">399: </a>        dw = PetscPowScalar(aw,beta);
407
 
<a name="line400">400: </a>        gw = coef*bw*(t0 - tw);
408
 
<a name="line401">401: </a>
409
 
<a name="line402">402: </a>        te = x[j][i + 1];
410
 
<a name="line403">403: </a>        ae = 0.5*(t0 + te);
411
 
<a name="line404">404: </a>        be = PetscPowScalar(ae,bm1);
412
 
<a name="line405">405: </a>        <font color="#B22222">/* de = be * ae; */</font>
413
 
<a name="line406">406: </a>        de = PetscPowScalar(ae,beta);
414
 
<a name="line407">407: </a>        ge = coef*be*(te - t0);
415
 
<a name="line408">408: </a>
416
 
<a name="line409">409: </a>        <font color="#B22222">/* left-hand bottom boundary */</font>
417
 
<a name="line410">410: </a>        <font color="#4169E1">if</font> (j == 0) {
418
 
 
419
 
<a name="line412">412: </a>          tn = x[j+1][i];
420
 
<a name="line413">413: </a>          an = 0.5*(t0 + tn);
421
 
<a name="line414">414: </a>          bn = PetscPowScalar(an,bm1);
422
 
<a name="line415">415: </a>          <font color="#B22222">/* dn = bn * an; */</font>
423
 
<a name="line416">416: </a>          dn = PetscPowScalar(an,beta);
424
 
<a name="line417">417: </a>          gn = coef*bn*(tn - t0);
425
 
<a name="line418">418: </a>
426
 
<a name="line419">419: </a>          v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;
427
 
<a name="line420">420: </a>          v[1] = - hydhx*(de + ge);                           col[1].j = j;         col[1].i = i+1;
428
 
<a name="line421">421: </a>          v[2] = - hxdhy*(dn + gn);                           col[2].j = j+1;       col[2].i = i;
429
 
<a name="line422">422: </a>          <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,3,col,v,INSERT_VALUES);
430
 
<a name="line423">423: </a>
431
 
<a name="line424">424: </a>        <font color="#B22222">/* left-hand interior boundary */</font>
432
 
<a name="line425">425: </a>        } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (j &lt; my-1) {
433
 
 
434
 
<a name="line427">427: </a>          ts = x[j-1][i];
435
 
<a name="line428">428: </a>          as = 0.5*(t0 + ts);
436
 
<a name="line429">429: </a>          bs = PetscPowScalar(as,bm1);
437
 
<a name="line430">430: </a>          <font color="#B22222">/* ds = bs * as; */</font>
438
 
<a name="line431">431: </a>          ds = PetscPowScalar(as,beta);
439
 
<a name="line432">432: </a>          gs = coef*bs*(t0 - ts);
440
 
<a name="line433">433: </a>
441
 
<a name="line434">434: </a>          tn = x[j+1][i];
442
 
<a name="line435">435: </a>          an = 0.5*(t0 + tn);
443
 
<a name="line436">436: </a>          bn = PetscPowScalar(an,bm1);
444
 
<a name="line437">437: </a>          <font color="#B22222">/* dn = bn * an; */</font>
445
 
<a name="line438">438: </a>          dn = PetscPowScalar(an,beta);
446
 
<a name="line439">439: </a>          gn = coef*bn*(tn - t0);
447
 
<a name="line440">440: </a>
448
 
<a name="line441">441: </a>          v[0] = - hxdhy*(ds - gs);                                      col[0].j = j-1;       col[0].i = i;
449
 
<a name="line442">442: </a>          v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
450
 
<a name="line443">443: </a>          v[2] = - hydhx*(de + ge);                                      col[2].j = j;         col[2].i = i+1;
451
 
<a name="line444">444: </a>          v[3] = - hxdhy*(dn + gn);                                      col[3].j = j+1;       col[3].i = i;
452
 
<a name="line445">445: </a>          <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,4,col,v,INSERT_VALUES);
453
 
<a name="line446">446: </a>        <font color="#B22222">/* left-hand top boundary */</font>
454
 
<a name="line447">447: </a>        } <font color="#4169E1">else</font> {
455
 
 
456
 
<a name="line449">449: </a>          ts = x[j-1][i];
457
 
<a name="line450">450: </a>          as = 0.5*(t0 + ts);
458
 
<a name="line451">451: </a>          bs = PetscPowScalar(as,bm1);
459
 
<a name="line452">452: </a>          <font color="#B22222">/* ds = bs * as; */</font>
460
 
<a name="line453">453: </a>          ds = PetscPowScalar(as,beta);
461
 
<a name="line454">454: </a>          gs = coef*bs*(t0 - ts);
462
 
<a name="line455">455: </a>
463
 
<a name="line456">456: </a>          v[0] = - hxdhy*(ds - gs);                            col[0].j = j-1;       col[0].i = i;
464
 
<a name="line457">457: </a>          v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
465
 
<a name="line458">458: </a>          v[2] = - hydhx*(de + ge);                            col[2].j = j;         col[2].i = i+1;
466
 
<a name="line459">459: </a>          <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,3,col,v,INSERT_VALUES);
467
 
<a name="line460">460: </a>        }
468
 
 
469
 
<a name="line462">462: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (i == mx-1) {
 
71
<a name="line63"> 63: </a><strong><font color="#4169E1"><a name="main"></a>int main(int argc,char **argv)</font></strong>
 
72
<a name="line64"> 64: </a>{
 
73
<a name="line65"> 65: </a>  <A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A>      *dmmg;
 
74
<a name="line66"> 66: </a>  <A href="../../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A>      snes;
 
75
<a name="line67"> 67: </a>  AppCtx    user;
 
76
<a name="line68"> 68: </a>  int       ierr,its,lits;
 
77
<a name="line69"> 69: </a>  <A href="../../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A> litspit;
 
78
<a name="line70"> 70: </a>  <A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>        da;
 
79
 
 
80
<a name="line72"> 72: </a>  <A href="../../../../docs/manualpages/Sys/PetscInitialize.html#PetscInitialize">PetscInitialize</A>(&amp;argc,&amp;argv,<A href="../../../../docs/manualpages/Sys/PETSC_NULL.html#PETSC_NULL">PETSC_NULL</A>,help);
 
81
 
 
82
<a name="line74"> 74: </a>  <font color="#B22222">/* set problem parameters */</font>
 
83
<a name="line75"> 75: </a>  user.tleft  = 1.0;
 
84
<a name="line76"> 76: </a>  user.tright = 0.1;
 
85
<a name="line77"> 77: </a>  user.beta   = 2.5;
 
86
<a name="line78"> 78: </a>  <A href="../../../../docs/manualpages/Sys/PetscOptionsGetReal.html#PetscOptionsGetReal">PetscOptionsGetReal</A>(<A href="../../../../docs/manualpages/Sys/PETSC_NULL.html#PETSC_NULL">PETSC_NULL</A>,<font color="#666666">"-tleft"</font>,&amp;user.tleft,<A href="../../../../docs/manualpages/Sys/PETSC_NULL.html#PETSC_NULL">PETSC_NULL</A>);
 
87
<a name="line79"> 79: </a>  <A href="../../../../docs/manualpages/Sys/PetscOptionsGetReal.html#PetscOptionsGetReal">PetscOptionsGetReal</A>(<A href="../../../../docs/manualpages/Sys/PETSC_NULL.html#PETSC_NULL">PETSC_NULL</A>,<font color="#666666">"-tright"</font>,&amp;user.tright,<A href="../../../../docs/manualpages/Sys/PETSC_NULL.html#PETSC_NULL">PETSC_NULL</A>);
 
88
<a name="line80"> 80: </a>  <A href="../../../../docs/manualpages/Sys/PetscOptionsGetReal.html#PetscOptionsGetReal">PetscOptionsGetReal</A>(<A href="../../../../docs/manualpages/Sys/PETSC_NULL.html#PETSC_NULL">PETSC_NULL</A>,<font color="#666666">"-beta"</font>,&amp;user.beta,<A href="../../../../docs/manualpages/Sys/PETSC_NULL.html#PETSC_NULL">PETSC_NULL</A>);
 
89
<a name="line81"> 81: </a>  user.bm1  = user.beta - 1.0;
 
90
<a name="line82"> 82: </a>  user.coef = user.beta/2.0;
 
91
 
 
92
 
 
93
<a name="line85"> 85: </a>  <font color="#B22222">/*</font>
 
94
<a name="line86"> 86: </a><font color="#B22222">      Create the multilevel <A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A> data structure </font>
 
95
<a name="line87"> 87: </a><font color="#B22222">  */</font>
 
96
<a name="line88"> 88: </a>  <A href="../../../../docs/manualpages/KSP/DMMGCreate.html#DMMGCreate">DMMGCreate</A>(<A href="../../../../docs/manualpages/Sys/PETSC_COMM_WORLD.html#PETSC_COMM_WORLD">PETSC_COMM_WORLD</A>,3,&amp;user,&amp;dmmg);
 
97
 
 
98
<a name="line90"> 90: </a>  <font color="#B22222">/*</font>
 
99
<a name="line91"> 91: </a><font color="#B22222">      Set the <A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A> (grid structure) for the grids.</font>
 
100
<a name="line92"> 92: </a><font color="#B22222">  */</font>
 
101
<a name="line93"> 93: </a>  <A href="../../../../docs/manualpages/DA/DACreate2d.html#DACreate2d">DACreate2d</A>(<A href="../../../../docs/manualpages/Sys/PETSC_COMM_WORLD.html#PETSC_COMM_WORLD">PETSC_COMM_WORLD</A>,DA_NONPERIODIC,DA_STENCIL_STAR,5,5,<A href="../../../../docs/manualpages/Sys/PETSC_DECIDE.html#PETSC_DECIDE">PETSC_DECIDE</A>,<A href="../../../../docs/manualpages/Sys/PETSC_DECIDE.html#PETSC_DECIDE">PETSC_DECIDE</A>,1,1,0,0,&amp;da);
 
102
<a name="line94"> 94: </a>  <A href="../../../../docs/manualpages/KSP/DMMGSetDM.html#DMMGSetDM">DMMGSetDM</A>(dmmg,(<A href="../../../../docs/manualpages/DA/DM.html#DM">DM</A>)da);
 
103
<a name="line95"> 95: </a>  <A href="../../../../docs/manualpages/DA/DADestroy.html#DADestroy">DADestroy</A>(da);
 
104
 
 
105
<a name="line97"> 97: </a>  <font color="#B22222">/*</font>
 
106
<a name="line98"> 98: </a><font color="#B22222">     Create the nonlinear solver, and tell the <A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A> structure to use it</font>
 
107
<a name="line99"> 99: </a><font color="#B22222">  */</font>
 
108
<a name="line100">100: </a>  <A href="../../../../docs/manualpages/SNES/DMMGSetSNES.html#DMMGSetSNES">DMMGSetSNES</A>(dmmg,FormFunction,FormJacobian);
 
109
 
 
110
<a name="line102">102: </a>  <font color="#B22222">/*</font>
 
111
<a name="line103">103: </a><font color="#B22222">      <A href="../../../../docs/manualpages/Profiling/PreLoadBegin.html#PreLoadBegin">PreLoadBegin</A>() means that the following section of code is run twice. The first time</font>
 
112
<a name="line104">104: </a><font color="#B22222">     through the flag PreLoading is on this the nonlinear solver is only run for a single step.</font>
 
113
<a name="line105">105: </a><font color="#B22222">     The second time through (the actually timed code) the maximum iterations is set to 10</font>
 
114
<a name="line106">106: </a><font color="#B22222">     Preload of the executable is done to eliminate from the timing the time spent bring the </font>
 
115
<a name="line107">107: </a><font color="#B22222">     executable into memory from disk (paging in).</font>
 
116
<a name="line108">108: </a><font color="#B22222">  */</font>
 
117
<a name="line109">109: </a>  <A href="../../../../docs/manualpages/Profiling/PreLoadBegin.html#PreLoadBegin">PreLoadBegin</A>(PETSC_TRUE,<font color="#666666">"Solve"</font>);
 
118
<a name="line110">110: </a>    <A href="../../../../docs/manualpages/SNES/DMMGSetInitialGuess.html#DMMGSetInitialGuess">DMMGSetInitialGuess</A>(dmmg,FormInitialGuess);
 
119
<a name="line111">111: </a>    <A href="../../../../docs/manualpages/KSP/DMMGSolve.html#DMMGSolve">DMMGSolve</A>(dmmg);
 
120
<a name="line112">112: </a>  <A href="../../../../docs/manualpages/Profiling/PreLoadEnd.html#PreLoadEnd">PreLoadEnd</A>();
 
121
<a name="line113">113: </a>  snes = DMMGGetSNES(dmmg);
 
122
<a name="line114">114: </a>  <A href="../../../../docs/manualpages/SNES/SNESGetIterationNumber.html#SNESGetIterationNumber">SNESGetIterationNumber</A>(snes,&amp;its);
 
123
<a name="line115">115: </a>  <A href="../../../../docs/manualpages/SNES/SNESGetNumberLinearIterations.html#SNESGetNumberLinearIterations">SNESGetNumberLinearIterations</A>(snes,&amp;lits);
 
124
<a name="line116">116: </a>  litspit = ((<A href="../../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A>)lits)/((<A href="../../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A>)its);
 
125
<a name="line117">117: </a>  <A href="../../../../docs/manualpages/Sys/PetscPrintf.html#PetscPrintf">PetscPrintf</A>(<A href="../../../../docs/manualpages/Sys/PETSC_COMM_WORLD.html#PETSC_COMM_WORLD">PETSC_COMM_WORLD</A>,<font color="#666666">"Number of Newton iterations = %d\n"</font>,its);
 
126
<a name="line118">118: </a>  <A href="../../../../docs/manualpages/Sys/PetscPrintf.html#PetscPrintf">PetscPrintf</A>(<A href="../../../../docs/manualpages/Sys/PETSC_COMM_WORLD.html#PETSC_COMM_WORLD">PETSC_COMM_WORLD</A>,<font color="#666666">"Number of Linear iterations = %d\n"</font>,lits);
 
127
<a name="line119">119: </a>  <A href="../../../../docs/manualpages/Sys/PetscPrintf.html#PetscPrintf">PetscPrintf</A>(<A href="../../../../docs/manualpages/Sys/PETSC_COMM_WORLD.html#PETSC_COMM_WORLD">PETSC_COMM_WORLD</A>,<font color="#666666">"Average Linear its / Newton = %e\n"</font>,litspit);
 
128
 
 
129
<a name="line121">121: </a>  <A href="../../../../docs/manualpages/KSP/DMMGDestroy.html#DMMGDestroy">DMMGDestroy</A>(dmmg);
 
130
<a name="line122">122: </a>  <A href="../../../../docs/manualpages/Sys/PetscFinalize.html#PetscFinalize">PetscFinalize</A>();
 
131
 
 
132
<a name="line124">124: </a>  <font color="#4169E1">return</font> 0;
 
133
<a name="line125">125: </a>}
 
134
<a name="line126">126: </a><font color="#B22222">/* --------------------  Form initial approximation ----------------- */</font>
 
135
<a name="line129">129: </a><strong><font color="#4169E1"><a name="FormInitialGuess"></a>int FormInitialGuess(<A href="../../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A> snes,<A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A> X,void *ptr)</font></strong>
 
136
<a name="line130">130: </a>{
 
137
<a name="line131">131: </a>  <A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A>        dmmg = (<A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A>)ptr;
 
138
<a name="line132">132: </a>  AppCtx      *user = (AppCtx*)dmmg-&gt;user;
 
139
<a name="line133">133: </a>  int         i,j,ierr,xs,ys,xm,ym;
 
140
<a name="line134">134: </a>  <A href="../../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A>   tleft = user-&gt;tleft;
 
141
<a name="line135">135: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A> **x;
 
142
 
 
143
 
 
144
<a name="line139">139: </a>  <font color="#B22222">/* Get ghost points */</font>
 
145
<a name="line140">140: </a>  <A href="../../../../docs/manualpages/DA/DAGetCorners.html#DAGetCorners">DAGetCorners</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,&amp;xs,&amp;ys,0,&amp;xm,&amp;ym,0);
 
146
<a name="line141">141: </a>  <A href="../../../../docs/manualpages/DA/DAVecGetArray.html#DAVecGetArray">DAVecGetArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,X,&amp;x);
 
147
 
 
148
<a name="line143">143: </a>  <font color="#B22222">/* Compute initial guess */</font>
 
149
<a name="line144">144: </a>  <font color="#4169E1">for</font> (j=ys; j&lt;ys+ym; j++) {
 
150
<a name="line145">145: </a>    <font color="#4169E1">for</font> (i=xs; i&lt;xs+xm; i++) {
 
151
<a name="line146">146: </a>      x[j][i] = tleft;
 
152
<a name="line147">147: </a>    }
 
153
<a name="line148">148: </a>  }
 
154
<a name="line149">149: </a>  <A href="../../../../docs/manualpages/DA/DAVecRestoreArray.html#DAVecRestoreArray">DAVecRestoreArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,X,&amp;x);
 
155
<a name="line150">150: </a>  <font color="#4169E1">return</font>(0);
 
156
<a name="line151">151: </a>}
 
157
<a name="line152">152: </a><font color="#B22222">/* --------------------  Evaluate Function F(x) --------------------- */</font>
 
158
<a name="line155">155: </a><strong><font color="#4169E1"><a name="FormFunction"></a>int FormFunction(<A href="../../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A> snes,<A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A> X,<A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A> F,void* ptr)</font></strong>
 
159
<a name="line156">156: </a>{
 
160
<a name="line157">157: </a>  <A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A>    dmmg = (<A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A>)ptr;
 
161
<a name="line158">158: </a>  AppCtx  *user = (AppCtx*)dmmg-&gt;user;
 
162
<a name="line159">159: </a>  int     ierr,i,j,mx,my,xs,ys,xm,ym;
 
163
<a name="line160">160: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  zero = 0.0,one = 1.0;
 
164
<a name="line161">161: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  hx,hy,hxdhy,hydhx;
 
165
<a name="line162">162: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw,fn = 0.0,fs = 0.0,fe =0.0,fw = 0.0;
 
166
<a name="line163">163: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  tleft,tright,beta;
 
167
<a name="line164">164: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  **x,**f;
 
168
<a name="line165">165: </a>  <A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A>     localX;
 
169
 
 
170
<a name="line168">168: </a>  <A href="../../../../docs/manualpages/DA/DAGetLocalVector.html#DAGetLocalVector">DAGetLocalVector</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,&amp;localX);
 
171
<a name="line169">169: </a>  <A href="../../../../docs/manualpages/DA/DAGetInfo.html#DAGetInfo">DAGetInfo</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,<A href="../../../../docs/manualpages/Sys/PETSC_NULL.html#PETSC_NULL">PETSC_NULL</A>,&amp;mx,&amp;my,0,0,0,0,0,0,0,0);
 
172
<a name="line170">170: </a>  hx    = one/(<A href="../../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A>)(mx-1);  hy    = one/(<A href="../../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A>)(my-1);
 
173
<a name="line171">171: </a>  hxdhy = hx/hy;               hydhx = hy/hx;
 
174
<a name="line172">172: </a>  tleft = user-&gt;tleft;         tright = user-&gt;tright;
 
175
<a name="line173">173: </a>  beta  = user-&gt;beta;
 
176
<a name="line174">174: </a>
 
177
<a name="line175">175: </a>  <font color="#B22222">/* Get ghost points */</font>
 
178
<a name="line176">176: </a>  <A href="../../../../docs/manualpages/DA/DAGlobalToLocalBegin.html#DAGlobalToLocalBegin">DAGlobalToLocalBegin</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,X,INSERT_VALUES,localX);
 
179
<a name="line177">177: </a>  <A href="../../../../docs/manualpages/DA/DAGlobalToLocalEnd.html#DAGlobalToLocalEnd">DAGlobalToLocalEnd</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,X,INSERT_VALUES,localX);
 
180
<a name="line178">178: </a>  <A href="../../../../docs/manualpages/DA/DAGetCorners.html#DAGetCorners">DAGetCorners</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,&amp;xs,&amp;ys,0,&amp;xm,&amp;ym,0);
 
181
<a name="line179">179: </a>  <A href="../../../../docs/manualpages/DA/DAVecGetArray.html#DAVecGetArray">DAVecGetArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,localX,&amp;x);
 
182
<a name="line180">180: </a>  <A href="../../../../docs/manualpages/DA/DAVecGetArray.html#DAVecGetArray">DAVecGetArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,F,&amp;f);
 
183
 
 
184
<a name="line182">182: </a>  <font color="#B22222">/* Evaluate function */</font>
 
185
<a name="line183">183: </a>  <font color="#4169E1">for</font> (j=ys; j&lt;ys+ym; j++) {
 
186
<a name="line184">184: </a>    <font color="#4169E1">for</font> (i=xs; i&lt;xs+xm; i++) {
 
187
<a name="line185">185: </a>      t0 = x[j][i];
 
188
 
 
189
<a name="line187">187: </a>      <font color="#4169E1">if</font> (i &gt; 0 &amp;&amp; i &lt; mx-1 &amp;&amp; j &gt; 0 &amp;&amp; j &lt; my-1) {
 
190
 
 
191
<a name="line189">189: </a>        <font color="#B22222">/* general interior volume */</font>
 
192
 
 
193
<a name="line191">191: </a>        tw = x[j][i-1];
 
194
<a name="line192">192: </a>        aw = 0.5*(t0 + tw);
 
195
<a name="line193">193: </a>        dw = PetscPowScalar(aw,beta);
 
196
<a name="line194">194: </a>        fw = dw*(t0 - tw);
 
197
 
 
198
<a name="line196">196: </a>        te = x[j][i+1];
 
199
<a name="line197">197: </a>        ae = 0.5*(t0 + te);
 
200
<a name="line198">198: </a>        de = PetscPowScalar(ae,beta);
 
201
<a name="line199">199: </a>        fe = de*(te - t0);
 
202
 
 
203
<a name="line201">201: </a>        ts = x[j-1][i];
 
204
<a name="line202">202: </a>        as = 0.5*(t0 + ts);
 
205
<a name="line203">203: </a>        ds = PetscPowScalar(as,beta);
 
206
<a name="line204">204: </a>        fs = ds*(t0 - ts);
 
207
<a name="line205">205: </a>
 
208
<a name="line206">206: </a>        tn = x[j+1][i];
 
209
<a name="line207">207: </a>        an = 0.5*(t0 + tn);
 
210
<a name="line208">208: </a>        dn = PetscPowScalar(an,beta);
 
211
<a name="line209">209: </a>        fn = dn*(tn - t0);
 
212
 
 
213
<a name="line211">211: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (i == 0) {
 
214
 
 
215
<a name="line213">213: </a>        <font color="#B22222">/* left-hand boundary */</font>
 
216
<a name="line214">214: </a>        tw = tleft;
 
217
<a name="line215">215: </a>        aw = 0.5*(t0 + tw);
 
218
<a name="line216">216: </a>        dw = PetscPowScalar(aw,beta);
 
219
<a name="line217">217: </a>        fw = dw*(t0 - tw);
 
220
 
 
221
<a name="line219">219: </a>        te = x[j][i+1];
 
222
<a name="line220">220: </a>        ae = 0.5*(t0 + te);
 
223
<a name="line221">221: </a>        de = PetscPowScalar(ae,beta);
 
224
<a name="line222">222: </a>        fe = de*(te - t0);
 
225
 
 
226
<a name="line224">224: </a>        <font color="#4169E1">if</font> (j &gt; 0) {
 
227
<a name="line225">225: </a>          ts = x[j-1][i];
 
228
<a name="line226">226: </a>          as = 0.5*(t0 + ts);
 
229
<a name="line227">227: </a>          ds = PetscPowScalar(as,beta);
 
230
<a name="line228">228: </a>          fs = ds*(t0 - ts);
 
231
<a name="line229">229: </a>        } <font color="#4169E1">else</font> {
 
232
<a name="line230">230: </a>           fs = zero;
 
233
<a name="line231">231: </a>        }
 
234
 
 
235
<a name="line233">233: </a>        <font color="#4169E1">if</font> (j &lt; my-1) {
 
236
<a name="line234">234: </a>          tn = x[j+1][i];
 
237
<a name="line235">235: </a>          an = 0.5*(t0 + tn);
 
238
<a name="line236">236: </a>          dn = PetscPowScalar(an,beta);
 
239
<a name="line237">237: </a>          fn = dn*(tn - t0);
 
240
<a name="line238">238: </a>        } <font color="#4169E1">else</font> {
 
241
<a name="line239">239: </a>          fn = zero;
 
242
<a name="line240">240: </a>        }
 
243
 
 
244
<a name="line242">242: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (i == mx-1) {
 
245
 
 
246
<a name="line244">244: </a>        <font color="#B22222">/* right-hand boundary */</font>
 
247
<a name="line245">245: </a>        tw = x[j][i-1];
 
248
<a name="line246">246: </a>        aw = 0.5*(t0 + tw);
 
249
<a name="line247">247: </a>        dw = PetscPowScalar(aw,beta);
 
250
<a name="line248">248: </a>        fw = dw*(t0 - tw);
 
251
<a name="line249">249: </a>
 
252
<a name="line250">250: </a>        te = tright;
 
253
<a name="line251">251: </a>        ae = 0.5*(t0 + te);
 
254
<a name="line252">252: </a>        de = PetscPowScalar(ae,beta);
 
255
<a name="line253">253: </a>        fe = de*(te - t0);
 
256
<a name="line254">254: </a>
 
257
<a name="line255">255: </a>        <font color="#4169E1">if</font> (j &gt; 0) {
 
258
<a name="line256">256: </a>          ts = x[j-1][i];
 
259
<a name="line257">257: </a>          as = 0.5*(t0 + ts);
 
260
<a name="line258">258: </a>          ds = PetscPowScalar(as,beta);
 
261
<a name="line259">259: </a>          fs = ds*(t0 - ts);
 
262
<a name="line260">260: </a>        } <font color="#4169E1">else</font> {
 
263
<a name="line261">261: </a>          fs = zero;
 
264
<a name="line262">262: </a>        }
 
265
<a name="line263">263: </a>
 
266
<a name="line264">264: </a>        <font color="#4169E1">if</font> (j &lt; my-1) {
 
267
<a name="line265">265: </a>          tn = x[j+1][i];
 
268
<a name="line266">266: </a>          an = 0.5*(t0 + tn);
 
269
<a name="line267">267: </a>          dn = PetscPowScalar(an,beta);
 
270
<a name="line268">268: </a>          fn = dn*(tn - t0);
 
271
<a name="line269">269: </a>        } <font color="#4169E1">else</font> {
 
272
<a name="line270">270: </a>          fn = zero;
 
273
<a name="line271">271: </a>        }
 
274
 
 
275
<a name="line273">273: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (j == 0) {
 
276
 
 
277
<a name="line275">275: </a>        <font color="#B22222">/* bottom boundary,and i &lt;&gt; 0 or mx-1 */</font>
 
278
<a name="line276">276: </a>        tw = x[j][i-1];
 
279
<a name="line277">277: </a>        aw = 0.5*(t0 + tw);
 
280
<a name="line278">278: </a>        dw = PetscPowScalar(aw,beta);
 
281
<a name="line279">279: </a>        fw = dw*(t0 - tw);
 
282
 
 
283
<a name="line281">281: </a>        te = x[j][i+1];
 
284
<a name="line282">282: </a>        ae = 0.5*(t0 + te);
 
285
<a name="line283">283: </a>        de = PetscPowScalar(ae,beta);
 
286
<a name="line284">284: </a>        fe = de*(te - t0);
 
287
 
 
288
<a name="line286">286: </a>        fs = zero;
 
289
 
 
290
<a name="line288">288: </a>        tn = x[j+1][i];
 
291
<a name="line289">289: </a>        an = 0.5*(t0 + tn);
 
292
<a name="line290">290: </a>        dn = PetscPowScalar(an,beta);
 
293
<a name="line291">291: </a>        fn = dn*(tn - t0);
 
294
 
 
295
<a name="line293">293: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (j == my-1) {
 
296
 
 
297
<a name="line295">295: </a>        <font color="#B22222">/* top boundary,and i &lt;&gt; 0 or mx-1 */</font>
 
298
<a name="line296">296: </a>        tw = x[j][i-1];
 
299
<a name="line297">297: </a>        aw = 0.5*(t0 + tw);
 
300
<a name="line298">298: </a>        dw = PetscPowScalar(aw,beta);
 
301
<a name="line299">299: </a>        fw = dw*(t0 - tw);
 
302
 
 
303
<a name="line301">301: </a>        te = x[j][i+1];
 
304
<a name="line302">302: </a>        ae = 0.5*(t0 + te);
 
305
<a name="line303">303: </a>        de = PetscPowScalar(ae,beta);
 
306
<a name="line304">304: </a>        fe = de*(te - t0);
 
307
 
 
308
<a name="line306">306: </a>        ts = x[j-1][i];
 
309
<a name="line307">307: </a>        as = 0.5*(t0 + ts);
 
310
<a name="line308">308: </a>        ds = PetscPowScalar(as,beta);
 
311
<a name="line309">309: </a>        fs = ds*(t0 - ts);
 
312
 
 
313
<a name="line311">311: </a>        fn = zero;
 
314
 
 
315
<a name="line313">313: </a>      }
 
316
 
 
317
<a name="line315">315: </a>      f[j][i] = - hydhx*(fe-fw) - hxdhy*(fn-fs);
 
318
 
 
319
<a name="line317">317: </a>    }
 
320
<a name="line318">318: </a>  }
 
321
<a name="line319">319: </a>  <A href="../../../../docs/manualpages/DA/DAVecRestoreArray.html#DAVecRestoreArray">DAVecRestoreArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,localX,&amp;x);
 
322
<a name="line320">320: </a>  <A href="../../../../docs/manualpages/DA/DAVecRestoreArray.html#DAVecRestoreArray">DAVecRestoreArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,F,&amp;f);
 
323
<a name="line321">321: </a>  <A href="../../../../docs/manualpages/DA/DARestoreLocalVector.html#DARestoreLocalVector">DARestoreLocalVector</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,&amp;localX);
 
324
<a name="line322">322: </a>  <A href="../../../../docs/manualpages/Profiling/PetscLogFlops.html#PetscLogFlops">PetscLogFlops</A>((22 + 4*POWFLOP)*ym*xm);
 
325
<a name="line323">323: </a>  <font color="#4169E1">return</font>(0);
 
326
<a name="line324">324: </a>}
 
327
<a name="line325">325: </a><font color="#B22222">/* --------------------  Evaluate Jacobian F(x) --------------------- */</font>
 
328
<a name="line328">328: </a><strong><font color="#4169E1"><a name="FormJacobian"></a>int FormJacobian(<A href="../../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A> snes,<A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A> X,<A href="../../../../docs/manualpages/Mat/Mat.html#Mat">Mat</A> *J,<A href="../../../../docs/manualpages/Mat/Mat.html#Mat">Mat</A> *B,<A href="../../../../docs/manualpages/Mat/MatStructure.html#MatStructure">MatStructure</A> *flg,void *ptr)</font></strong>
 
329
<a name="line329">329: </a>{
 
330
<a name="line330">330: </a>  <A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A>         dmmg = (<A href="../../../../docs/manualpages/DA/DMMG.html#DMMG">DMMG</A>)ptr;
 
331
<a name="line331">331: </a>  AppCtx       *user = (AppCtx*)dmmg-&gt;user;
 
332
<a name="line332">332: </a>  <A href="../../../../docs/manualpages/Mat/Mat.html#Mat">Mat</A>          jac = *J;
 
333
<a name="line333">333: </a>  int          ierr,i,j,mx,my,xs,ys,xm,ym;
 
334
<a name="line334">334: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
 
335
<a name="line335">335: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
 
336
<a name="line336">336: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  tleft,tright,beta,bm1,coef;
 
337
<a name="line337">337: </a>  <A href="../../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A>  v[5],**x;
 
338
<a name="line338">338: </a>  <A href="../../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A>          localX;
 
339
<a name="line339">339: </a>  <A href="../../../../docs/manualpages/Mat/MatStencil.html#MatStencil">MatStencil</A>   col[5],row;
 
340
 
 
341
<a name="line342">342: </a>  <A href="../../../../docs/manualpages/DA/DAGetLocalVector.html#DAGetLocalVector">DAGetLocalVector</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,&amp;localX);
 
342
<a name="line343">343: </a>  *flg = SAME_NONZERO_PATTERN;
 
343
<a name="line344">344: </a>  <A href="../../../../docs/manualpages/DA/DAGetInfo.html#DAGetInfo">DAGetInfo</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,<A href="../../../../docs/manualpages/Sys/PETSC_NULL.html#PETSC_NULL">PETSC_NULL</A>,&amp;mx,&amp;my,0,0,0,0,0,0,0,0);
 
344
<a name="line345">345: </a>  hx    = one/(<A href="../../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A>)(mx-1);  hy     = one/(<A href="../../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A>)(my-1);
 
345
<a name="line346">346: </a>  hxdhy = hx/hy;               hydhx  = hy/hx;
 
346
<a name="line347">347: </a>  tleft = user-&gt;tleft;         tright = user-&gt;tright;
 
347
<a name="line348">348: </a>  beta  = user-&gt;beta;               bm1    = user-&gt;bm1;                coef = user-&gt;coef;
 
348
 
 
349
<a name="line350">350: </a>  <font color="#B22222">/* Get ghost points */</font>
 
350
<a name="line351">351: </a>  <A href="../../../../docs/manualpages/DA/DAGlobalToLocalBegin.html#DAGlobalToLocalBegin">DAGlobalToLocalBegin</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,X,INSERT_VALUES,localX);
 
351
<a name="line352">352: </a>  <A href="../../../../docs/manualpages/DA/DAGlobalToLocalEnd.html#DAGlobalToLocalEnd">DAGlobalToLocalEnd</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,X,INSERT_VALUES,localX);
 
352
<a name="line353">353: </a>  <A href="../../../../docs/manualpages/DA/DAGetCorners.html#DAGetCorners">DAGetCorners</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,&amp;xs,&amp;ys,0,&amp;xm,&amp;ym,0);
 
353
<a name="line354">354: </a>  <A href="../../../../docs/manualpages/DA/DAVecGetArray.html#DAVecGetArray">DAVecGetArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,localX,&amp;x);
 
354
 
 
355
<a name="line356">356: </a>  <font color="#B22222">/* Evaluate Jacobian of function */</font>
 
356
<a name="line357">357: </a>  <font color="#4169E1">for</font> (j=ys; j&lt;ys+ym; j++) {
 
357
<a name="line358">358: </a>    <font color="#4169E1">for</font> (i=xs; i&lt;xs+xm; i++) {
 
358
<a name="line359">359: </a>      t0 = x[j][i];
 
359
 
 
360
<a name="line361">361: </a>      <font color="#4169E1">if</font> (i &gt; 0 &amp;&amp; i &lt; mx-1 &amp;&amp; j &gt; 0 &amp;&amp; j &lt; my-1) {
 
361
 
 
362
<a name="line363">363: </a>        <font color="#B22222">/* general interior volume */</font>
 
363
 
 
364
<a name="line365">365: </a>        tw = x[j][i-1];
 
365
<a name="line366">366: </a>        aw = 0.5*(t0 + tw);
 
366
<a name="line367">367: </a>        bw = PetscPowScalar(aw,bm1);
 
367
<a name="line368">368: </a>        <font color="#B22222">/* dw = bw * aw */</font>
 
368
<a name="line369">369: </a>        dw = PetscPowScalar(aw,beta);
 
369
<a name="line370">370: </a>        gw = coef*bw*(t0 - tw);
 
370
 
 
371
<a name="line372">372: </a>        te = x[j][i+1];
 
372
<a name="line373">373: </a>        ae = 0.5*(t0 + te);
 
373
<a name="line374">374: </a>        be = PetscPowScalar(ae,bm1);
 
374
<a name="line375">375: </a>        <font color="#B22222">/* de = be * ae; */</font>
 
375
<a name="line376">376: </a>        de = PetscPowScalar(ae,beta);
 
376
<a name="line377">377: </a>        ge = coef*be*(te - t0);
 
377
 
 
378
<a name="line379">379: </a>        ts = x[j-1][i];
 
379
<a name="line380">380: </a>        as = 0.5*(t0 + ts);
 
380
<a name="line381">381: </a>        bs = PetscPowScalar(as,bm1);
 
381
<a name="line382">382: </a>        <font color="#B22222">/* ds = bs * as; */</font>
 
382
<a name="line383">383: </a>        ds = PetscPowScalar(as,beta);
 
383
<a name="line384">384: </a>        gs = coef*bs*(t0 - ts);
 
384
<a name="line385">385: </a>
 
385
<a name="line386">386: </a>        tn = x[j+1][i];
 
386
<a name="line387">387: </a>        an = 0.5*(t0 + tn);
 
387
<a name="line388">388: </a>        bn = PetscPowScalar(an,bm1);
 
388
<a name="line389">389: </a>        <font color="#B22222">/* dn = bn * an; */</font>
 
389
<a name="line390">390: </a>        dn = PetscPowScalar(an,beta);
 
390
<a name="line391">391: </a>        gn = coef*bn*(tn - t0);
 
391
 
 
392
<a name="line393">393: </a>        v[0] = - hxdhy*(ds - gs);                                      col[0].j = j-1;       col[0].i = i;
 
393
<a name="line394">394: </a>        v[1] = - hydhx*(dw - gw);                                      col[1].j = j;         col[1].i = i-1;
 
394
<a name="line395">395: </a>        v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
 
395
<a name="line396">396: </a>        v[3] = - hydhx*(de + ge);                                      col[3].j = j;         col[3].i = i+1;
 
396
<a name="line397">397: </a>        v[4] = - hxdhy*(dn + gn);                                      col[4].j = j+1;       col[4].i = i;
 
397
<a name="line398">398: </a>        <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,5,col,v,INSERT_VALUES);
 
398
 
 
399
<a name="line400">400: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (i == 0) {
 
400
 
 
401
<a name="line402">402: </a>        <font color="#B22222">/* left-hand boundary */</font>
 
402
<a name="line403">403: </a>        tw = tleft;
 
403
<a name="line404">404: </a>        aw = 0.5*(t0 + tw);
 
404
<a name="line405">405: </a>        bw = PetscPowScalar(aw,bm1);
 
405
<a name="line406">406: </a>        <font color="#B22222">/* dw = bw * aw */</font>
 
406
<a name="line407">407: </a>        dw = PetscPowScalar(aw,beta);
 
407
<a name="line408">408: </a>        gw = coef*bw*(t0 - tw);
 
408
<a name="line409">409: </a>
 
409
<a name="line410">410: </a>        te = x[j][i + 1];
 
410
<a name="line411">411: </a>        ae = 0.5*(t0 + te);
 
411
<a name="line412">412: </a>        be = PetscPowScalar(ae,bm1);
 
412
<a name="line413">413: </a>        <font color="#B22222">/* de = be * ae; */</font>
 
413
<a name="line414">414: </a>        de = PetscPowScalar(ae,beta);
 
414
<a name="line415">415: </a>        ge = coef*be*(te - t0);
 
415
<a name="line416">416: </a>
 
416
<a name="line417">417: </a>        <font color="#B22222">/* left-hand bottom boundary */</font>
 
417
<a name="line418">418: </a>        <font color="#4169E1">if</font> (j == 0) {
 
418
 
 
419
<a name="line420">420: </a>          tn = x[j+1][i];
 
420
<a name="line421">421: </a>          an = 0.5*(t0 + tn);
 
421
<a name="line422">422: </a>          bn = PetscPowScalar(an,bm1);
 
422
<a name="line423">423: </a>          <font color="#B22222">/* dn = bn * an; */</font>
 
423
<a name="line424">424: </a>          dn = PetscPowScalar(an,beta);
 
424
<a name="line425">425: </a>          gn = coef*bn*(tn - t0);
 
425
<a name="line426">426: </a>
 
426
<a name="line427">427: </a>          v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;
 
427
<a name="line428">428: </a>          v[1] = - hydhx*(de + ge);                           col[1].j = j;         col[1].i = i+1;
 
428
<a name="line429">429: </a>          v[2] = - hxdhy*(dn + gn);                           col[2].j = j+1;       col[2].i = i;
 
429
<a name="line430">430: </a>          <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,3,col,v,INSERT_VALUES);
 
430
<a name="line431">431: </a>
 
431
<a name="line432">432: </a>        <font color="#B22222">/* left-hand interior boundary */</font>
 
432
<a name="line433">433: </a>        } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (j &lt; my-1) {
 
433
 
 
434
<a name="line435">435: </a>          ts = x[j-1][i];
 
435
<a name="line436">436: </a>          as = 0.5*(t0 + ts);
 
436
<a name="line437">437: </a>          bs = PetscPowScalar(as,bm1);
 
437
<a name="line438">438: </a>          <font color="#B22222">/* ds = bs * as; */</font>
 
438
<a name="line439">439: </a>          ds = PetscPowScalar(as,beta);
 
439
<a name="line440">440: </a>          gs = coef*bs*(t0 - ts);
 
440
<a name="line441">441: </a>
 
441
<a name="line442">442: </a>          tn = x[j+1][i];
 
442
<a name="line443">443: </a>          an = 0.5*(t0 + tn);
 
443
<a name="line444">444: </a>          bn = PetscPowScalar(an,bm1);
 
444
<a name="line445">445: </a>          <font color="#B22222">/* dn = bn * an; */</font>
 
445
<a name="line446">446: </a>          dn = PetscPowScalar(an,beta);
 
446
<a name="line447">447: </a>          gn = coef*bn*(tn - t0);
 
447
<a name="line448">448: </a>
 
448
<a name="line449">449: </a>          v[0] = - hxdhy*(ds - gs);                                      col[0].j = j-1;       col[0].i = i;
 
449
<a name="line450">450: </a>          v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
 
450
<a name="line451">451: </a>          v[2] = - hydhx*(de + ge);                                      col[2].j = j;         col[2].i = i+1;
 
451
<a name="line452">452: </a>          v[3] = - hxdhy*(dn + gn);                                      col[3].j = j+1;       col[3].i = i;
 
452
<a name="line453">453: </a>          <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,4,col,v,INSERT_VALUES);
 
453
<a name="line454">454: </a>        <font color="#B22222">/* left-hand top boundary */</font>
 
454
<a name="line455">455: </a>        } <font color="#4169E1">else</font> {
 
455
 
 
456
<a name="line457">457: </a>          ts = x[j-1][i];
 
457
<a name="line458">458: </a>          as = 0.5*(t0 + ts);
 
458
<a name="line459">459: </a>          bs = PetscPowScalar(as,bm1);
 
459
<a name="line460">460: </a>          <font color="#B22222">/* ds = bs * as; */</font>
 
460
<a name="line461">461: </a>          ds = PetscPowScalar(as,beta);
 
461
<a name="line462">462: </a>          gs = coef*bs*(t0 - ts);
470
462
<a name="line463">463: </a>
471
 
<a name="line464">464: </a>        <font color="#B22222">/* right-hand boundary */</font>
472
 
<a name="line465">465: </a>        tw = x[j][i-1];
473
 
<a name="line466">466: </a>        aw = 0.5*(t0 + tw);
474
 
<a name="line467">467: </a>        bw = PetscPowScalar(aw,bm1);
475
 
<a name="line468">468: </a>        <font color="#B22222">/* dw = bw * aw */</font>
476
 
<a name="line469">469: </a>        dw = PetscPowScalar(aw,beta);
477
 
<a name="line470">470: </a>        gw = coef*bw*(t0 - tw);
 
463
<a name="line464">464: </a>          v[0] = - hxdhy*(ds - gs);                            col[0].j = j-1;       col[0].i = i;
 
464
<a name="line465">465: </a>          v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[1].j = row.j = j; col[1].i = row.i = i;
 
465
<a name="line466">466: </a>          v[2] = - hydhx*(de + ge);                            col[2].j = j;         col[2].i = i+1;
 
466
<a name="line467">467: </a>          <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,3,col,v,INSERT_VALUES);
 
467
<a name="line468">468: </a>        }
 
468
 
 
469
<a name="line470">470: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (i == mx-1) {
478
470
<a name="line471">471: </a>
479
 
<a name="line472">472: </a>        te = tright;
480
 
<a name="line473">473: </a>        ae = 0.5*(t0 + te);
481
 
<a name="line474">474: </a>        be = PetscPowScalar(ae,bm1);
482
 
<a name="line475">475: </a>        <font color="#B22222">/* de = be * ae; */</font>
483
 
<a name="line476">476: </a>        de = PetscPowScalar(ae,beta);
484
 
<a name="line477">477: </a>        ge = coef*be*(te - t0);
485
 
<a name="line478">478: </a>
486
 
<a name="line479">479: </a>        <font color="#B22222">/* right-hand bottom boundary */</font>
487
 
<a name="line480">480: </a>        <font color="#4169E1">if</font> (j == 0) {
488
 
 
489
 
<a name="line482">482: </a>          tn = x[j+1][i];
490
 
<a name="line483">483: </a>          an = 0.5*(t0 + tn);
491
 
<a name="line484">484: </a>          bn = PetscPowScalar(an,bm1);
492
 
<a name="line485">485: </a>          <font color="#B22222">/* dn = bn * an; */</font>
493
 
<a name="line486">486: </a>          dn = PetscPowScalar(an,beta);
494
 
<a name="line487">487: </a>          gn = coef*bn*(tn - t0);
495
 
<a name="line488">488: </a>
496
 
<a name="line489">489: </a>          v[0] = - hydhx*(dw - gw);                           col[0].j = j;         col[0].i = i-1;
497
 
<a name="line490">490: </a>          v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
498
 
<a name="line491">491: </a>          v[2] = - hxdhy*(dn + gn);                           col[2].j = j+1;       col[2].i = i;
499
 
<a name="line492">492: </a>          <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,3,col,v,INSERT_VALUES);
500
 
<a name="line493">493: </a>
501
 
<a name="line494">494: </a>        <font color="#B22222">/* right-hand interior boundary */</font>
502
 
<a name="line495">495: </a>        } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (j &lt; my-1) {
503
 
 
504
 
<a name="line497">497: </a>          ts = x[j-1][i];
505
 
<a name="line498">498: </a>          as = 0.5*(t0 + ts);
506
 
<a name="line499">499: </a>          bs = PetscPowScalar(as,bm1);
507
 
<a name="line500">500: </a>          <font color="#B22222">/* ds = bs * as; */</font>
508
 
<a name="line501">501: </a>          ds = PetscPowScalar(as,beta);
509
 
<a name="line502">502: </a>          gs = coef*bs*(t0 - ts);
510
 
<a name="line503">503: </a>
511
 
<a name="line504">504: </a>          tn = x[j+1][i];
512
 
<a name="line505">505: </a>          an = 0.5*(t0 + tn);
513
 
<a name="line506">506: </a>          bn = PetscPowScalar(an,bm1);
514
 
<a name="line507">507: </a>          <font color="#B22222">/* dn = bn * an; */</font>
515
 
<a name="line508">508: </a>          dn = PetscPowScalar(an,beta);
516
 
<a name="line509">509: </a>          gn = coef*bn*(tn - t0);
517
 
<a name="line510">510: </a>
518
 
<a name="line511">511: </a>          v[0] = - hxdhy*(ds - gs);                                      col[0].j = j-1;       col[0].i = i;
519
 
<a name="line512">512: </a>          v[1] = - hydhx*(dw - gw);                                      col[1].j = j;         col[1].i = i-1;
520
 
<a name="line513">513: </a>          v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
521
 
<a name="line514">514: </a>          v[3] = - hxdhy*(dn + gn);                                      col[3].j = j+1;       col[3].i = i;
522
 
<a name="line515">515: </a>          <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,4,col,v,INSERT_VALUES);
523
 
<a name="line516">516: </a>        <font color="#B22222">/* right-hand top boundary */</font>
524
 
<a name="line517">517: </a>        } <font color="#4169E1">else</font> {
525
 
 
526
 
<a name="line519">519: </a>          ts = x[j-1][i];
527
 
<a name="line520">520: </a>          as = 0.5*(t0 + ts);
528
 
<a name="line521">521: </a>          bs = PetscPowScalar(as,bm1);
529
 
<a name="line522">522: </a>          <font color="#B22222">/* ds = bs * as; */</font>
530
 
<a name="line523">523: </a>          ds = PetscPowScalar(as,beta);
531
 
<a name="line524">524: </a>          gs = coef*bs*(t0 - ts);
532
 
<a name="line525">525: </a>
533
 
<a name="line526">526: </a>          v[0] = - hxdhy*(ds - gs);                            col[0].j = j-1;       col[0].i = i;
534
 
<a name="line527">527: </a>          v[1] = - hydhx*(dw - gw);                            col[1].j = j;         col[1].i = i-1;
535
 
<a name="line528">528: </a>          v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
536
 
<a name="line529">529: </a>          <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,3,col,v,INSERT_VALUES);
537
 
<a name="line530">530: </a>        }
538
 
 
539
 
<a name="line532">532: </a>      <font color="#B22222">/* bottom boundary,and i &lt;&gt; 0 or mx-1 */</font>
540
 
<a name="line533">533: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (j == 0) {
541
 
 
542
 
<a name="line535">535: </a>        tw = x[j][i-1];
543
 
<a name="line536">536: </a>        aw = 0.5*(t0 + tw);
544
 
<a name="line537">537: </a>        bw = PetscPowScalar(aw,bm1);
545
 
<a name="line538">538: </a>        <font color="#B22222">/* dw = bw * aw */</font>
546
 
<a name="line539">539: </a>        dw = PetscPowScalar(aw,beta);
547
 
<a name="line540">540: </a>        gw = coef*bw*(t0 - tw);
548
 
 
549
 
<a name="line542">542: </a>        te = x[j][i+1];
550
 
<a name="line543">543: </a>        ae = 0.5*(t0 + te);
551
 
<a name="line544">544: </a>        be = PetscPowScalar(ae,bm1);
552
 
<a name="line545">545: </a>        <font color="#B22222">/* de = be * ae; */</font>
553
 
<a name="line546">546: </a>        de = PetscPowScalar(ae,beta);
554
 
<a name="line547">547: </a>        ge = coef*be*(te - t0);
555
 
 
556
 
<a name="line549">549: </a>        tn = x[j+1][i];
557
 
<a name="line550">550: </a>        an = 0.5*(t0 + tn);
558
 
<a name="line551">551: </a>        bn = PetscPowScalar(an,bm1);
559
 
<a name="line552">552: </a>        <font color="#B22222">/* dn = bn * an; */</font>
560
 
<a name="line553">553: </a>        dn = PetscPowScalar(an,beta);
561
 
<a name="line554">554: </a>        gn = coef*bn*(tn - t0);
562
 
<a name="line555">555: </a>
563
 
<a name="line556">556: </a>        v[0] = - hydhx*(dw - gw);                           col[0].j = j;         col[0].i = i-1;
564
 
<a name="line557">557: </a>        v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
565
 
<a name="line558">558: </a>        v[2] = - hydhx*(de + ge);                           col[2].j = j;         col[2].i = i+1;
566
 
<a name="line559">559: </a>        v[3] = - hxdhy*(dn + gn);                           col[3].j = j+1;       col[3].i = i;
567
 
<a name="line560">560: </a>        <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,4,col,v,INSERT_VALUES);
568
 
<a name="line561">561: </a>
569
 
<a name="line562">562: </a>      <font color="#B22222">/* top boundary,and i &lt;&gt; 0 or mx-1 */</font>
570
 
<a name="line563">563: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (j == my-1) {
571
 
<a name="line564">564: </a>
572
 
<a name="line565">565: </a>        tw = x[j][i-1];
573
 
<a name="line566">566: </a>        aw = 0.5*(t0 + tw);
574
 
<a name="line567">567: </a>        bw = PetscPowScalar(aw,bm1);
575
 
<a name="line568">568: </a>        <font color="#B22222">/* dw = bw * aw */</font>
576
 
<a name="line569">569: </a>        dw = PetscPowScalar(aw,beta);
577
 
<a name="line570">570: </a>        gw = coef*bw*(t0 - tw);
578
 
 
579
 
<a name="line572">572: </a>        te = x[j][i+1];
580
 
<a name="line573">573: </a>        ae = 0.5*(t0 + te);
581
 
<a name="line574">574: </a>        be = PetscPowScalar(ae,bm1);
582
 
<a name="line575">575: </a>        <font color="#B22222">/* de = be * ae; */</font>
583
 
<a name="line576">576: </a>        de = PetscPowScalar(ae,beta);
584
 
<a name="line577">577: </a>        ge = coef*be*(te - t0);
585
 
 
586
 
<a name="line579">579: </a>        ts = x[j-1][i];
587
 
<a name="line580">580: </a>        as = 0.5*(t0 + ts);
588
 
<a name="line581">581: </a>        bs = PetscPowScalar(as,bm1);
589
 
<a name="line582">582: </a>         <font color="#B22222">/* ds = bs * as; */</font>
590
 
<a name="line583">583: </a>        ds = PetscPowScalar(as,beta);
591
 
<a name="line584">584: </a>        gs = coef*bs*(t0 - ts);
592
 
 
593
 
<a name="line586">586: </a>        v[0] = - hxdhy*(ds - gs);                            col[0].j = j-1;       col[0].i = i;
594
 
<a name="line587">587: </a>        v[1] = - hydhx*(dw - gw);                            col[1].j = j;         col[1].i = i-1;
595
 
<a name="line588">588: </a>        v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
596
 
<a name="line589">589: </a>        v[3] = - hydhx*(de + ge);                            col[3].j = j;         col[3].i = i+1;
597
 
<a name="line590">590: </a>        <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,4,col,v,INSERT_VALUES);
598
 
<a name="line591">591: </a>
599
 
<a name="line592">592: </a>      }
600
 
<a name="line593">593: </a>    }
601
 
<a name="line594">594: </a>  }
602
 
<a name="line595">595: </a>  <A href="../../../../docs/manualpages/Mat/MatAssemblyBegin.html#MatAssemblyBegin">MatAssemblyBegin</A>(jac,MAT_FINAL_ASSEMBLY);
603
 
<a name="line596">596: </a>  <A href="../../../../docs/manualpages/DA/DAVecRestoreArray.html#DAVecRestoreArray">DAVecRestoreArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,localX,(void**)&amp;x);
604
 
<a name="line597">597: </a>  <A href="../../../../docs/manualpages/Mat/MatAssemblyEnd.html#MatAssemblyEnd">MatAssemblyEnd</A>(jac,MAT_FINAL_ASSEMBLY);
605
 
<a name="line598">598: </a>  <A href="../../../../docs/manualpages/DA/DARestoreLocalVector.html#DARestoreLocalVector">DARestoreLocalVector</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,&amp;localX);
606
 
 
607
 
<a name="line600">600: </a>  <A href="../../../../docs/manualpages/Profiling/PetscLogFlops.html#PetscLogFlops">PetscLogFlops</A>((41 + 8*POWFLOP)*xm*ym);
608
 
<a name="line601">601: </a>  <font color="#4169E1">return</font>(0);
609
 
<a name="line602">602: </a>}
 
471
<a name="line472">472: </a>        <font color="#B22222">/* right-hand boundary */</font>
 
472
<a name="line473">473: </a>        tw = x[j][i-1];
 
473
<a name="line474">474: </a>        aw = 0.5*(t0 + tw);
 
474
<a name="line475">475: </a>        bw = PetscPowScalar(aw,bm1);
 
475
<a name="line476">476: </a>        <font color="#B22222">/* dw = bw * aw */</font>
 
476
<a name="line477">477: </a>        dw = PetscPowScalar(aw,beta);
 
477
<a name="line478">478: </a>        gw = coef*bw*(t0 - tw);
 
478
<a name="line479">479: </a>
 
479
<a name="line480">480: </a>        te = tright;
 
480
<a name="line481">481: </a>        ae = 0.5*(t0 + te);
 
481
<a name="line482">482: </a>        be = PetscPowScalar(ae,bm1);
 
482
<a name="line483">483: </a>        <font color="#B22222">/* de = be * ae; */</font>
 
483
<a name="line484">484: </a>        de = PetscPowScalar(ae,beta);
 
484
<a name="line485">485: </a>        ge = coef*be*(te - t0);
 
485
<a name="line486">486: </a>
 
486
<a name="line487">487: </a>        <font color="#B22222">/* right-hand bottom boundary */</font>
 
487
<a name="line488">488: </a>        <font color="#4169E1">if</font> (j == 0) {
 
488
 
 
489
<a name="line490">490: </a>          tn = x[j+1][i];
 
490
<a name="line491">491: </a>          an = 0.5*(t0 + tn);
 
491
<a name="line492">492: </a>          bn = PetscPowScalar(an,bm1);
 
492
<a name="line493">493: </a>          <font color="#B22222">/* dn = bn * an; */</font>
 
493
<a name="line494">494: </a>          dn = PetscPowScalar(an,beta);
 
494
<a name="line495">495: </a>          gn = coef*bn*(tn - t0);
 
495
<a name="line496">496: </a>
 
496
<a name="line497">497: </a>          v[0] = - hydhx*(dw - gw);                           col[0].j = j;         col[0].i = i-1;
 
497
<a name="line498">498: </a>          v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
 
498
<a name="line499">499: </a>          v[2] = - hxdhy*(dn + gn);                           col[2].j = j+1;       col[2].i = i;
 
499
<a name="line500">500: </a>          <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,3,col,v,INSERT_VALUES);
 
500
<a name="line501">501: </a>
 
501
<a name="line502">502: </a>        <font color="#B22222">/* right-hand interior boundary */</font>
 
502
<a name="line503">503: </a>        } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (j &lt; my-1) {
 
503
 
 
504
<a name="line505">505: </a>          ts = x[j-1][i];
 
505
<a name="line506">506: </a>          as = 0.5*(t0 + ts);
 
506
<a name="line507">507: </a>          bs = PetscPowScalar(as,bm1);
 
507
<a name="line508">508: </a>          <font color="#B22222">/* ds = bs * as; */</font>
 
508
<a name="line509">509: </a>          ds = PetscPowScalar(as,beta);
 
509
<a name="line510">510: </a>          gs = coef*bs*(t0 - ts);
 
510
<a name="line511">511: </a>
 
511
<a name="line512">512: </a>          tn = x[j+1][i];
 
512
<a name="line513">513: </a>          an = 0.5*(t0 + tn);
 
513
<a name="line514">514: </a>          bn = PetscPowScalar(an,bm1);
 
514
<a name="line515">515: </a>          <font color="#B22222">/* dn = bn * an; */</font>
 
515
<a name="line516">516: </a>          dn = PetscPowScalar(an,beta);
 
516
<a name="line517">517: </a>          gn = coef*bn*(tn - t0);
 
517
<a name="line518">518: </a>
 
518
<a name="line519">519: </a>          v[0] = - hxdhy*(ds - gs);                                      col[0].j = j-1;       col[0].i = i;
 
519
<a name="line520">520: </a>          v[1] = - hydhx*(dw - gw);                                      col[1].j = j;         col[1].i = i-1;
 
520
<a name="line521">521: </a>          v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
 
521
<a name="line522">522: </a>          v[3] = - hxdhy*(dn + gn);                                      col[3].j = j+1;       col[3].i = i;
 
522
<a name="line523">523: </a>          <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,4,col,v,INSERT_VALUES);
 
523
<a name="line524">524: </a>        <font color="#B22222">/* right-hand top boundary */</font>
 
524
<a name="line525">525: </a>        } <font color="#4169E1">else</font> {
 
525
 
 
526
<a name="line527">527: </a>          ts = x[j-1][i];
 
527
<a name="line528">528: </a>          as = 0.5*(t0 + ts);
 
528
<a name="line529">529: </a>          bs = PetscPowScalar(as,bm1);
 
529
<a name="line530">530: </a>          <font color="#B22222">/* ds = bs * as; */</font>
 
530
<a name="line531">531: </a>          ds = PetscPowScalar(as,beta);
 
531
<a name="line532">532: </a>          gs = coef*bs*(t0 - ts);
 
532
<a name="line533">533: </a>
 
533
<a name="line534">534: </a>          v[0] = - hxdhy*(ds - gs);                            col[0].j = j-1;       col[0].i = i;
 
534
<a name="line535">535: </a>          v[1] = - hydhx*(dw - gw);                            col[1].j = j;         col[1].i = i-1;
 
535
<a name="line536">536: </a>          v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
 
536
<a name="line537">537: </a>          <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,3,col,v,INSERT_VALUES);
 
537
<a name="line538">538: </a>        }
 
538
 
 
539
<a name="line540">540: </a>      <font color="#B22222">/* bottom boundary,and i &lt;&gt; 0 or mx-1 */</font>
 
540
<a name="line541">541: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (j == 0) {
 
541
 
 
542
<a name="line543">543: </a>        tw = x[j][i-1];
 
543
<a name="line544">544: </a>        aw = 0.5*(t0 + tw);
 
544
<a name="line545">545: </a>        bw = PetscPowScalar(aw,bm1);
 
545
<a name="line546">546: </a>        <font color="#B22222">/* dw = bw * aw */</font>
 
546
<a name="line547">547: </a>        dw = PetscPowScalar(aw,beta);
 
547
<a name="line548">548: </a>        gw = coef*bw*(t0 - tw);
 
548
 
 
549
<a name="line550">550: </a>        te = x[j][i+1];
 
550
<a name="line551">551: </a>        ae = 0.5*(t0 + te);
 
551
<a name="line552">552: </a>        be = PetscPowScalar(ae,bm1);
 
552
<a name="line553">553: </a>        <font color="#B22222">/* de = be * ae; */</font>
 
553
<a name="line554">554: </a>        de = PetscPowScalar(ae,beta);
 
554
<a name="line555">555: </a>        ge = coef*be*(te - t0);
 
555
 
 
556
<a name="line557">557: </a>        tn = x[j+1][i];
 
557
<a name="line558">558: </a>        an = 0.5*(t0 + tn);
 
558
<a name="line559">559: </a>        bn = PetscPowScalar(an,bm1);
 
559
<a name="line560">560: </a>        <font color="#B22222">/* dn = bn * an; */</font>
 
560
<a name="line561">561: </a>        dn = PetscPowScalar(an,beta);
 
561
<a name="line562">562: </a>        gn = coef*bn*(tn - t0);
 
562
<a name="line563">563: </a>
 
563
<a name="line564">564: </a>        v[0] = - hydhx*(dw - gw);                           col[0].j = j;         col[0].i = i-1;
 
564
<a name="line565">565: </a>        v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
 
565
<a name="line566">566: </a>        v[2] = - hydhx*(de + ge);                           col[2].j = j;         col[2].i = i+1;
 
566
<a name="line567">567: </a>        v[3] = - hxdhy*(dn + gn);                           col[3].j = j+1;       col[3].i = i;
 
567
<a name="line568">568: </a>        <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,4,col,v,INSERT_VALUES);
 
568
<a name="line569">569: </a>
 
569
<a name="line570">570: </a>      <font color="#B22222">/* top boundary,and i &lt;&gt; 0 or mx-1 */</font>
 
570
<a name="line571">571: </a>      } <font color="#4169E1">else</font> <font color="#4169E1">if</font> (j == my-1) {
 
571
<a name="line572">572: </a>
 
572
<a name="line573">573: </a>        tw = x[j][i-1];
 
573
<a name="line574">574: </a>        aw = 0.5*(t0 + tw);
 
574
<a name="line575">575: </a>        bw = PetscPowScalar(aw,bm1);
 
575
<a name="line576">576: </a>        <font color="#B22222">/* dw = bw * aw */</font>
 
576
<a name="line577">577: </a>        dw = PetscPowScalar(aw,beta);
 
577
<a name="line578">578: </a>        gw = coef*bw*(t0 - tw);
 
578
 
 
579
<a name="line580">580: </a>        te = x[j][i+1];
 
580
<a name="line581">581: </a>        ae = 0.5*(t0 + te);
 
581
<a name="line582">582: </a>        be = PetscPowScalar(ae,bm1);
 
582
<a name="line583">583: </a>        <font color="#B22222">/* de = be * ae; */</font>
 
583
<a name="line584">584: </a>        de = PetscPowScalar(ae,beta);
 
584
<a name="line585">585: </a>        ge = coef*be*(te - t0);
 
585
 
 
586
<a name="line587">587: </a>        ts = x[j-1][i];
 
587
<a name="line588">588: </a>        as = 0.5*(t0 + ts);
 
588
<a name="line589">589: </a>        bs = PetscPowScalar(as,bm1);
 
589
<a name="line590">590: </a>         <font color="#B22222">/* ds = bs * as; */</font>
 
590
<a name="line591">591: </a>        ds = PetscPowScalar(as,beta);
 
591
<a name="line592">592: </a>        gs = coef*bs*(t0 - ts);
 
592
 
 
593
<a name="line594">594: </a>        v[0] = - hxdhy*(ds - gs);                            col[0].j = j-1;       col[0].i = i;
 
594
<a name="line595">595: </a>        v[1] = - hydhx*(dw - gw);                            col[1].j = j;         col[1].i = i-1;
 
595
<a name="line596">596: </a>        v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge);  col[2].j = row.j = j; col[2].i = row.i = i;
 
596
<a name="line597">597: </a>        v[3] = - hydhx*(de + ge);                            col[3].j = j;         col[3].i = i+1;
 
597
<a name="line598">598: </a>        <A href="../../../../docs/manualpages/Mat/MatSetValuesStencil.html#MatSetValuesStencil">MatSetValuesStencil</A>(jac,1,&amp;row,4,col,v,INSERT_VALUES);
 
598
<a name="line599">599: </a>
 
599
<a name="line600">600: </a>      }
 
600
<a name="line601">601: </a>    }
 
601
<a name="line602">602: </a>  }
 
602
<a name="line603">603: </a>  <A href="../../../../docs/manualpages/Mat/MatAssemblyBegin.html#MatAssemblyBegin">MatAssemblyBegin</A>(jac,MAT_FINAL_ASSEMBLY);
 
603
<a name="line604">604: </a>  <A href="../../../../docs/manualpages/DA/DAVecRestoreArray.html#DAVecRestoreArray">DAVecRestoreArray</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,localX,&amp;x);
 
604
<a name="line605">605: </a>  <A href="../../../../docs/manualpages/Mat/MatAssemblyEnd.html#MatAssemblyEnd">MatAssemblyEnd</A>(jac,MAT_FINAL_ASSEMBLY);
 
605
<a name="line606">606: </a>  <A href="../../../../docs/manualpages/DA/DARestoreLocalVector.html#DARestoreLocalVector">DARestoreLocalVector</A>((<A href="../../../../docs/manualpages/DA/DA.html#DA">DA</A>)dmmg-&gt;dm,&amp;localX);
 
606
 
 
607
<a name="line608">608: </a>  <A href="../../../../docs/manualpages/Profiling/PetscLogFlops.html#PetscLogFlops">PetscLogFlops</A>((41 + 8*POWFLOP)*xm*ym);
 
608
<a name="line609">609: </a>  <font color="#4169E1">return</font>(0);
 
609
<a name="line610">610: </a>}
610
610
 
611
611
</pre>
612
612
</body>