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+unsubscr...@googlegroups.com
<mailto:dealii+unsubscr...@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: bange...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/220e3da0-74af-db2f-8e9f-e7f422e1cd3c%40colostate.edu.
