
Re: Optimization of a finite volume differencing scheme for multispecies transport problem
Posted:
Aug 2, 2008 4:59 PM


On 1 août, 22:09, boulou...@gmail.com wrote: > On 1 août, 12:40, Olin Perry Norton <mylastn...@icet.msstate.invalid> > wrote: > > > > > boulou...@gmail.com wrote: > > > I am working on a 3d finite volume scheme for an advectiondiffusion > > > reaction problem involving a large number of chemical species (more > > > than 60) and a large domain (an big lake for example). Since this > > > scheme will be used on large problem, I want it to be as efficient as > > > possible. The linear operators are splitted in 2 : > > > > (1) advectiondiffusion is solved using a fully implicit finite volume > > > discretisation with a multigrid method for solving the linear system > > > of equations > > > (2) chemistry is solved using a RungeKuttaRosenbrock solver for > > > stiff ODE. > > > > The transport (1) actually have the following form > > > > foreach specie in speciesList { > > > construct_matrix(); > > > solve_linear_system(); > > > } > > > > and takes a lot of time on the computer. > > > > Assuming that diffusion coefficients are the same for all species, the > > > whole fluid (including all species) should follow the same path during > > > the transport. I wonder if it really necessairy to loop over all > > > species and compute the transport several time. It is possible to > > > compute the transport of the fluid once, and after reuse this > > > calculation to the different species ? > > > > I would really appreciate suggestion or reference on this. > > > What you suggest is indeed possible. > > > Suppose that C(n) is the vector of concentrations of a certain > > species at timestep n. This vector will have as many dimensions > > as there are cells in your grid, and there will be a similar > > vector for each species. I have used just one index, "n", to > > indicate the timestep, but we could clearly add others to > > indicate species and to identify grid location. > > > The advectiondiffusion equation is linear, so, unless you > > are doing something like fluxlimiting, this concentration > > vector at the next timestep, n+1, will be a linear function > > of C(n), i.e., > > > C(n+1) = T * C(n) > > > where T is a matrix. It is a square matrix  the number of > > entries in this matrix is the square of the number of grid points > > you have. The matrix T takes a concentration profile > > and diffuses and convects it by one timestep. > > > Clearly, since the advection and diffusion process is the same > > for all species, once you have computed T for one species, > > you can use the same T for all species. > > > Also note that, unless the diffusion coefficient or the flow > > velocity in your lake change with time, you can continue to use > > the same T for every timestep. > > > Mathematically, this is based on the fact that the advectiondiffusion > > equation is linear in the concentration variable  if you're familiar > > with Green's functions, this is basically a discretized version of > > a Green's function. > > > The drawback I see is that ithis method would require solving for > > and storing a large matrix  if M is the number of gridpoints > > in your problem, then T would be an MxM matrix. Perhaps there is > > a clever way around this, or maybe it is manageable. Most of the > > entries in T should be very close to zero. > > > Olin Perry Norton > > Thank you very much for these informations, really appreciated. > > If I understand, what you suggest is to do the following : > > 1) Suppose I use a fully implicit (backward Euler for example) or semi > implicit (Cranknicholson) time integrating scheme. I will use a > finite volume discretisation to create a linear system in the form > > A * c(n+1) = c(n) > > 2) Compute > > T = inverse(A) > > Using a finite volume discretisation, with an hybrid differencing > scheme (central/upwind depending on the Peclet number) for the > advective part, the matrix A will be diagonnally dominant and all > entries are positive. Hence, the inverse of A exists and could be > computer with an appropriate method (do you have any suggestions on > this ?) > > 3) loop over all species and solve > > C(n+1) = T * C(n) > > for each one. > > If my understanding of your suggestion is correct, this would be a > very elegant solution to my problem ! > > If you have reference on this method or it's application in CFD, I > would be interested to read them. > > Again thank you very much.
The problem with this approach, as you mentionned, is that if A is sparse then A^1 will almost surely be dense. Since I will need to solve this system more than 60 times on a grid of 500x500x50 the inverse matrix will probably takes a lot of memory.
Maybe saving the A=LU factorization for further use could be a better idea. It is certainly less computationally expensive, but I am not sure if L and U will be sparse.
If both matrices are sparce, I could save them and once I have solved
(LU) * c(n+1) = c(n)
for one specie, I could use them for another specie.

