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/
<http://www.math.colostate.edu/~bangerth/>
--
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
<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 dea...@googlegroups.com <javascript:>.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/4f313231-dbb3-445f-923c-9eaff17ab783o%40googlegroups.com
<https://groups.google.com/d/msgid/dealii/4f313231-dbb3-445f-923c-9eaff17ab783o%40googlegroups.com?utm_medium=email&utm_source=footer>.
--
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
<mailto:dealii+unsubscr...@googlegroups.com>.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/b764b4b7-02d2-4139-95d9-68c30ad4f2a9o%40googlegroups.com
<https://groups.google.com/d/msgid/dealii/b764b4b7-02d2-4139-95d9-68c30ad4f2a9o%40googlegroups.com?utm_medium=email&utm_source=footer>.
--
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/84b533d2-c55a-2e96-9e99-7ecd9899a784%40gmail.com.
