[petsc-users] Jacobian matrix construction when dof > 1
Smith, Barry F.
bsmith at mcs.anl.gov
Thu Jan 25 10:08:14 CST 2018
> 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.
More information about the petsc-users
mailing list