Dear All,
I am solving a cylinder symmetric Poisson problem in deal.ii.
My program is already working properly, but I found the solution by
trial&error, and my hope is to get the final understanding of what I did
via the user group.
My question is closely related to this stack exchange question
<https://scicomp.stackexchange.com/questions/41715/how-do-you-handle-the-singularity-in-polar-or-cylindrical-coordinates>
, wherethe math is done for an equation independent on the z coordinate,
but depending on theta.
The Poisson equation (strong form) in cylinder coordinates with cylinder
symmetry is
1r∂∂r(r∂∂rϕ(r,z))+∂2dz2ϕ(r,z)=−ρ(r,z)
In order to avoid the sinularity at r = 0, Prof. Bangerth suggested in
the SE link to multiply with the proper area element 2 p r dr dz in
order to get the weak form\frac{1}{r} \frac{∂}{∂r}(r \frac{∂}{∂r}
\phi(r,z)) + \frac{∂^2}{dz^2}\phi(r,z) = - \rho(r,z)
∫0R∫−∞∞φ(r,z){∂∂r(r∂∂rϕ(r,z))+r∂2∂z2ϕ(r,z)}drdz= -∫0R∫−∞∞φ(r,z)ρ(r,z)rdrdz
After integration by parts, I get the following equation (where I am not
so sure anymore if it is correct):
∫0R∫−∞∞r∂φ∂r∂ϕ∂r+r∂φ∂z∂ϕdzdrdz=∫0R∫−∞∞ϕρrdrdz\int_0^R
\int_{-\infinity}^\infinity r \frac{∂\varphi}{∂r} \frac{∂\phi}{∂r} +
r \frac{∂\varphi}{∂z}\frac{∂\phi}{dz} dr dz = \int_0^R
\int_{-\infinity}^\infinity \phi \rho r dr dz
My questions now are
a) Is the weak form above correct for the cylinder symmetric case?
b) How do I implement this in deal.ii?
c) When I look at the first term in the PDE's strong form, I can carry
out the first term r derivative by applying the derivative product rule,
which results in terms having a first and a second derivative in r,
respectively. Therefore, the solution should be in any case a linear
combination of first derivative and second derivative values.
I used tutorial step-63 with success as a starting point, where the
advection-diffusion equation has both a second derivative (diffusion)
term, and a first derivative (=advection) term. Further, I used an
analytic known solution of my problem to compare the results, and I get
reasonable solutions for e = 1 and b = 2. But why beta = 2 and not beta = 1?
Another way of looking at the problem: I found just by experimenting,
that the following to deal.ii implementations give the same (reasonable)
results for my problem:
//--- first implementation ----
double advection_direction[0] = *2.0*, double advection_direction[1] = 0;
cell_matrix(i, j) +=
fe_values.JxW(q_point)* (
(*radius* * fe_values.shape_grad(i, q_point) *
fe_values.shape_grad(j, q_point)) +
(fe_values.shape_value(i, q_point) *
(advection_direction * fe_values.shape_grad(j, q_point) );
cell_rhs(i) += *radius * radius* * some_value;
//---- second implementation, but identical results ---
double advection_direction[0] = *1.0*, double advection_direction[1] = 0;
cell_matrix(i, j) +=
fe_values.JxW(q_point)* (
(fe_values.shape_grad(i, q_point) * fe_values.shape_grad(j, q_point)) +
(*1.0 / radius* * fe_values.shape_value(i, q_point) *
(advection_factor * fe_values.shape_grad(j, q_point) );
cell_rhs(i) += some_value;
The difference between both versions seems to be a factor r² in the weak
form. But I don't understand why the advection term needs a different
value ( 1 versus 2) for both versions. My impression is that the factor
of 2 is a result of the integration by parts, when including a factor r²
(its derivative is 2r), but I don't see the actual connection between my
experimental results and the math above.
Any help or hint into the right direction is appreciated
Regards
Andreas
--
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/2cab4609-9570-4026-9355-483ada91fa1b%40gmail.com.
