Solving nonlinear heat equation with radiation

[Sylvain Mathonnière]

Dear all,

Some background :

I have been working for some time now on a project which is to solve the radiative transfer equation using the discrete ordinate method coupled with the nonlinear heat equation which include radiative terms. I am following this approach here

https://hal.archives-ouvertes.fr/hal-01273062/document

At each time step, I am solving first the radiative transfer equation and then solving the heat equation and then moving to the next time step.

I am mainly following tutorial 31 since it seems close to my problem (replacing Navier-Stokes with Radiative transfer equation and adding the nonlinear term in the heat equation), I used it to see how to couple equations in the solver and how to obtain result for previous time step.

Among the particularities of my program is the use of discontinuous galerkin method for the radiative transfer equation (RTE) and more standard galerkin method for the heat equation (HE), my constructor looks like this :

template <int dim> DG_FEM<dim>::DG_FEM ()

:fe_HE (FE_Q<dim>(1), 1), dof_handler_HE (triangulation),

fe_RTE (FE_DGQ<dim>(1), N_dir), dof_handler_RTE (triangulation)

{}

With fe_HE and fe_RTE as FEsystem. Here N_dir represents the direction discretisation used in the discrete ordinate method and for my case I have N_dir = 40. This is the best way I could come up to for this problem.

I use two dof_handler since I uses FE_Q and FE_DGQ like in tutorial 31 and so it is not too confusing for me.

I hope this part is relatively clear, if not, or if it is clear but something looks horrifiyingly wrong with the way I am doing, please do not hesitate to inquire more from me.

Getting to my question :

When time comes to solve the heat equation I use the Crank-Nicolson scheme (like in the paper, starting equation 44), but at the end, I have trouble to obtain a linear system like AT = B (A=total matrix, T=temperature, B= RHS). The paper is really vague on this part and I believe the notation is misleading from equation 56.

The way it is done is by linearising the nonlinear term T^4 using Taylor Young approximation :

T_{n+1}^4 = T_{n}^4 + 4 T_{n}^3 * (T_{n+1} - T_{n})

Keep in mind T_{n+1} are vectors and the "*" represents term to term multiplication.

Following the article I end up with the heat equation looking like :

[ M/k + A/2] T_{n+1} = [ M/k - A/2] T_{n} - theta * M * (-T_{n}^4 + 2 *T_{n}^3 * T_{n+1}) + extra_RHS

Where M is the mass matrix, A is the stiffness matrix, k the time step and everything else is constant not relevant to my question.

The annoying part of this equation is that I cannot factor all the term T_{n+1} easily on the left side since "T_{n}^3 * T_{n+1}" is actually a vector (T_{n}_i^3*T_{n+1}_i) for i €[1;total_dof]. The only way I found (and I hope it is correct) is to rewrite :

(T_{n}^3 * T_{n+1})_{1<=i<=total_dof} = Diag(T_{n}^3) * T_{n+1}.

A 2D version of what I would like to do is shown below:

I am basically creating a diagonal matrix with (T_{n}^3)_i on the diagonal; The product of this diagonal matrix with the vector T_{n+1} should give me back the vector I want.

If I can do that then I can factor the vector T_{n+1} and rewrite my equation as :

[ M/k + A/2 + 2 * Diag(T_{n}^3) ] T_{n+1} = Function_of( T_{n})

and I am happy.

My question is : How can I do this Diagonal matrix creation in a setting similar to tutorial 31 where I am looping through the cells and filling a local_matrix everytime (so I do not have access to the full T_{n} vector but just the local cell T_{n} vector) ?

Finally, generally speaking, does this approach sounds "healthy" to you ?  or should I use some Newton's approach like in tutorial 15 ? or is there something obvious I am missing, like my math is wrong or something else ?

Any help would be greatly appreciated.

Best regards,

Sylvain Mathonnière

[Sylvain Mathonnière]

It seems my image of the 2D version did not get through. It was something like  (matrix * vector = vector):

| T_{n}_1            0  |  |T_{n+1}_1|   =       |T_{n}_1 * T_{n+1}_1|

| 0           T_{n}_2   |  |T_{n+1}_2|            |T_{n}_2 * T_{n+1}_2 |

I hope it goes through nicely...

[Wolfgang Bangerth]

Sylvain,
I have not tried to follow your equations in detail, but let me point out how
this is generally done. The "right" approach is to start with the PDE, then to
use a time discretization *of the PDE* that leads to a linear occurrence of
T^{n+1} on the left hand side (possibly multiplied by factors that involve
T^n), and only then think about the spatial discretization. For the latter, as
usual you multiply by a test function and integrate by parts, then approximate
the integrals by quadrature.

From your description, it seems to me that you are trying to force a special
structure of your linear system, and what you arrive at may or may not be
correct, but you don't know why because you are skipping steps of the outline
above. For example, if you have a term of the form

T_n^3 T_{n+1}

then in general this will *not* lead to multiplying the vector of unknowns
\vec T_{n+1} by a diagonal matrix with entries equal to [T_n]_i^3 *unless you
use a particular quadrature formula*. Rather, if you use a Gauss quadrature
formula, you will end up multiplying [T_{n+1}]_i by factors that will also
include [T_n]_j for degrees of freedom j that are neighbors of i.

Best
W.


On 6/28/21 8:11 AM, Sylvain Mathonnière wrote:
>
> Dear all,
>
> *Some background :*
> **

> I have been working for some time now on a project which is to solve the

> *radiative transfer equation* using the discrete ordinate method coupled with
> the *nonlinear heat equation* which include radiative terms. I am following

> this approach here
>
> https://hal.archives-ouvertes.fr/hal-01273062/document
>
> At each time step, I am solving first the radiative transfer equation and then
> solving the heat equation and then moving to the next time step.
>
> I am mainly following tutorial 31 since it seems close to my problem
> (replacing Navier-Stokes with Radiative transfer equation and adding the
> nonlinear term in the heat equation), I used it to see how to couple equations
> in the solver and how to obtain result for previous time step.
>
> Among the particularities of my program is the use of discontinuous galerkin
> method for the radiative transfer equation (RTE) and more standard galerkin
> method for the heat equation (HE), my constructor looks like this :
>

> template <intdim> DG_FEM<dim>::DG_FEM ()

> :fe_HE (FE_Q<dim>(1), 1), dof_handler_HE (triangulation),
> fe_RTE (FE_DGQ<dim>(1), N_dir), dof_handler_RTE (triangulation)
> {}
> With fe_HE and fe_RTE as FEsystem. Here N_dir represents the direction
> discretisation used in the discrete ordinate method and for my case I have
> N_dir = 40. This is the best way I could come up to for this problem.
>
> I use two dof_handler since I uses FE_Q and FE_DGQ like in tutorial 31 and so
> it is not too confusing for me.
>
> I hope this part is relatively clear, if not, or if it is clear but something
> looks horrifiyingly wrong with the way I am doing, please do not hesitate to
> inquire more from me.
>
>

> *Getting to my question :*
> *
> *

> When time comes to solve the heat equation I use the Crank-Nicolson scheme
> (like in the paper, starting equation 44), but at the end, I have trouble to
> obtain a linear system like AT = B (A=total matrix, T=temperature, B= RHS).
> The paper is really vague on this part and I believe the notation is
> misleading from equation 56.
>
> The way it is done is by linearising the nonlinear term T^4 using Taylor Young
> approximation :
> T_{n+1}^4 = T_{n}^4 + 4 T_{n}^3 * (T_{n+1} - T_{n})
> Keep in mind T_{n+1} are vectors and the "*" represents term to term
> multiplication.
>
> Following the article I end up with the heat equation looking like :
>
> [ M/k + A/2] T_{n+1} = [ M/k - A/2] T_{n} - theta * M * (-T_{n}^4 + 2

> **T_{n}^3 * T_{n+1}*) + extra_RHS

>
> Where M is the mass matrix, A is the stiffness matrix, k the time step and
> everything else is constant not relevant to my question.
>
> The annoying part of this equation is that I cannot factor all the term
> T_{n+1} easily on the left side since "T_{n}^3 * T_{n+1}" is actually a vector
> (T_{n}_i^3*T_{n+1}_i) for i €[1;total_dof]. The only way I found (and I hope
> it is correct) is to rewrite :
>

> *(T_{n}^3 * T_{n+1})_{1<=i<=total_dof} = Diag(T_{n}^3) * T_{n+1}.*

> A 2D version of what I would like to do is shown below:

> **
> *
> *
> **
> *
> *

> I am basically creating a diagonal matrix with (T_{n}^3)_i on the diagonal;
> The product of this diagonal matrix with the vector T_{n+1} should give me
> back the vector I want.
> If I can do that then I can factor the vector T_{n+1} and rewrite my equation as :
>
> [ M/k + A/2 + 2 * Diag(T_{n}^3) ] T_{n+1} = Function_of( T_{n})
>
> and I am happy.
>

> *My question is :* How can I do this Diagonal matrix creation in a setting

> similar to tutorial 31 where I am looping through the cells and filling a
> local_matrix everytime (so I do not have access to the full T_{n} vector but
> just the local cell T_{n} vector) ?
>
> Finally, generally speaking, does this approach sounds "healthy" to you ?  or
> should I use some Newton's approach like in tutorial 15 ? or is there
> something obvious I am missing, like my math is wrong or something else ?
>
> Any help would be greatly appreciated.
>
> Best regards,
>
> Sylvain Mathonnière

> **
>
>
> --
> The deal.II project is located at http://www.dealii.org/
> <https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.dealii.org%2F&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Ce345d6aef515473cad2d08d93a3eb17f%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637604863211988910%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=j2acVzXklf8OCYAWbSysPrrH%2Ffl0UpC6MW160ru3nGM%3D&reserved=0>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fforum%2Fdealii%3Fhl%3Den&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Ce345d6aef515473cad2d08d93a3eb17f%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637604863211988910%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=T%2FlnlY2icDLxZsa12n4z3%2FHXNzTdLBRHCMlfDRY0x5g%3D&reserved=0>
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com
> <mailto:dealii+un...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/bdcb220c-251a-486d-a047-b7030597383dn%40googlegroups.com
> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fdealii%2Fbdcb220c-251a-486d-a047-b7030597383dn%2540googlegroups.com%3Futm_medium%3Demail%26utm_source%3Dfooter&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Ce345d6aef515473cad2d08d93a3eb17f%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637604863211998906%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=XU5kLsJaryhDRV0CvR118q%2BrfgU0CjorvGPLFu3xhn4%3D&reserved=0>.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Sylvain Mathonnière]

Thank you for your answer ! I was trying to follow the paper where they do spatial discretisation before time discretisation, so that is why.

Why would one typically prefer to do it the other way (time discretisation before spatial) ? I checked and it is indeed the case in the deal.II tutorials I looked at.

I realised now that my problem is not completely linked with deal.II and so my questions might seem inapproriate on such forum but I saw in the guideline that this forum is " is also open to some broader discussions on the finite element method, numerical methods, C++ and similar". So I will give my questions a shot ^^. (don't hesitate to politely send me back to another place if this is getting too far for the purpose of this forum).

So I applied the Crank Nicolson discretisation first before  multiplying by test functions and integrating. But then I am however confused how to handle a term like (T_{n}^3 * T_{n+1} , w_i)_omega
(where (A,B)_omega is the usual notation for the integral over the spatial region omega of the product A*B and w_i is my test function).

If I write T_{n} and T_{n+1} as a finite linear combinaison of w_i such that (T_{n} = sum_i([T_{n}]_i * w_i)) then T_{n}^3 is monstruous unless I used an orthogonal basis of test function (w_i,w_j) =delta_{i,j} and then I end up with T_{n}^3 = sum([T_{n}]_i^3 * w_i). If I follow my orthogonal basis for test function (is it even legal ?) I end up with:

for 1<= j <= total_dof

(T_{n}^3*T_{n+1}, w_j)_omega = sum_l( T_{n+1}_l (sum_k( T_{n}_k * w_k^3*w_l, w_j))_omega), which is not as bad. So I would have a system like A*T_{n+1} where A is a matrix with element (l,j) : sum_k( T_{n}_k * w_k^3*w_l, w_j))_omega.

If I follow this method, I end up with a linear system. However I then have 2 questions :

1- Is this even mathematically sound ? I should be able to choose my basis orthogonal as far as I can think. Can I then do that with deal.II ?

2 - Is this the approach I should be following for handling such term or should I rethink the time discretisation scheme ?

Finally, I want to thank you for your time as I notice you personally answer many questions on this forum and many others in other forums (I found you in two more forums !)

Best,

Sylvain Mathonnière

[Sylvain Mathonnière]

I realised I made some errors in my previous post. I need an orthonormal basis of test function (not just orthogonal) to do what I claimed and even then I have T_{n}^3 = sum([T_{n}]_i^3 * w_i) and then  (T_{n}^3*T_{n+1}, w_j)_omega = sum_l( T_{n+1}_l (sum_k( T_{n}_k^3 * w_k*w_l, w_j))_omega).

[Wolfgang Bangerth]

> So I applied the Crank Nicolson discretisation first before  multiplying by
> test functions and integrating. But then I am however confused how to handle a

> term like *(T_{n}^3 * T_{n+1} , w_i)_omega*
> (where *(A,B)_omega* is the usual notation for the integral over the spatial
> region omega of the product A*B and *w_i* is my test function).
>
> If I write *T_{n}* and *T_{n+1}* as a finite linear combinaison of w_i such
> that (*T_{n} = sum_i([T_{n}]_i * w_i)*) then *T_{n}^3* is monstruous unless I

You're thinking of T_n^3 as a polynomial of particular high degree, but we
don't actually care about that. We just need to evaluate it at quadrature
points, and to do that you evaluate
T_n(x_q) = \sum [T_n]_j \varphi_j(x_q)
and then take whatever value that gives you to the third power. You never need
T_n^3 as a *function*, you just need to be able to evaluate it at quadrature
points, and so questions such as whether functions are orthogonal don't
actually matter.

> 1- Is this even mathematically sound ? I should be able to choose my basis
> orthogonal as far as I can think. Can I then do that with deal.II ?

Think about it this way: You want to compute an integral
\int c(x) w_i(x) w_j(x)
where c(x) is a coefficient that in your case involves T_n^3. You won't in
general be able to find orthogonal basis functions w_k for these kind of
weights c(x). You just have to approximate the integral by quadrature.

> 2 - Is this the approach I should be following for handling such term or
> should I rethink the time discretisation scheme ?

No. Just think of the old solution as a coefficient. You probably want to take
a look at step-15, which shares many of the characteristics of what you want
to do.

Best
W.

[Sylvain Mathonnière]

Thank you for the clear answer, I understand my mistake now, it is much simpler like that.

On a side note (irrelevant now), when I was talking about orthonormal basis, I was refering to something like \int w_i(x) w_j(x)= kronecker_{i,j}, which should be doable and in which basis T_n^3 would look like T_{n}^3 = sum([T_{n}]_i^3 * w_i).

I was also considering to mimic tutorial 15 with the Newton method if this strategy (linearising T^4) was not performing well, but with your reply now I will stick to my method.

By the way, Isn't there a mistake in tutorial 15 when it says F'(u_n, δu_n) = -F(u_n) (just after "...use a damping parameter αn to get better global convergence behavior: "). Shouldn't we read F'(u_n, δu_n) *  δu_n = -F(u_n) otherwise it is not homogeneous ?

Best,

Sylvain

[Wolfgang Bangerth]

On 7/5/21 3:27 AM, Sylvain Mathonnière wrote:
>
> On a side note (irrelevant now), when I was talking about orthonormal basis, I
> was refering to something like \int w_i(x) w_j(x)= kronecker_{i,j}, which

> should be doable and in which basis T_n^3 would look like *T_{n}^3 =
> sum([T_{n}]_i^3 * w_i).*

No, that's not right. You don't have w_i w_k = w_i. The kronecker property
only holds if you integrate over things.

> By the way, Isn't there a mistake in tutorial 15 when it says F'(u_n, δu_n) =

> -F(u_n) (just after "...use a damping parameter αnto get better global

> convergence behavior: "). Shouldn't we read F'(u_n, δu_n) * δu_n = -F(u_n)
> otherwise it is not homogeneous ?

That's just a question of notation. When we say
F'(un, delta un)
then this quantity is linear in delta u_n. Other people might have written this as
J(un) delta u_n