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 -- 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/bdcb220c-251a-486d-a047-b7030597383dn%40googlegroups.com.
