# Re: [deal.II] Geometric Conservation Law

Dear Martin,

I would like to start by thanking you for all of your help. Here is a link
to my fork: https://github.com/AlexanderCicchino/dealii
It is not fully working yet, I need your help for the "evaluate" as you
mentioned earlier. I have been trying multiple different ways but cannot
seem to get it to work properly.
First, I have a test setup in:
tests/mappings/mapping_q_generic_GCL_curvilinear.cc
where currently I am outputting everything, this will change in the future
and I will later on add a more complicated curvilinear mesh to test. But,
importantly, the current implementation fails the test.
Also, I made changes to source/fe/mapping_q_generic.cc as you suggested in
the function mentioned above. The last part, which is including the
evaluate is in:

https://github.com/AlexanderCicchino/dealii/blob/master/source/fe/mapping_q_generic.cc

line 1622. We have X_l*grad(X_m) evaluated at each quadrature point, now we
need to somehow have the gradient of that evaluated at the quadrature
points. I assumed after line 1622 that that gradient is written in
cyclically through it for the conservative curl form.

Please let me know how you suggest I should proceed/setup the evaluate call
at line 1624.
Also, I noticed that the conservative curl form from Kopriva is not well
posed for 2D. In the past, for 2D we would extend the grid by unit 1 in the
z direction to properly evaluate the metric terms since, for example we
need: d/(d \zeta) ( z* dy / (d\eta) ). Any suggestions on how to implement
the 2D version in mapping Q generic?

Thank you,
Alex

On Wednesday, June 24, 2020 at 9:40:55 AM UTC-4, Martin Kronbichler wrote:
>
> Dear Alex,
>
> Great! I would suggest to start by simply adding new code to the
> maybe_update_q_points_Jacobians_... function with the option to turn it off
> or on. Depending on how the final implementation will look like we might
> want to move that to a separate place, but I think it will be less
> repetitive if we use a single place.
>
> Best,
> Martin
> On 22.06.20 19:59, Alexander Cicchino wrote:
>
> Dear Martin,
>
> Thank you very much! I have been working on making the test case not
> depend on our in house flowsolver's functions.
> I think that implementing Eq. 36 the "conservative curl" form would be
> sufficient.
> Yes this procedure sounds perfect to me, and I agree with the dimension of
> the object described. I have been going through the source code that you
> sent to familiarize myself with the objects. Should I be adding to the
> function maybe_update_q_points_Jacobians_and_grads_tensor or should I
> create a new function for it?
>
> Thank you,
> Alex
>
> On Friday, June 19, 2020 at 5:09:14 AM UTC-4, Martin Kronbichler wrote:
>>
>> Dear Alex,
>>
>> Great! The first thing we need to know is the equation. I had a quick
>> look in the paper by Kopriva and I think we want to use either equation
>> (36) or (37), depending on whether we consider the conservative or
>> invariant curl form, respectively. In either case, it appears that we need
>> to do this in a two-step procedure. The first step is to compute X_l and
>> \nabla_\xi X_m, which in deal.II speak are the "q_points" and "Jacobians".
>> The implementation in mapping_q_generic.cc is a bit involved because we
>> have a slow algorithm (working for arbitrary quadrature rules) and a fast
>> one for tensor product quadrature rules. Let us consider the fast one
>> because I think we have most ingredients available, whereas we would need
>> to fill additional fields for the slow algorithm. The source code for those
>> parts is here:
>>
>>
>> https://github.com/dealii/dealii/blob/9e05a87db802ecd073bf7567d77f3491170d84b4/source/fe/mapping_q_generic.cc#L1463-L1592
>>
>> I skipped the part on the Hessians (second derivative of transformation)
>> because we won't need them. The important parts here are the extractions of
>> the positions in line 1554 and the extraction of the gradients
>> (contravariant transformations) in line 1575. Those two parts will be the
>> starting point for the second phase we need to do in addition: According to
>> the algorithm by Kopriva, we need to define this as the curl of the
>> discrete interpolation of X_l \nabla_\xi X_m. To get the curl, we need
>> another round through the SelectEvaluator::evaluate() call in that function
>> to get the reference-cell gradient of that object, from which we can then
>> collect the entries of the curl. To call into evaluate one more time, we
>> also need a new data.shape_info object that does the collocation evaluation
>> of derivatives. That should only be two lines that I can show you how and
>> where to add, so let us not worry about that part now. What is important to
>> understand first (in terms of index notation vs tensor notation) is the
>> dimension of the object. I believe that X_l \nabla_\xi X_m is a rank-two
>> tensor, so it has dim*dim components, and we compute the gradient that
>> gives us a dim * dim * dim tensor. Taking the curl in the derivative and
>> inner tensor dimension space, we get rid of one component and up with a dim
>> * dim tensor as expected. The final step we need to do is to divide by the
>> determinant of the Jacobian (contravariant vectors), because the inverse
>> Jacobian in deal.II does not yet pre-multiply with the determinant.
>>
>> Does that procedure sound reasonable to you? If yes, we could start
>> putting together the ingredients. It would be good to have a mesh (the
>> curvilinear case you were mentioning) where we can test those formulas.
>>
>> Best,
>> Martin
>> On 17.06.20 18:37, Alexander Cicchino wrote:
>>
>> Dear Martin,
>>
>> Thank you for your response. Yes I agree that only some local
>> computations are necessary to implement the identities.
>> Yes I would be interested in this feature and trying to implement it. Do
>> you have any suggestions on where I should start and overall practices I
>> should follow?
>>
>> Thank you,
>> Alex
>>
>> On Wednesday, June 17, 2020 at 1:19:29 AM UTC-4, Martin Kronbichler
>> wrote:
>>>
>>> Dear Alex,
>>>
>>> This has been on my list of things to implement and verify with deal.II
>>> over a range of examples for quite a while, so I'm glad you bringing the
>>> topic up. It is definitely true that our way to define Jacobians does not
>>> take those identities into account, but I believe we should add support for
>>> them. The nice thing is that only some local computations are necessary, so
>>> having the option to use it in the polynomial mapping classes would be
>>> great. If you would be interested in this feature and trying to implement
>>> things, I'd be happy to guide you to the right places in the code.
>>>
>>> Best,
>>> Martin
>>> On 17.06.20 06:02, Alexander Cicchino wrote:
>>>
>>> Thank you for responding Wolfgang Bangerth.
>>>
>>> The GCL condition comes from the discretized scheme satisfying
>>> free-stream preservation. I will demonstrate this for 2D below, (can be
>>> interpreted for spectral, DG, finite difference, finite volume etc):
>>> Consider the conservation law: \frac{\partial W}{\partial t} +
>>> \frac{\partial F}{\partial x} +\frac{\partial G}{\partial y} =0
>>> Transforming this to the reference computational space (x,y)->(\xi,
>>> \eta):
>>> J*\frac{\partial W}{\partial t} + J*\frac{ \partial \xi}{\partial x} *
>>> \frac{\partial F}{\partial \xi} + J * \frac{ \partial \eta}{\partial x}*
>>> \frac{\partial F}{\partial \eta} + J * \frac{ \partial \xi}{\partial y} *
>>> \frac{\partial G}{\partial \xi} + J*\frac{ \partial \eta}{\partial
>>> y}*\frac{\partial G}{\partial \eta}
>>> Putting this in conservative form results in:
>>> J\frac{\partial W}{\partial t} + \frac{\partial}{\partial \xi} (
>>> J*F*\frac{\partial \xi}{\partial x} +J*G*\frac{\partial \xi}{\partial y} )
>>> + \frac{\partial}{\partial \eta} ( J*F*\frac{\partial \eta}{\partial x}
>>> +J*G*\frac{\partial \eta}{\partial y} ) - F*( GCL in x) - G*(GCL in y) =0
>>>
>>> where GCL in x = \frac{\partial }{\partial \xi} ( det(J)* \frac{\partial
>>> \xi
>>>  }{\partial x}) + \frac{\partial }{\partial \eta}( det(J)*
>>> \frac{\partial
>>> \eta}{\partial x} )
>>> similarly for y.
>>>
>>> So for the conservative numerical scheme to satisfy free stream
>>> preservation, the GCL conditions must go to zero.
>>> For linear grids, there are no issues with the classical definition for
>>> the inverse of the Jacobian, but what Kopriva had shown (before him Thomas
>>> and Lombard), was that the metric Jacobian has to be calculated in either a
>>> "conservative curl form" or an "invariant curl form" since it reduces the
>>> GCL condition to the divergence of a curl, which is always discretely
>>> satisfied. In the paper by Kopriva, he shows this, an example in 3D:
>>>  Analytically
>>> J*\frac{\partial \xi}{\partial x} = \frac{\partial z}{\partial \zeta} *
>>> \frac{\partial y}{\partial \eta} - \frac{\partial z}{\partial \eta} *
>>> \frac{\partial y}{\partial \zeta}
>>>
>>> but the primer doesn't satisfy free-stream preservation while the latter
>>> ("conservative curl form") does.
>>>
>>> I will put together a unit test for a curvilinear grid.
>>>
>>> Thank you,
>>> Alex
>>>
>>> On Tuesday, June 16, 2020 at 10:24:59 PM UTC-4, Wolfgang Bangerth wrote:
>>>>
>>>>
>>>> Alexander,
>>>>
>>>> > I am wondering if anybody has also found that the inverse of the
>>>> Jacobian from
>>>> > FE Values, with MappingQGeneric does not satisfy the Geometric
>>>> Conservation
>>>> > Law (GCL), in the sense of:
>>>> >
>>>> > Kopriva, David A. "Metric identities and the discontinuous spectral
>>>> element
>>>> > method on curvilinear meshes." /Journal of Scientific Computing/ 26.3
>>>> (2006): 301.
>>>> >
>>>> > on curvilinear elements/manifolds in 3D.
>>>> > That is:
>>>> > \frac{\partial }{\partial \hat{x}_1} *det(J)* \frac{\partial
>>>> \hat{x}_1
>>>> > }{\partial x_1} + \frac{\partial }{\partial \hat{x}_2} *det(J)*
>>>> \frac{\partial
>>>> > \hat{x}_2}{\partial x} + \frac{\partial }{\partial \hat{x}_3} *
>>>> > det(J)*\frac{\partial \hat{x}_3 }{\partial x_1} != 0 (GCL says it
>>>> should =0,
>>>> > similarly for x_2 and x_3)
>>>> >
>>>> > If so or if not, also, has anybody found a remedy to have the inverse
>>>> of the
>>>> > Jacobian from FE Values with MappingQGeneric to satisfy the GCL.
>>>>
>>>> I'm not sure any of us have ever thought about it. (I haven't -- but I
>>>> really
>>>> shouldn't speak for anyone else.) Can you explain what this equality
>>>> represents? Why should it hold?
>>>>
>>>> I'm also unsure whether we've ever checked whether it holds (exactly or
>>>> approximately). Can you create a small test program that illustrates
>>>> the
>>>> behavior you are seeing?
>>>>
>>>> Best
>>>>   W.
>>>>
>>>> --
>>>> ------------------------------------------------------------------------
>>>>
>>>> Wolfgang Bangerth          email:                 bang...@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
>>> ---
>>> 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
>>> To view this discussion on the web visit
>>>
>>> .
>>>
>>> --
>> The deal.II project is located at http://www.dealii.org/
>> For mailing list/forum options, see
>> ---
>> 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
>> To view this discussion on the web visit
>>
>> .
>>
>> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> ---
> 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
> To view this discussion on the web visit
>
> .
>
>

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see