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).*
El martes, 29 de junio de 2021 a la(s) 13:24:31 UTC+2, Sylvain Mathonnière
escribió:
> 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
> El martes, 29 de junio de 2021 a la(s) 03:58:58 UTC+2, Wolfgang Bangerth
> escribió:
>
>>
>> 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 [email protected]
>> > <mailto:[email protected]>.
>> > 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: [email protected]
>> 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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/859824ee-f897-4ff3-949b-c934a0619d34n%40googlegroups.com.
