# Re: [deal.II] Geometric Conservation Law

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
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 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