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
<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.
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
<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/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
<mailto:dealii+unsubscr...@googlegroups.com>.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/cf020345-b304-45d2-a228-4081da2d4effo%40googlegroups.com
<https://groups.google.com/d/msgid/dealii/cf020345-b304-45d2-a228-4081da2d4effo%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/07648a56-570c-dc31-ccab-3c42596b8d45%40gmail.com.
