[aspect-devel] questions regarding accuracy of no-heat flux boundary condition and help needed for the right entropy viscosity parameters

[Wolfgang Bangerth]

>>> What is the size of the heat flux through the left/right and how does
>>> it compare to the top/bottom heat flux? If this is non-zero but
>>> several orders of magnitude smaller, then that is as good as it can
>>> get. If not, then I don't know.
>> Hi Timo, thanks for you response. Yes, the left and right flux (O(1))
>> are comparable small to the top and bottom (O(10^4)). However, I still
>> think if we assigned the zero non-flux boundary, then it as least should
>> be in some lower order accuracy near zero, say (10^{-1})... But perhaps
>> you are right, considering the error from the matrix iteration solver.

As well as the discretization. In the finite element context, the heat 
flux is not zero in a pointwise sense, but only in the sense of some 
weighted integral. A heat flux of ~1 on the sides is relatively small 
compared to the heat flux of ~10^4 at the top/bottom boundaries. If you 
want to know whether this is a systematic error or just a discretization 
artifact, you can run on a sequence of finer and finer meshes. The 
top/bottom heat flux should converge to something that does not change 
with the mesh size, whereas the heat flux at the left/right should go to 
zero with mesh refinement.


>> Also I am more concerned about the right formula to compute the heat
>> flux in the  statistics postprocessor. Numerically, we are solving
>> transport equation
>>
>> \rho dT/dt - \grad\cdot (k \grad T) + \grad\cdot (\nu \grad T) + v\cdot
>> \grad T = f
>>
>> where \grad\cdot (\nu \grad T) comes from the artificial entropy
>> viscosity. Therefore, in the weak formulation, we assume
>>
>>   (-k+\nu) dT/dn = 0
>>
> I just realized that the correctly formula should be
>
> \rho dT/dt - \grad\cdot (k \grad T) - \grad\cdot (\nu \grad T) + v\cdot
> \grad T = f
>
> therefore we imposed
>
> -(k+\nu) dT/dn = 0
>
> and we only assume k and \nu positive numbers, so it could be the
> magnitude of (k +\nu) affecting the final flux calculation. I'll check
> the value dT/dn instead.

Conceptually, the boundary conditions we impose is indeed
   -(k+\nu) dT/dn = 0
In practice, you of course only want
   -k dT/dn = 0
but that is not what is imposed. If you want, the boundary condition you 
get is really
   -k dT/dn = nu dT/dn
and so the total heat flux across the boundary is not going to be zero 
but proportional to the artificial viscosity. However, because nu goes 
to zero as the mesh is refined, you will converge to the correct 
solution. This is all you can probably hope.


>> which is the no-heat flux condition. However, in the statistics
>> postprocessor, we only compute  (-k) dT/dn, so perhaps this is the
>> reason why we couldn't have a zero heat flux?

Yes, this seems correct.

Best
  W.

-- 
------------------------------------------------------------------------
Wolfgang Bangerth               email:            bangerth at math.tamu.edu
                                 www: http://www.math.tamu.edu/~bangerth/