(1) advection-diffusion 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 Runge-Kutta-Rosenbrock 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 advection-diffusion equation is linear, so, unless you
are doing something like flux-limiting, 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 advection-diffusion
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 (Crank-nicholson) 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.
snip of a large portion of text
>
>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=3DLU 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) =3D c(n)
>
>for one specie, I could use them for another specie.
yes, but why not trying a sparse lu code? there are lots
available and if the situation with your matrices is not
quite unhappy sparsity in the matrix will result in a not too strong fill in
such that you can take full advantage of this.
umfpack for example
but this trick will indeed only word if the transmission speed of the different
species is not too different.
hth
peter