I am working on solving complex system of 2D PDEs governing the behaviour of plasma in a gas lamp during discharge. Recent tests have shown that because of steep gradients in temperature field and overall explosive nature of a process I won't go anywhere further without adding gas dynamics to the model. As far as I know from my mentor, there will be supersonic shocks traveling across the 5mm lamp. As far as I understand, to even think about somewhat accurate solution I first have to implement some kind of mesh refinement algorithm.

I have read Berger and Colella's LAMR for Shock Hydrodynamics, Colella's short course of block-structured AMR at Berkley, Structured AMR by Vakili and Martin, Colella, Anghel and Alexander's AMR for Multiscale Nonequilibrium physics. So I understand that there have to be performed somewhat untrivial operations during the interpolation from coarse grid onto finer ones, synchronisations between halfsteps of finer grids, averaging of finer grids onto coarse ones and some kind of flux correction to preserve conservation laws.

That brings to the question: What exactly happens during the time step in Adaptive Mesh Refinement algorithm? I would appreciate full explanation, but here are some specific questions that are likely to arise either way:

Consider parabolic problem on two-level grid, $L_0$ & $L_1$, $L_0$ being the coarse one.

1. Is it correct that to get least amount of errors, I have to first advance L0 problem by time step $\tau$ as if there was no mesh refinement at all, then interpolate resulting $L_0$ data to get boundary conditions for $L_1$ problem at $t+\frac{\tau}{2}$ and advance $L_1$ for $\frac{\tau}{2}$, then interpolate L0 data for $L_1$ problem at $t + \tau$ and advance $L_1$ fol half-step again, then average $L_1$ data onto L0 grid and correct $L_0$ points adjacent to L1 boundary for flux conservation using $L_1$ data? Do I understand correctly that those flux fixups do not propagate further into coarse grid?

1.a. During correction of $L_0$ data using $L_1$, do I recalculate all $L_1$-bordering L0 points at the same time, or am I allowed to use corrected data on $L_0$ $X$-point alongside with data averaged from $L_1$ grid during my correcting of $Y$-point adjacent to $X$? (To illustrate: I am talking about something similar to difference between synchronous and asynchronous Hopfield networks).

1.b. If it is even possible, how do I synchronize grids for elliptic task? Radiative transfer is a large part of my system, too.

2. How do I generate boundary conditions on fine-fine grid boundary?

2.a. In some of the papers I've read I noticed the phrase similar to "If possible, boundary conditions are taken from the same-level grid or finer". How can this be, if at the moment of me updating the first of $L_k$ grids I do not yet have any data from the other $L_k$ or $L_{k+1}$ grids?

2.b. If the very concept of fine-fine boundary is not possible or leads to significant errors, then do I have to merge those contacting fine grids and solve the problem on the non-rectangular grid?

2.c. If I am to solve for the non-rectangular grid, what SLAE solvers besides iterative are good for sparse matrices? I was using Thomas algorithm modified for block-tridiagonal matrices, but as far as I understand, matrix for generic non-rectangular grid won't be tridiagonal.