Your procedure looks correct, so I suggest you have a bug in the code. You might check this by first seeing if your program maintains a constant solution, then checking if a second order method maintains a linear solution exactly.

 @Marsha J. Berger: Prof. Berger, thanks a lot for your reply. In the attached image, you can see the PDE that I want to solve and also the flux definition that I am using. I am not very sure if this is correct for MR/AMR scheme(I am currently reading the book Finite Volume Methods for Hyperbolic Problems by Prof. Leveque to figure this out). However, when I use a uniform mesh (without refinement) with these flux definition, I tend to get correct solutions.

Also, I would like to emphasize here that, I am using this local refinement to just check if my procedure is correct or not. After this, I plan to use error estimator and grid generator to use AMR.

Last thing I would like to mention is that my error (sum of difference between new value and old value at all grid points) without refinement goes below 10-6. But when I put this fine mesh on three points, error does not go beyond some 10-3.

Once again, I appreciate your reply.

Hi Tanmay, the discretizations seems to be OK, although I do not use this notation. But the algebraic system will be the same. I will take  in assumption that you code is right and there is not bug. Fist, this error calculation is not the better way to evaluate if the  system are going to convergence. You can have a wrong solution that does not changes, and the difference  will decrease to zero.  This procedure is  suitable to  evaluate the changing in the field such next to a steady state solution, for example. You can think in your system as  [A]{x}={b} where [A] is your coefficient matrix, {x} your field variable, and {b} a vector with boundaries, sources, and initial condition (depends upon your numerical approach). With a iteractive solver, you will have  first solution {x'}, and your residue is equal to {R'} ={b} - [A]{x'}. The absolute error is the euclidian norm of {R'} : error = ||R'||, that is,  you must verify if the solution respects the  algebraic equations, or  answer "how far is my algebraic equation from zero??"  You can normalize, or even take the maximum value in {R'} and do calculate again, and again: {R''} ={b'} - [A]{x''};{R'''} ={b'} - [A]{x'''} , etc. Your boundary can be other trouble,  check if after refining the flux in the boundaries have changed. It is a Finite Volume Method, the fluxes must be conserved, if not, all procedure is no sense, because this is the basement philosophy of the method, conservation of fluxes. Think in in your interpolation in the finer mesh, some errors can be also introduced and   information be lost. 

I hope to help you in some way. Good look!

Hello Mr. Ricardo. Thanks for your reply. I have modified the picture in my comment to make it more clear. I am using explicit scheme. Without the mesh refinement, this error criterion, when converges gives satisfactory results, so in my opinion, it should not create much problem here.

This is my algorithm in detail. Let 0 be coarse grid and 1 be fine.

# Initialization on the coarse grid + Physical B.C.

# Integration on the coarse grid

    # Bilinear interpolation in space to find the initial values for fine mesh (interior plus ghost points)

    # One integration on fine mesh using the flux calculated from these initial values to get the (half) updated interior fine mesh points

    # Bilinear interpolation in space and linear interpolation in time to get the (half) updated ghost points

    # Second integration on the fine mesh to catch up with the coarse mesh time.

    # Replace coarse mesh point covered by fine mesh points by their average

    # Conservation fix-up at the coarse points nearest to interface (2 on right hand and left hand each, one on bottom and top in this case)

# Repeat

There is one point I would like to ask here. Since this method was proposed for hyperbolic PDEs. As far as I could know, PDE that I am solving is not really hyperbolic. Or is it? 

Prof. Berger, since my governing PDE is parabolic (and not hyperbolic), can I use your mesh refinement procedure? Can there be problems because of this fact?

Dear Tanmay Agrawal,

Your equation is unsteady convection-diffusion equation, which of parabolic type. Since you use explicit method (which is conditionally stable) you should satisfy a stability condition. A stability criterion of the diffusion equations in 2D case says

            dt < dx^2*dy^2/(2*lambda*(dx^2+ dy^2)

If you are at the edge of stability, be aware that if you reduce your space steps dx and dy by factor 2, you should reduce your time step by factor 4.

Also, the central differencing scheme (you use it) may introduce nonphysical oscillations into your solution if the local Peclet number is greater than two.