Diffusion Solver (Overview)
The Diffusion-Solver simulates the diffusive Influence on the System. Normally, Diffusion is time-dependent, but in the steady-state-solver, it is assumed, that the diffusion is an order faster than the cell-reactions and, that a steady-state is reached in each time-step, where the uptake-reaction of the cell and the local diffusion are in equilibrium. In that case, the Diffusivity mainly influences the uptake-radius of the cell, since a cell technically has access to a larger amount of substrate, if it diffuses fast.
The Simulation is split in time-steps. During one time-step, the medium diffuses and the uptake-rate of the cells is computed. After one time-step, the cells consume medium, grow and generate pressure (by simply taking up more space now). So, a small time-step can cause errors, because the equilibrium might not already be reached by the end of the time-step, while a large time-step causes an error, because the cells are not growing often enough.
The Diffusion-Solver works with a multigrid-algorithm, which means, that the domain is stored into multiple grids of different resolution. If you have a 9x9 domain, your solver will generate a 5x5 and a 3x3-grid from it. This is necessary, because solving the problem on the fine grid would take a lot of iterations. Solving on the coarse grid however, takes less iterations, but produces a more rough solution. Combining these two options is a good way of solving the system in an adequate time with adequate precision.
But how precise can we actually be? That is a good question, because you actually never know what the real solution should look like. All you can do is estimate the reachable precision and then iterate until you are satisfyingly close. To do this, we use a Truncation-Error and the Residual. Those two are computed pretty similar, but the Truncation-Error is our breakup-criterion, and therefore only gets computed once for a loop. The Residual however is computed after each loop and compared to the Truncation-Error. Is it low enough? If yes, success. If not, give it another run.
The Solver is based on looped processes. At the start, it creates the coarsest grid it can, using the values from the last time-step, respectively the initial conditions. How coarse this grid is depends on your domain. If you have 9x9, the coarsest will be 3x3. If you have 17x9, the coarsest will be 5x3. Every created dimension has (2^x)+1 grid-points (because is has 2^x Elements and Points=Elements+1). How coarse the coarsest grid will be is calculated by making the logarithm to base 2 of the smallest dimension of your domain. The result is the amount of grids that can be produced (called orders in the program). The maxOrder is your finest grid (your basic domain), while order 0 will be your coarsest.
The most important part is: you have your coarse grid now, and it's filled with the old values (respectively those of the bulk-layer, because those never change). This coarse grid gets relaxed (matrix-relaxation-method), which basically means, its getting iterated. After a fixed number of iterations, the grid gets interpolated to the next finer order. This new grid then gets relaxed again and some error-estimation is done in the background (for example setting the TruncationError). Now, the Solver overwrites the coarse grid, by filling in the last values again and doing the same iterations. This is important, because your previous grid (that is now stored in order 1), has gone through additional iteration and comparing these two solutions can give you the means to estimate the error (the Residual). So now, you have an iterated fine grid of order 1 and have calculated an error of the solver. If this is sufficiently low, the solver will break this cycle and go to the next solving-step. If not, it gives it another go.
Most of the time, your grid will have more than order 0 and 1, so the produced grid is not fine enough yet. That's why the Solver will now make a second loop with order 0,1 and 2. It does basically the same again, just on 3 Levels now, but this time, it is important to compute a coarse grid out of a fine one (for error-computation). This process is called restriction. It's quite similar to interpolation, but this time, you actually want to make fewer points out of many. Using an interpolation would always give you the value of the middle-node, because it's weighted much higher than the surrounding values. Therefore, there are specific weighting-algorithms, that are used to “blur” your grid and produce a coarser version, that nevertheless incorporates every value. This loop will also set a new Truncation-Error and compare the recent error to it at the end. If the precision is sufficient, the loop is done and if this fine order is actually your basic domain, you are done. If not, the loop now makes the same process with the additional order 3 and so on.
The Solver is finished, your error is ignorably small and the concentration-grid seems fine. The last thing to do now is copying this solution in your real grid. The Solver operates on computational grids, because he constantly overwrites these and applies computations to them. The real grid however stores the actually computed values and gets overwritten only, when the solver is finished. This is also the grid your values are taken from, when referring to the former time-step.
The Solver basically works with 6 internal grids. The first and most obvious one is the conc.-Grid, that stores our actual values.
The Next two are reac and diffreac. Reac is the Raction-Grid, that gets filled by updateReacRateAndDiffRate and holds the values necessary to compute the Uptake of the cells during the calculation (used in computeLop). Diffreac is the derivative of that grid and necessary, because we also compute a differential L-Operator, and therefore need this value too.
Itau and itemp are two very similar grids. They get filled by simply computing a residual during downward 1 and 2. Itau is then used to set the Truncation-error, while the norm of itemp (for temporal) is the actual residual, that is used for the break-criterion.
During one vCycle, and after the truncError has been set, itemp gets substracted from itau. This itau is then written in rhs, the last grid we use. Rhs is used during computeLop, because the change of Concentration is computed by: (lop-rhs)/dlop. While the other grids are only set to 0 in ResetMultiGrids, rhs also gets set back during initLoop, which gets called after every seperate vCycle and while the Other Grids constantly overwrite themselves, rhs is preserved during one vCycle, making it the sum of all “itau - itemp” in this vCycle.
We assume that (conc(t+1)-conc(t))/delta(t) approximates dc/dt. This is equal to the System-Constant Ddt/h² (Diffusivitytime/space²) times the gradient of the System. The best approximation for the gradient is, that grad(t+1/2) is a good estimation for the variable gradient in time. For small timesteps, grad(t+1/2) is close to grad(t), which gives us the forward-Euler formula. This might be inaccurate, but it is a good estimation, that does not require any knowledge of t+1. We use this as a Predictor. Next, we can use a Corrector to earn a more accurate Solution. Since we now have values for t+1, we can make the assumption, that grad(t+1/2) = 1/2*(grad(t)+grad(t+1)) and finally get the Crank Nicholson Formula: (conc(t+1)-conc(t))/delta(t)=Diffusivitytimestep/(space)²1/2*(gradient(t+1)+gradient(t)).
After this short introduction, you might want to have a look at the code itself, because I didn't go into detail about all the operations the Solver does. Below, you will find explanations of the most important methods inside the code.
following are the methods that are invoked by stepSolveDiffusion and make up the solver
this method sets every computational Grid to 0 (Grids, that are not used over the whole simulation, but just refilled and computed for interation-steps)
solves the Grid on the coarsest level by resetting the values to sBulk(out of computation) or the former timestep-value (always saved at the end of Solving) and then relaxing (Matrix-relaxation-method) it.
sets the values to those of the last timestep for computed points (boundary layer and biomass) and to bulk for non-computed points (perfectly mixed medium above the boundary-layer).
(confusing, because there are 2 methods called relax) this method basically just updates UpdateReacAandDiffRate and then invokes the "real" relax.
sets the values of the solutes in the neccessary computational grids to their real value. (Overall, or computed during last step)
invokes ComputeUptakeRate
solves the uptake-rate according to the kinetics specified in the solver. Influenced by specific growrthrate, substrate-concentration and existing biomass.
This is actually relaxing -> "solving" the current grid.
Calculates h (total size of the System in µm / number of grid-elements) and h2i -> for computing L-Operator. Takes an element of the grid, that is in boundarylayer or biomass-domain and sets up the neighbouring elemts, (fillDiff) using h2i for computeLop (L-Operator) and computeDiffLop (derivate)
(lop & dlop are basically the diffusive computations)
computes the residual (res) out of lop, _rhs and dlop and substracts the residual from the conc-Grid. negative concentrations get set to zero.
Formula: 0,5x(diff[point]+diff[below]) x (Value[point]-Value[below]) + 0,5x(diff[point]+diff[above]) x (Value[point]-Value[above]))-Reaction[point]
Basically makes the Average of the Diffusivity between two neighbouring Points and multiplies it by the Concentration-Difference of the two points. The Reaction is substracted from that value, because that would be the uptake of the cells at that particular point
diffLop is just the derivate of that formula, derived by the Concentration
interpolates from solved coarsest grid to one grid finer and sets _rhs to 0 (so _rhs is only preserved inside one V-Cycle)
[V-Cycle begins]
relaxes the fine grid, created by initLoop
retricts the fine grid back to the coarse Level. Computes the residual of itemp on fine order and restricts it to coarse order
Does Error-Calculation on fine grids (_itemp & _itau) and restricts them to one level coarser. Sets the truncation-error (only on first call)
computes the L-operators for _itau and fills the grid with them
substract _itemp (computed in downward 1) from itau
Restricts _rhs to coarse order (sould still be 0, so its just for the dimensions).
Now add the _itau-Grid to _rhs
Grid gets resettet and solved on coarsest Level another time. This overwrites the coarse grid, previously generated. (explained later)
Restricts the fine grid (from first solved coarse grid) and writes it into itemp. This itemp gets substracted from the second solved coarse Grid, which makes it something like a delta-conc of the two solved coarse grids. This delta-conc gehts interpolated and written into itau. Itau gets added to the fine solution of the conc-grid, so its basically: "conc+"interpolated"(delta-conc)"
all negative Concentrations get set to 0 (should be very close to zero anyway, so negative Concentrations are just a matter of accuracy)
relaxes fine Grid again
for next timestep
computes the residual ("rest-diffusion") of itemp and substracts the _rhs "actual values" from it. Basically making it a delta of another Diffusion-step with current values. Then, it computes the (euclidian) norm of that matrix for the residual of it, that gets compared to the tuncationError (set by downard 2). If res is sufficiently low, the loop breaks and the solver finishes.
bulk-grid-points get set back to bulk-value
Reminder: maxOrder would be the finest Grid, while order=0 is the coarsest
(The residual is computed servereal times, while the quite similarly calculated truncation-error is only computed once. The Truncation-Error gives a raw value for the reachable precision (epsilon) and once the change of the concentration (the residual of itemp) is smaller than that error, the solution is sufficiently convergent)
(the addition ...BoundaryLayer in a Method just tells, that the Boundary Layer is recognizes before computing. One should always use these methods, because otherwise the programm always computes the whole domain, but for some cases, the Definition of the bl-Grid might be worth looking at)
loop: the vCycle solves the Equation on two orders. _Itemp, _itau and _rhs are getting preserved, when the Grid is solved in the same level multiple times. When solution isn't convergent enough (res/truncerror), the Solver goes one step higher -> solves one additional finer order.
order: inner counter of the methods in the loop outer: the upper border of the loop maxorder: the absolute border for the solution
- The loop sets in after solveCoarsest.
- Initloop with one order higher -> order 2 gets set.
- Relax order 2.
- Restrict to order 1 (makes sense here)
- Order -1 (order now 1)
- Relax order 1
- Restrict to order 0 (doesn't make sense here, but the computation is so small, that the if-statement would cause more confusion, than it would benefit the code)
- Order -1 (order now 0) -> while-statement not longer met -> break
- solveCoarsest
- upward with order++, until order = outer