static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";

#include <petscdmda.h>
#include <petscksp.h>
#include <petsctime.h>

typedef struct {
    PassiveReal mu;         // viscosity
    Vec         phi;        // level set function
} Stokesparas;

typedef struct {
    PassiveReal x_min,x_max,y_min,y_max;
} resolpara;

typedef struct {
    PetscScalar u,v,p;
} Field;

typedef struct {
    KSP         ksp;        // linear Krylov solver
    Stokesparas para;
} Stokes2D;

extern PetscErrorCode Stokes2D_initial(Stokes2D *s, DM da, resolpara *rp);
extern PetscErrorCode Stokes2D_solve(Stokes2D *s);
extern PetscErrorCode Stokes2D_finalize(Stokes2D *s);
PetscErrorCode printSol(Stokes2D* s, resolpara *rp);
extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,MatStructure*,void*);
extern PetscErrorCode ComputeRHS(KSP,Vec,void*);

int main(int argc,char **argv)
{
    Stokes2D       stokes;
    DM             da;
    PetscErrorCode ierr;
    PetscReal       t0,t1,t2,t3;
    resolpara       rp;
    
    ierr = PetscTime(&t0);CHKERRQ(ierr);
    PetscInitialize(&argc,&argv,(char*)0,help);
    ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-3,-3,PETSC_DECIDE,PETSC_DECIDE,3,3,0,0,&da);CHKERRQ(ierr);
    Stokes2D_initial(&stokes, da, &rp);
    ierr = PetscTime(&t1);CHKERRQ(ierr);
    Stokes2D_solve(&stokes);
    ierr = PetscTime(&t2);CHKERRQ(ierr);
    printSol(&stokes, &rp);
    ierr = PetscTime(&t3);CHKERRQ(ierr);
    
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Time cost for creating solver context %g s, and for solving %g s, and for printing %g s.\n",t1-t0,t2-t1,t3-t2);CHKERRQ(ierr);
    ierr = Stokes2D_finalize(&stokes);CHKERRQ(ierr);
    ierr = PetscFinalize();CHKERRQ(ierr);
    
    return 0;
}

PetscScalar phiEvalAt(PetscScalar x, PetscScalar y)
{
    return PetscSqrtScalar(PetscPowScalar(x-.5,2.0)+PetscPowScalar(y-.5,2.0))-.2;
}

PetscScalar phixEvalAt(PetscScalar x, PetscScalar y)
{
    PetscScalar h = 1e-6;
    return (phiEvalAt(x+h/2,y)-phiEvalAt(x-h/2,y))/h;
}

PetscScalar phiyEvalAt(PetscScalar x, PetscScalar y)
{
    PetscScalar h = 1e-6;
    return (phiEvalAt(x,y+h/2)-phiEvalAt(x,y-h/2))/h;
}

PetscErrorCode normal_phiEvalAt(PetscScalar x, PetscScalar y, PetscScalar *n1, PetscScalar *n2)
{
    PetscScalar px = phixEvalAt(x,y), py = phiyEvalAt(x,y);
    
    PetscFunctionBeginUser;
    *n1 = px/PetscSqrtScalar(px*px+py*py);
    *n2 = py/PetscSqrtScalar(px*px+py*py);
    PetscFunctionReturn(0);
}

PetscScalar u_exact(PetscScalar x, PetscScalar y)
{
    return PetscSinScalar(2*PETSC_PI*x)*PetscCosScalar(2*PETSC_PI*y);
}

PetscScalar v_exact(PetscScalar x, PetscScalar y)
{
    return -PetscCosScalar(2*PETSC_PI*x)*PetscSinScalar(2*PETSC_PI*y);
}

PetscScalar p_exact(PetscScalar x, PetscScalar y)
{
    return -PetscSinScalar(2*PETSC_PI*x)*PetscSinScalar(2*PETSC_PI*y);
}

PetscScalar f1_exact(PetscScalar x, PetscScalar y)
{
    PetscScalar h = 1e-6;
    return -(u_exact(x+h,y)+u_exact(x-h,y)+u_exact(x,y+h)+u_exact(x,y-h)-4*u_exact(x,y))/h/h+(p_exact(x+h/2.0,y)-p_exact(x-h/2.0,y))/h;
}

PetscScalar f2_exact(PetscScalar x, PetscScalar y)
{
    PetscScalar h = 1e-6;
    return -(v_exact(x+h,y)+v_exact(x-h,y)+v_exact(x,y+h)+v_exact(x,y-h)-4*v_exact(x,y))/h/h+(p_exact(x,y+h/2.0)-p_exact(x,y-h/2.0))/h;
}

PetscErrorCode bdry_force(PetscScalar x, PetscScalar y, PetscScalar n1, PetscScalar n2, PetscScalar *g1, PetscScalar *g2)
{
    PetscScalar h = 1e-6, ux, uy, vx, vy;
    
    PetscFunctionBeginUser;
    ux  = (u_exact(x+h/2,y)-u_exact(x-h/2,y))/h;
    uy  = (u_exact(x,y+h/2)-u_exact(x,y-h/2))/h;
    vx  = (v_exact(x+h/2,y)-v_exact(x-h/2,y))/h;
    vy  = (v_exact(x,y+h/2)-v_exact(x,y-h/2))/h;
    *g1 = (n1*n1-n2*n2)*(uy+vx)+2*n1*n2*(vy-ux);
    *g2 = 2*(n1*n1-n2*n2)*ux+2*n1*n2*(uy+vx)-p_exact(x,y);
    PetscFunctionReturn(0);
}

PetscErrorCode x_exact(KSP ksp,Vec x, resolpara *rp)
{
    PetscErrorCode ierr;
    PetscInt       i,j,mx,my,xm,ym,xs,ys;
    PetscScalar    hx,hy;
    Field          **array;
    DM             da;
    
    PetscFunctionBeginUser;
    ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
    ierr = DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
    hx   = (rp->x_max-rp->x_min) / (PetscReal)(mx-1);
    hy   = (rp->y_max-rp->y_min) / (PetscReal)(my-1);
    ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
    ierr = DMDAVecGetArray(da, x, &array);CHKERRQ(ierr);
    for (j=ys; j<ys+ym; j++) {
        for (i=xs; i<xs+xm; i++) {
            array[j][i].u = u_exact(i*hx,(j+.5)*hy);
            array[j][i].v = v_exact((i+.5)*hx,j*hy);
            array[j][i].p = p_exact((i+.5)*hx,(j+.5)*hy);
        }
    }
    ierr = DMDAVecRestoreArray(da, x, &array);CHKERRQ(ierr);
    ierr = VecAssemblyBegin(x);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(x);CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

PetscErrorCode Stokes2D_initial(Stokes2D *s, DM da, resolpara *rp)
{
    PetscErrorCode ierr;
    
    PetscFunctionBeginUser;
    rp->x_min = 0;
    rp->x_max = 1;
    rp->y_min = 0;
    rp->y_max = 1;
    ierr = PetscOptionsGetReal(NULL,"-x_min",&(rp->x_min),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL,"-x_max",&(rp->x_max),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL,"-y_min",&(rp->y_min),NULL);CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(NULL,"-y_max",&(rp->y_max),NULL);CHKERRQ(ierr);
    ierr = KSPCreate(PETSC_COMM_WORLD,&(s->ksp));CHKERRQ(ierr);
    ierr = KSPSetComputeRHS(s->ksp,ComputeRHS,rp);CHKERRQ(ierr);
    ierr = KSPSetComputeOperators(s->ksp,ComputeMatrix,rp);CHKERRQ(ierr);
    ierr = KSPSetDM(s->ksp,da);CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

PetscErrorCode Stokes2D_solve(Stokes2D *s)
{
    PetscErrorCode ierr;
    
    PetscFunctionBeginUser;
    ierr = KSPSetFromOptions(s->ksp);CHKERRQ(ierr);
    ierr = KSPSetUp(s->ksp);CHKERRQ(ierr);
    ierr = KSPSolve(s->ksp,NULL,NULL);CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

PetscErrorCode Stokes2D_finalize(Stokes2D *s)
{
    PetscErrorCode ierr;
    
    PetscFunctionBeginUser;
    ierr = KSPDestroy(&(s->ksp));CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

PetscErrorCode printSol(Stokes2D* s, resolpara *rp)
{
    PetscErrorCode ierr;
    PetscViewer    viewer;
    Vec            b,x,zero;
    Mat            A;
    PetscScalar    r,hx,hy;
    PetscInt       mx, my;
    DM             da;
    
    PetscFunctionBeginUser;
    ierr = KSPGetSolution(s->ksp,&x);CHKERRQ(ierr);
    ierr = KSPGetRhs(s->ksp,&b);CHKERRQ(ierr);
    ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "output/ksp_x.txt", &viewer);CHKERRQ(ierr);
    ierr = VecView(x, viewer);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
    ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "output/ksp_b.txt", &viewer);CHKERRQ(ierr);
    ierr = VecView(b, viewer);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
    
    ierr = x_exact(s->ksp, b, rp);CHKERRQ(ierr);
    ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "output/ksp_sol.txt", &viewer);CHKERRQ(ierr);
    ierr = VecView(b, viewer);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
    
    ierr = KSPGetOperators(s->ksp,&A,NULL,NULL);CHKERRQ(ierr);
    ierr = VecDuplicate(b,&zero);CHKERRQ(ierr);
    ierr = MatMult(A,b,zero);
    ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "output/ksp_zero.txt", &viewer);CHKERRQ(ierr);
    ierr = VecView(zero, viewer);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
    ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "output/ksp_A.txt", &viewer);CHKERRQ(ierr);
    ierr = MatView(A, viewer);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
    
    VecAXPY(x,-1,b);
    VecNorm(x,NORM_2,&r);
    ierr = KSPGetDM(s->ksp,&da);CHKERRQ(ierr);
    ierr = DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
    hx   = (rp->x_max-rp->x_min) / (PetscReal)(mx-1);
    hy   = (rp->y_max-rp->y_min) / (PetscReal)(my-1);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm error %g %g.\n",r,r*PetscSqrtScalar(hx*hy));CHKERRQ(ierr);
    
    PetscFunctionReturn(0);
}

PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
{
    PetscErrorCode ierr;
    PetscInt       i,j,mx,my,xm,ym,xs,ys;
    PetscScalar    hx,hy;
    Field          **array;
    DM             da;
    resolpara      *rp;
    
    PetscFunctionBeginUser;
    rp   = ctx;
    ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
    ierr = DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
    hx   = (rp->x_max-rp->x_min) / (PetscReal)(mx-1);
    hy   = (rp->y_max-rp->y_min) / (PetscReal)(my-1);
    ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
    ierr = DMDAVecGetArray(da, b, &array);CHKERRQ(ierr);
    for (j=ys; j<ys+ym; j++) {
        for (i=xs; i<xs+xm; i++) {
            if (i==mx-1 || j==my-1) {
                array[j][i].u = u_exact(i*hx,(j+.5)*hy);
                array[j][i].v = v_exact((i+.5)*hx,j*hy);
                array[j][i].p = p_exact((i+.5)*hx,(j+.5)*hy);
            } else if (i==0) {
                if (j==my-1) {
                    array[j][i].u = u_exact(i*hx,(j+.5)*hy);
                    array[j][i].v = v_exact((i+.5)*hx,j*hy);
                    array[j][i].p = p_exact((i+.5)*hx,(j+.5)*hy);
                } else if (j==0) {
                    array[j][i].u = u_exact(i*hx,(j+.5)*hy);
                    array[j][i].v = v_exact((i+.5)*hx,j*hy);
                    array[j][i].p = p_exact((i+.5)*hx,(j+.5)*hy)*hx;
                } else {
                    array[j][i].u = u_exact(i*hx,(j+.5)*hy);
                    array[j][i].v = (- v_exact((i-.5)*hx,j*hy)/hx/hx - f2_exact((i+.5)*hx,j*hy))*hy*hy;
                    array[j][i].p = 0;
                }
            } else if (j==0) {
                if (i==mx-1) {
                    array[j][i].u = u_exact(i*hx,(j+.5)*hy);
                    array[j][i].v = v_exact((i+.5)*hx,j*hy);
                    array[j][i].p = p_exact((i+.5)*hx,(j+.5)*hy);
                } else if (i==0) {
                    array[j][i].u = u_exact(i*hx,(j+.5)*hy);
                    array[j][i].v = v_exact((i+.5)*hx,j*hy);
                    array[j][i].p = p_exact((i+.5)*hx,(j+.5)*hy)*hx;
                } else {
                    array[j][i].u = (- u_exact(i*hx,(j-.5)*hy)/hy/hy - f1_exact(i*hx,(j+.5)*hy))*hx*hx;
                    array[j][i].v = v_exact((i+.5)*hx,j*hy);
                    array[j][i].p = 0;
                }
            } else {
                array[j][i].u = -f1_exact(i*hx,(j+.5)*hy)*hx*hx;
                array[j][i].v = -f2_exact((i+.5)*hx,j*hy)*hy*hy;
                array[j][i].p = 0;
            }
        }
    }
    ierr = DMDAVecRestoreArray(da, b, &array);CHKERRQ(ierr);
    ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,MatStructure *str,void *ctx)
{
    PetscErrorCode ierr;
    PetscInt       i,j,k,mx,my,xm,ym,xs,ys,s1,s2,s3,s4,s5,s6,s7;
    PetscScalar    v[15], alpha[6], beta[6], gamma[6], eta[6], xi[3], a[12], b[15];
    PetscScalar    hx, hy, dhx, dhy, dhx2, dhy2, hxdhy, hydhx;
    PetscScalar    y0, y1, y2, y3, y4, x_star, y_star, n1, n2, l1, l2;
    MatStencil     row, col[15];
    DM             da;
    resolpara      *rp;
    PetscInt       choice, num_points;
    
    PetscFunctionBeginUser;
    rp        = ctx;
    ierr      = KSPGetDM(ksp,&da);CHKERRQ(ierr);
    ierr      = DMDAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
    hx        = (rp->x_max-rp->x_min) / (PetscReal)(mx-1);
    hy        = (rp->y_max-rp->y_min) / (PetscReal)(my-1);
    dhx       = 1/hx;
    dhy       = 1/hy;
    dhx2      = 1/hx/hx;
    dhy2      = 1/hy/hy;
    hxdhy     = hx/hy;
    hydhx     = hy/hx;
    ierr      = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
    for (j=ys; j<ys+ym; j++) {
        for (i=xs; i<xs+xm; i++) {
            row.i = i; row.j = j;
            if (i==mx-1 || j==my-1) {
                // on the right or top bdry
                v[0]=1;   row.c = 0;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
                v[0]=1;   row.c = 1;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
                v[0]=1;   row.c = 2;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
            } else if (i==0) {
                // on the left bdry
                if (j==my-1) {
                    v[0]=1;   row.c = 0;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
                    v[0]=1;   row.c = 1;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
                    v[0]=1;   row.c = 2;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
                } else if (j==0) {
                    v[0]=1;   row.c = 0;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
                    v[0]=1;   row.c = 1;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
                    row.c = 2;
                    v[0] = -1;              col[0].i = i;   col[0].j = j;   col[0].c = 0;
                    v[1] = 1;               col[1].i = i+1; col[1].j = j;   col[1].c = 0;
                    v[2] = -hxdhy;          col[2].i = i;   col[2].j = j;   col[2].c = 1;
                    v[3] = hxdhy;           col[3].i = i;   col[3].j = j+1; col[3].c = 1;
                    v[4] = hx;              col[4].i = i;   col[4].j = j;   col[4].c = 2;
                    ierr = MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr);
                } else {
                    v[0]=1;   row.c = 0;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
                    row.c = 1;
                    v[0] = hydhx*hydhx;         col[0].i = i+1; col[0].j = j;   col[0].c = 1;
                    v[1] = (-2*hydhx*hydhx-2);  col[1].i = i;   col[1].j = j;   col[1].c = 1;
                    v[2] = 1;                   col[2].i = i;   col[2].j = j-1; col[2].c = 1;
                    v[3] = 1;                   col[3].i = i;   col[3].j = j+1; col[3].c = 1;
                    v[4] = hy;                  col[4].i = i;   col[4].j = j-1; col[4].c = 2;
                    v[5] = -hy;                 col[5].i = i;   col[5].j = j;   col[5].c = 2;
                    ierr = MatSetValuesStencil(jac,1,&row,6,col,v,INSERT_VALUES);CHKERRQ(ierr);
                    row.c = 2;
                    v[0] = -hydhx;          col[0].i = i;   col[0].j = j;   col[0].c = 0;
                    v[1] = hydhx;           col[1].i = i+1; col[1].j = j;   col[1].c = 0;
                    v[2] = -1;              col[2].i = i;   col[2].j = j;   col[2].c = 1;
                    v[3] = 1;               col[3].i = i;   col[3].j = j+1; col[3].c = 1;
                    ierr = MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);CHKERRQ(ierr);
                }
            } else if (j==0) {
                // on the bottom bdry
                if (i==mx-1) {
                    v[0]=1;   row.c = 0;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
                    v[0]=1;   row.c = 1;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
                    v[0]=1;   row.c = 2;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
                } else if (i==0) {
                    v[0]=1;   row.c = 0;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
                    v[0]=1;   row.c = 1;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
                    row.c = 2;
                    v[0] = -1;              col[0].i = i;   col[0].j = j;   col[0].c = 0;
                    v[1] = 1;               col[1].i = i+1; col[1].j = j;   col[1].c = 0;
                    v[2] = -hxdhy;          col[2].i = i;   col[2].j = j;   col[2].c = 1;
                    v[3] = hxdhy;           col[3].i = i;   col[3].j = j+1; col[3].c = 1;
                    v[4] = hx;              col[4].i = i;   col[4].j = j;   col[4].c = 2;
                    ierr = MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr);
                } else {
                    row.c = 0;
                    v[0] = 1;                       col[0].i = i-1; col[0].j = j;   col[0].c = 0;
                    v[1] = 1;                       col[1].i = i+1; col[1].j = j;   col[1].c = 0;
                    v[2] = (-2-2*hxdhy*hxdhy);      col[2].i = i;   col[2].j = j;   col[2].c = 0;
                    v[3] = hxdhy*hxdhy;             col[3].i = i;   col[3].j = j+1; col[3].c = 0;
                    v[4] = hx;                      col[4].i = i-1; col[4].j = j;   col[4].c = 2;
                    v[5] = -hx;                     col[5].i = i;   col[5].j = j;   col[5].c = 2;
                    ierr = MatSetValuesStencil(jac,1,&row,6,col,v,INSERT_VALUES);CHKERRQ(ierr);
                    v[0] = 1;   row.c = 1;  ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
                    row.c = 2;
                    v[0] = -1;              col[0].i = i;   col[0].j = j;   col[0].c = 0;
                    v[1] = 1;               col[1].i = i+1; col[1].j = j;   col[1].c = 0;
                    v[2] = -hxdhy;          col[2].i = i;   col[2].j = j;   col[2].c = 1;
                    v[3] = hxdhy;           col[3].i = i;   col[3].j = j+1; col[3].c = 1;
                    ierr = MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);CHKERRQ(ierr);
                }
            } else {
                // interior of the domain
                row.c = 0;
                v[0] = 1;                   col[0].i = i-1; col[0].j = j;   col[0].c = 0;
                v[1] = 1;                   col[1].i = i+1; col[1].j = j;   col[1].c = 0;
                v[2] = (-2-2*hxdhy*hxdhy);  col[2].i = i;   col[2].j = j;   col[2].c = 0;
                v[3] = hxdhy*hxdhy;         col[3].i = i;   col[3].j = j-1; col[3].c = 0;
                v[4] = hxdhy*hxdhy;         col[4].i = i;   col[4].j = j+1; col[4].c = 0;
                v[5] = hx;                  col[5].i = i-1; col[5].j = j;   col[5].c = 2;
                v[6] = -hx;                 col[6].i = i;   col[6].j = j;   col[6].c = 2;
                ierr = MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);CHKERRQ(ierr);
                row.c = 1;
                v[0] = hydhx*hydhx;         col[0].i = i-1; col[0].j = j;   col[0].c = 1;
                v[1] = hydhx*hydhx;         col[1].i = i+1; col[1].j = j;   col[1].c = 1;
                v[2] = (-2*hydhx*hydhx-2);  col[2].i = i;   col[2].j = j;   col[2].c = 1;
                v[3] = 1;                   col[3].i = i;   col[3].j = j-1; col[3].c = 1;
                v[4] = 1;                   col[4].i = i;   col[4].j = j+1; col[4].c = 1;
                v[5] = hy;                  col[5].i = i;   col[5].j = j-1; col[5].c = 2;
                v[6] = -hy;                 col[6].i = i;   col[6].j = j;   col[6].c = 2;
                ierr = MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);CHKERRQ(ierr);
                row.c = 2;
                v[0] = -1;              col[0].i = i;   col[0].j = j;   col[0].c = 0;
                v[1] = 1;               col[1].i = i+1; col[1].j = j;   col[1].c = 0;
                v[2] = -hxdhy;          col[2].i = i;   col[2].j = j;   col[2].c = 1;
                v[3] = hxdhy;           col[3].i = i;   col[3].j = j+1; col[3].c = 1;
                ierr = MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);CHKERRQ(ierr);
            }
        }
    }
    ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    PetscFunctionReturn(0);
}