[petsc-users] Jacobian matrix construction when dof > 1

[Smith, Barry F.]

> On Jan 25, 2018, at 8:38 AM, Rahul Samala <srahul_05 at yahoo.co.in> wrote:
> 
> Hello PetSc users,
> 
> I am a beginner to PetSc and I am trying to write Jacobian matrix for my problem involving 2 nonlinear equations in 2 unknown variables, solving in a 1D domain.
> Most of the petsc examples involve one variable and I am finding it difficult to understand how the Jacobian matrix is built when there are more than one variables.
> For my problem the DOF is 2 (since 2 unknown variables). I am using finite difference and my stencil is 3 (i-1,i,i+1).
> Jacobian matrix will have 2x2 submatrix elements and with stencil 3, we have to set 4x3=12 values for each grid point.
> Let the unknown variables are "u,v" and (F1_u)_i-1 represent "partial derivative of residual of equation 1 at grid point 'i' with respect to variable 'u' at i-1".
> 
> Is the following pseudo-code the correct way of filling jacobian matrix?
> 
> FormJacobianLocal(DMDALocalInfo *info,Field *sol,Mat J,Mat Jpre, AppCtx *user)
> {
>   .......
>     MatStencil     col[2][3],row;
>     PetscReal       v[2][3];
> .........
> for(i=xs;i<xs+xm;i++){
>     row.i = i; row.c=0;
>     { // Internal points
>         col[0][0].i = i-1;     col[0][0].c = 0;     v[0][0] = jacobian value =  (F1_u)_i-1
>         col[0][1].i = i;        col[0][1].c = 0;     v[0][1] = jacobian value = (F1_u)_i
>         col[0][2].i = i+1;    col[0][2].c = 0;    v[0][1] = jacobian value =  (F1_u)_i+1

   Final one should be v[0][2] same below
>        
>         col[1][0].i = i-1;     col[1][0].c = 1;     v[1][0] = jacobian value =  (F1_v)_i-1
>         col[1][1].i = i;        col[1][1].c = 1;     v[1][1] = jacobian value = (F1_v)_i
>         col[1][2].i = i+1;    col[1][2].c = 1;    v[1][1] = jacobian value =  (F1_v)_i+1
>     }
>   MatSetValuesStencil(Jpre,1,&row,6,col,v,INSERT_VALUES);
> row.i = i; row.c=1;
>     { // Internal points
>         col[0][0].i = i-1;     col[0][0].c = 0;     v[0][0] = jacobian value =  (F2_u)_i-1
>         col[0][1].i = i;        col[0][1].c = 0;     v[0][1] = jacobian value = (F2_u)_i
>         col[0][2].i = i+1;    col[0][2].c = 0;    v[0][1] = jacobian value =  (F2_u)_i+1
>        
>         col[1][0].i = i-1;     col[1][0].c = 1;     v[1][0] = jacobian value =  (F2_v)_i-1
>         col[1][1].i = i;        col[1][1].c = 1;     v[1][1] = jacobian value = (F2_v)_i
>         col[1][2].i = i+1;    col[1][2].c = 1;    v[1][1] = jacobian value =  (F2_v)_i+1
>     }
>   MatSetValuesStencil(Jpre,1,&row,6,col,v,INSERT_VALUES);
> }

   Note you can run with -mat_view for a problem with say 4 grid points and then manually check that the entires are correct, I usually do this when implementing a new derivative.

  Barry

> .......
> }
> 
> Thank you,
> Rahul.