Axisymmetric problems

[Andreas Kyritsakis]

Dear all,

I have been using Deal.II for quite a time, as the basis of the multi-scale multi-physics code FEMOCS we have been developing in our groups. Until now our code has been strictly 3D, but now we need to solve some problems with rotational symmetry and we want to adapt the code for 2D axisymmetric problems in order to save computational power.

I have been trying to figure out how to set up such a system, but to be honest I am completely lost. I understand that I have to use a different mapping than the default one when definig the FEValues object (probably a mapping that is associated with a CylindricalManifold, but I cannot figure out how to do it, especially how to associate a certain manifold to a mapping. 

In general, I think the webpage is really missing a tutorial or an example where a simple 2D axisymmetric problem is solved (e.g. just the laplace/poisson equation in such a domain). Such problems are very common and a lot of people would benefit from such a tutorial. 

Can anyone help in this? 

Thank you very much,

Andreas Kyritsakis

[Wolfgang Bangerth]

Andreas,

> I have been using Deal.II for quite a time, as the basis of the multi-scale
> multi-physics code FEMOCS

> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fveskem%2Ffemocs&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cd7d8d7ed6a644b803f8908d7639a4fb1%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C1%7C637087386726073424&sdata=oXStQieIrHcNMPw4sO%2FfnCF%2B8mxdRgGXSrprxUrMo0w%3D&reserved=0>

> we have been developing in our groups. Until now our code has been strictly
> 3D, but now we need to solve some problems with rotational symmetry and we
> want to adapt the code for 2D axisymmetric problems in order to save
> computational power.
>
> I have been trying to figure out how to set up such a system, but to be honest
> I am completely lost. I understand that I have to use a different mapping than
> the default one when definig the FEValues object (probably a mapping that is
> associated with a CylindricalManifold, but I cannot figure out how to do it,
> especially how to associate a certain manifold to a mapping.

No, if you use cylindrical coordinates, then you work in r-z space and your
domain is probably just a rectangle. But your PDE changes.

As for resources: Yours is a rather frequently asked question. Rather than
trying to repeat what I've been writing on many other occasions, have you
tried to search through the mailing list archives to find some starting point?

> In general, I think the webpage is really missing a tutorial or an example
> where a simple 2D axisymmetric problem is solved (e.g. just the
> laplace/poisson equation in such a domain). Such problems are very common and
> a lot of people would benefit from such a tutorial.

Yes, we should eventually write such a tutorial...

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Andreas Kyritsakis]

Wolfgang, 

Thank you very much for your prompt response.

You are right that one could just leave the mapping as a simple 2D and change the weak formulation. However I think that it is possible to solve an axisymmetric problem without touching the weak formulation. I have been working for quite some time with GetDP, and there that is the standard way. For axisymmetric problems one just needs to change the "Jacobian" which in GetDP stands for the mapping between real cells and the reference/unit cell. I think the two approaches in the end boil down in the same equations and system assembly. In one case the coordinate transform scaling factors are considered as a part of the PDE and need to be explicitly inserted in the weak form, while in the other they are considered as part of the cell mapping and would appear implicitly in the jacobian of the transfrom (the JxW factors in the Deal.II framework). 

From a programming point of view, I think that the latter approach would produce more generic and reusable code, that is why I first tried to pursue it. I thought that this could be done in Deal.II, similarly to GetDP, by using an appropriate Mapping object. I guess such a Mapping class is not implemented though. In case I am wrong and it is (maybe a MappingManifold associated with a CylindricalManifold), would it be possible to help on how to implement it? 

In general, how can I associate a MappingManifold object with a certain Manifold object? In the documentation, such an association is mentioned several times, but it is not specified how it can be implemented.

Cheers,

Andreas

[Wolfgang Bangerth]

Andreas,

> You are right that one could just leave the mapping as a simple 2D and change
> the weak formulation. However I think that it is possible to solve an
> axisymmetric problem without touching the weak formulation.

You can, in principle, but it requires some tricks. If you want to use an
axisymmetric formulation, then in essence you formulate your shape functions as
phi_i(r,phi,z)
so that they are only defined on a r/z mesh, with no dependence on phi -- but
nevertheless defined in 3d. The point is that you have to have a 2d mesh for
the cross section in 3d. It's not a convenient set up. It's easier to recall
that if you have

\int_\Omega F(r(x,y,z),z) dx dy dz

(i.e., with an integrand that formally depends on x,y,z, but in reality only
on r,z) can be expressed as

\int_a^b \int_0^{2pi} \int_0^R F(r,z) r dr dphi dz
= 2pi \int_a^b \int_0^R F(r,z) r dr dz

In other words, the only change to the bilinear form is the addition of a
factor of 2*pi*r to all integrals that might appear in the formulation -- not
a bothersome task.

> I have been
> working for quite some time with GetDP, and there that is the standard way.

I'm curious how that works there.

> For axisymmetric problems one just needs to change the "Jacobian" which in
> GetDP stands for the mapping between real cells and the reference/unit cell.

So maybe the only difference there is the factor of 2*pi*r to |det J| (which
is part of the JxW factor that appears in every deal.II integral).

> I
> think the two approaches in the end boil down in the same equations and system
> assembly. In one case the coordinate transform scaling factors are considered
> as a part of the PDE and need to be explicitly inserted in the weak form,
> while in the other they are considered as part of the cell mapping and would
> appear implicitly in the jacobian of the transfrom (the JxW factors in the
> Deal.II framework).

Yes, that seems reasonable.

> From a programming point of view, I think that the latter approach would
> produce more generic and reusable code, that is why I first tried to pursue
> it. I thought that this could be done in Deal.II, similarly to GetDP, by using
> an appropriate Mapping object. I guess such a Mapping class is not implemented
> though. In case I am wrong and it is (maybe a MappingManifold associated with
> a CylindricalManifold), would it be possible to help on how to implement it?

You're correct -- the mapping classes don't do that. I also think that they
shouldn't: They have a different purpose and at least at the moment I don't
quite see how that could be put into what the mappings are supposed to do. The
place to do it might be where the FEValues class combines quadrature weights
with the determinant of the Jacobian -- though I might still think that it's
not so cumbersome to just keep everything in user code -- it'll be at most a
dozen lines of code for most applications.

> In general, how can I associate a MappingManifold object with a certain
> Manifold object? In the documentation, such an association is mentioned
> several times, but it is not specified how it can be implemented.

Every piece of a triangulation (cell, face, edge) has a manifold associated
with it. MappingManifold will simply use whatever it finds on a given object
-- you don't need to explicitly associate anything.

Can you point out where the places you mention about the documentation are?
Maybe we can clarify the text there. (In fact: If you wanted to, that would
make for a nice pull request to get used to contributing to the library! ;-) )

[Andreas Kyritsakis]

Hi Wolfgang,

Thank you very much. Your answer was very helpful. 

> You're correct -- the mapping classes don't do that. I also think that they
> shouldn't: They have a different purpose and at least at the moment I don't
> quite see how that could be put into what the mappings are supposed to do. The

> place to do it might be where the FEValues class combines quadrature weights
> with the determinant of the Jacobian -- though I might still think that it's
> not so cumbersome to just keep everything in user code -- it'll be at most a
> dozen lines of code for most applications.

It is not cumbersome at all to implement it in the assembly function. I implemented it, and it is just one extra line of code in the assembly function. The big advantage I find in somehow including that 2πr factor in the JxW(), applies when having a big code that solves several PDEs on the same domain, using 3D, 2D or axisymmetric geometries in general. In that case, it would be advantageous not to need to change every single assembly function depending on what kind of geometry you have, just use a different FEValues and implicitly include that factor on the JxW(). This would be in line with the generic programming approach of Deal.II for dimension-independent implementations, using the <dim> templates. 

> Can you point out where the places you mention about the documentation are?
> Maybe we can clarify the text there.

I think it was my big misunderstanding. I misinterpreted some points in this documentation page, but it was completely my mistake and I think the page is well-written. I could not understand how to connect the manifold to the Mapping, but now it is clear to me. Thank you very much for clarifying.

> If you wanted to, that would
> make for a nice pull request to get used to contributing to the library! ;-) )

Thank you for the invitation. I will try to see whether and what I can contribute.

[Wolfgang Bangerth]

> It is not cumbersome at all to implement it in the assembly function. I
> implemented it, and it is just one extra line of code in the assembly
> function. The big advantage I find in somehow including that 2πr factor in the
> JxW(), applies when having a big code that solves several PDEs on the same
> domain, using 3D, 2D or axisymmetric geometries in general. In that case, it
> would be advantageous not to need to change every single assembly function
> depending on what kind of geometry you have, just use a different FEValues and
> implicitly include that factor on the JxW(). This would be in line with the
> generic programming approach of Deal.II for dimension-independent
> implementations, using the <dim> templates.

How does GetDP do that? Is it a global flag that is set once? What if one
wanted to solve a coupled problem where one PDE is axisymmetric but another
one is fully 3d?

> I think it was my big misunderstanding. I misinterpreted some points in this

> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdealii.org%2Fdeveloper%2Fdoxygen%2Fdeal.II%2FclassMappingManifold_1_1InternalData.html&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7C553194e059364721c87308d768ec7f3e%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637093237276121270&sdata=Eke2ZtsNex86tnv%2FVM2gjgYiOYLEFrsklc5GSRnd6Tc%3D&reserved=0>

> documentation page, but it was completely my mistake and I think the page is
> well-written. I could not understand how to connect the manifold to the
> Mapping, but now it is clear to me. Thank you very much for clarifying.
>
> > If you wanted to, that would
> > make for a nice pull request to get used to contributing to the library! ;-) )
>
> Thank you for the invitation. I will try to see whether and what I can contribute.

OK, thank you in advance! We're looking forward to it ;-)

[Andreas Kyritsakis]

Wolfgang,

> How does GetDP do that? Is it a global flag that is set once? 

In GetDP the user defines in the input script an object called "Jacobian", which can be different for each domain region. When setting up an integral term in the weak formulation then, you have to specify which jacobian method will be used to evaluate the integrals. In analogy to Deal.II terminology, you define a certain mapping or FEValues object (I am not exactly sure which class calculates the jacobian in DealII), which applies to different parts of the triangulation (e.g. separated by material_id). Then whenever evaluating integrals in a certain part of the triangulation, the corresponding mapping/FEValue object would invoke the corresponding appropriate JxW that corresponds to that triangulartion object.

> What if one wanted to solve a coupled problem where one PDE is axisymmetric but another 
> one is fully 3d? 

In that case the domains/meshes should be different right? One domain is 2D the other is 3D. In that case one has just to define the desired jacobian to the domain that is axisymmetric, and declare that this jacobian is to be used when writing the corresponding weak form terms.

[Wolfgang Bangerth]

Andreas,

> > How does GetDP do that? Is it a global flag that is set once?
>
> In GetDP the user defines in the input script an object called "Jacobian",
> which can be different for each domain region. When setting up an integral
> term in the weak formulation then, you have to specify which jacobian method
> will be used to evaluate the integrals. In analogy to Deal.II terminology, you
> define a certain mapping or FEValues object (I am not exactly sure which class
> calculates the jacobian in DealII), which applies to different parts of the
> triangulation (e.g. separated by material_id). Then whenever evaluating
> integrals in a certain part of the triangulation, the corresponding
> mapping/FEValue object would invoke the corresponding appropriate JxW that
> corresponds to that triangulartion object.

Interesting. That would be some kind of metadata that the triangulation stores
-- "I am a cross section of a cylindrical domain". The FEValues class would
then query this flag when computing the JxW values and multiply them by
2*pi*r. It's an interesting concept.

I'm not sure it's ever going to be particularly high on my TODO list, but it's
something I wouldn't be opposed to either if anyone would be interested in
implementing it!