Regarding implementation of implicit thermal conduction in Joule miniapp

[Raghavendra0117]

Dear all,

I am new to MFEM and its usage. Presently, I am interested in solving a simple linear thermal diffusion problem. I have gone through some of the existing examples and miniapps that are already available with MFEM and I found the "joule" miniapp which use mixed FE spaces (L2, RT) suits to my requirement.

I was going through the joule mini app coding to understand the actual discrete eqns that are implemented. Presently, I am only interested in the thermal diffusion part of the "miniapp", neglecting E, B and joule heating (W). I got stuck at the thermal flux solver (--I assume it is an implicit solver--). I tried to derive the corresponding equation implemented in the solver with the help of some comments mentioned in joule.cpp & joule_solver.hpp (at line 85) and I found there are some differences. I have attached my derivations here (please refer to the attached file). Please correct me, if my derivations/assumptions are wrong. Also, It will be of great help If someone can suggest any Literature/docs/reports which describes the governing equations and its discrete form.

Thank you in advance
Raghavendra Kollipara

Attachment:
Thermal_conduction_implicit_formulation.pdf

[mlstowell]

Hi, @Raghavendra0117,
Your derivation appears to be correct but we do something a little different in the joule miniapp. Your derivation assumes the solver will compute $F^{n+1}$ which is quite natural. However, our ODE solvers compute $\frac{\partial F}{\partial t}$ and then use this to compute $F^{n+1}=F^n+dt \frac{\partial F}{\partial t}$ as a separate step. This step is hidden from the miniapp inside the ODE integrator's Step method. This miniapp uses one of MFEM's built in ODE integrators but I believe this time integration scheme follows the more sophisticated schemes offered in SUNDIALS' ARKODE library.

Best wishes,
Mark

[Raghavendra0117]

Thanks a lot, @mlstowell.

I went through the code related to solver (BackwardEuler in ode.cpp) and this is my understanding so far (pl correct me if i am wrong).

1. The solver is initialized with oper which is an object of MagneticDiffusionEOperator.
2. odesolver->step calls implicitsolve method of MagnetiDiffusionEOperator that outputs dX_dt = [dP_dt, dE_dt, dB_dt, dF_dt, dT_dt.]
3. This dX_dt is used to calculate the x^{n+1} as you mentioned in your comment.

Assuming my understanding is correct, the "implicit solve" routine will return/update dX_dt.

Inside the implicit solve routine, entries of dX_dt is defined as

dT.MakeRef(&L2FESpace, dX_dt,true_offset[0]);
dF.MakeRef(&HDivFESpace, dX_dt,true_offset[1]);

a2->RecoverFEMSolution(X2,v2,F);

weakDiv->AddMult(F, temp_lf, -1.0);

pcg_m3->Mult(temp_lf, dT);
There are no lines related to calculation of dF. Even if one makes thermal flux=0.0 inside the time loop before "odesolver->step()", there is no change in result of thermal flux^{n+1}. Thus, i came to conclusion that F^{n+1} is calculated and same goes to R.H.S of equation related to dT_dt calculation.

I appreciate your help if you can suggest any literature/publications/docs/reports that describes the discrete form of the implemented equations that help my understanding.

Thank you
Raghavendra Kollipara

[mlstowell]

Hi, @Raghavendra0117,
You've certainly found something odd here. It looks like both E and F are being treated in a non-standard manner. It seems that both of these fields are being solved for directly while dE/dt and dF/dt are both set to zero (so that the X+dt dX/dt step won't alter E and F). My guess is that this is done so that the boundary conditions on E and F can be applied more easily. There will be other ramifications of this choice having to do with the order in which the fields are updated and perhaps other things.

I think we agree that improved documentation is needed here.

Best wishes,
Mark

[Raghavendra0117]

Thank you @mlstowell for your feedback. I would appreciate your help to understand the equations implemented in the code and clear some of my doubts.

1. The main steps involved in the calculation of the thermal flux (F) in the code (implicit solve) are
   a) weakDiv->MultTranspose (T,*v2) ==> $v2 = \nabla T $
   b) v2 will go to R.H.S of equation [A2][X2]=[B2] and a2->RecoverFEMSolution() gives the flux F ^{n+1}.

Please refer the attached file for my derivations for implicit thermal diffusion formulations.
Thermal_diffusion_implicit_formulation_derivation.pdf

According to my derivations (implicit), there is a missing negative sign in the implementation of flux equation (please correct me, if I am wrong). Presently, it looks like a zero flux boundary condition is considered in the implementation (there is no boundary integrator). Please confirm the same.

2. The flux (F) calculated at (n+1) is used in the calculation of (dT/dt) --> variable "dT" in the implicit solve. This probably makes the entire calculation implicit.

However, the solution (Temperature) obtained for a simple test problem is going negative/oscillatory, when run with smaller time step. I am attaching an image of such a result.

It can be observed from the flux equation that when 'dt' tends to zero (refer to Eq .4 in the attached pdf file) becomes F^{n+1} = F^{n} = -k $\nabla T$^n. That is flux used in the dT/dt calculation is essentially at t^n, making the calculation purely explicit. This was okay, if the calculations are non-oscillatory despite being explicit in nature when dt tends to zero. However, for most of the test cases we tried, the solution tends to oscillate. The oscillations in the solution is minimized by using L2 space for the temperature field with basis::positive. However, the solution goes to negative for several cases tried even after that.

Please suggest me, how to proceed to fix the following issues.

1. Discrepancy in the derivation of the discrete equations in the implicit F solver (-ve sign)
2. Accounting for the boundary conditions for the F, in a general way. For example a neuman BC, applied, etc
3. Avoiding the solutions go oscillatory /negative.

Thank you
Raghavendra Kollipara