Re: [deal.II] Re: step-22 partial boundary conditions
Jane, I believe the formula is correct. the cubic term comes from p=z^3 being the pressure manufactured solution. so in (pI-2e) you get a z^3 term and indeed a linear term in the 2e portion. I see. I missed that there is of course also a pressure. This is how I have come to conclude that it is how i am implementing the Neumann conditions that are causing the problem, since with Dirichlet conditions the solver is good. All i can figure out is that perhaps we have calculated the vector g_N (Neumann condition) incorrectly, but that has been checked multiple times - if anyone notices something stupid we have done here then do let us know too. Well, then that is the point to debug the problem. You've established that the solution is wrong. Now you need to figure out how exactly is it wrong. Is the solution off by a constant? Is the profile wrong? If you look at the error as a cellwise constant function, is there a pattern? What happens if you multiply the prescribed boundary values by two -- is the error also multiplied by two? Does the pressure stay correct in that case, but the velocity change in some way? Is the sign of the velocity correct? If you collect enough clues, you can typically guess where the problem must lie. I will for now use vector extractors like you suggested and see if that helps at all? As you already found out, that's not helping with the bug. It just makes the code simpler. Best W. -- Wolfgang Bangerth email: bange...@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 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. For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] Re: step-22 partial boundary conditions
Following this, note that using the stress as a tensor function produced the same results/problems (same errors too as doing it with component_i) but wouldn't have thought that would have made a difference anyway... On Monday, March 5, 2018 at 6:59:43 PM UTC, Jane Lee wrote: > > Hi Wolfgang, > > I believe the formula is correct. the cubic term comes from p=z^3 being > the pressure manufactured solution. so in (pI-2e) you get a z^3 term and > indeed a linear term in the 2e portion. > > The code commpiles and the error analysis is correct with Dirichlet > conditions on the top and bottom (sides are no flux still). With the code > snippet where Neumann conditions are used at the top and bottom, the code > compiles and runs, and strangely, with the error analysis being correct for > the pressure, but not for the velocity. The velocity profile qualitatively > looks ok, but is quite off the exact solution and the error doesn't > converge with the refinement level. > > This is how I have come to conclude that it is how i am implementing the > Neumann conditions that are causing the problem, since with Dirichlet > conditions the solver is good. All i can figure out is that perhaps we have > calculated the vector g_N (Neumann condition) incorrectly, but that has > been checked multiple times - if anyone notices something stupid we have > done here then do let us know too. > > I will for now use vector extractors like you suggested and see if that > helps at all? > > Many thanks > -- 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. For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] Re: step-22 partial boundary conditions
Hi Wolfgang, I believe the formula is correct. the cubic term comes from p=z^3 being the pressure manufactured solution. so in (pI-2e) you get a z^3 term and indeed a linear term in the 2e portion. The code commpiles and the error analysis is correct with Dirichlet conditions on the top and bottom (sides are no flux still). With the code snippet where Neumann conditions are used at the top and bottom, the code compiles and runs, and strangely, with the error analysis being correct for the pressure, but not for the velocity. The velocity profile qualitatively looks ok, but is quite off the exact solution and the error doesn't converge with the refinement level. This is how I have come to conclude that it is how i am implementing the Neumann conditions that are causing the problem, since with Dirichlet conditions the solver is good. All i can figure out is that perhaps we have calculated the vector g_N (Neumann condition) incorrectly, but that has been checked multiple times - if anyone notices something stupid we have done here then do let us know too. I will for now use vector extractors like you suggested and see if that helps at all? Many thanks -- 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. For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] Re: step-22 partial boundary conditions
Jane,
1. I was wondering re dimensions because I could find a component mask
function or something similar to fe_values[velocities], eg, when using
fe_face_values which you need to apply the neumann conditions.
I see. Yes, if you apply a (vector) component mask to the FEFaceValues object,
then you will get dim-dimensional objects and only need dim-dimensional
function objects that give you the traction boundary values.
4. for the top and bottom boundaries, you are right, i do mean the stress as
mentioned in previous messages.
I see. But then the formula was still wrong, correct? (Your flow field is
quadratic in z, so if you compute the stress -- i.e., the derivative -- then
you should get something that is linear, not cubic in z.)
Either way, I am now implementing things in
the neumann BC way as in step-22. with the solution as before, i add the
following to my top boundary:
for (unsigned int face_number=0;
face_number::faces_per_cell; ++face_number)
if (cell->face(face_number)->at_boundary()
&&
(cell->face(face_number)->boundary_id() == 1))
{
fe_face_values.reinit (cell, face_number);
topstress.vector_value_list(fe_face_values.get_quadrature_points(),
topstress_values);
for (unsigned int q_point=0; q_pointwhere the top stress is g_N as in step-22 (where i calculate the normal as
plus or minus (0,1) depending on whether it's the top or bottom boundary).
This code looks correct at first glance. I think you can do better if you use
vector extractors because then the inner loop simply becomes like this:
for (unsigned int i=0; iThat is, you can do away with the component_i thing. You'll have to convert
your topstress function to be a TensorFunction<1,dim>, though, and call
topstress.value_list (..., ...)
which returns a vector of Tensor<1,dim> objects that you can use as the first
factor in the computation of local_rhs above.
I am running this off exactly using the equatoins in step-22 with the right
hand sides modified and using a direct solver instead.
The thing i have worked out is that as my program runs currently, but with
DIRICHLET conditions on top and bottom, it works. but when i try to use the
nuemann conditions as i have shown in the code snippet above does not (with
everything else kept the same). So i am wondering whether the snippet i have
shown here contains an error I have not figured out.
Nothing obvious at least. The question I will always ask when someone says "it
does not run" is "what exactly does this mean?". Does your code not compile?
Does it crash? Does it produce an error message? Does it compute a solution
that "looks wrong"? If so, how exactly does it look wrong and what do you
observe how it changes if you change things such as the magnitude of your
prescribes top stress? Etc. There is much to learn from analyzing how exactly
something is wrong.
Best
W.
--
Wolfgang Bangerth email: bange...@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
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.
For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] Re: step-22 partial boundary conditions
Hi Wolfgang, thanks so much for getting back to me
1. I was wondering re dimensions because I could find a component mask
function or something similar to fe_values[velocities], eg, when using
fe_face_values which you need to apply the neumann conditions.
2. This is my fault. I meant working in 2d with a solution that only acts
vertically, hence choosing a solution that is only dependent on y and that
is 0 in the xdirection.
3. indeed i was showing the full system right hand side component
4. for the top and bottom boundaries, you are right, i do mean the stress
as mentioned in previous messages. Either way, I am now implementing things
in the neumann BC way as in step-22. with the solution as before, i add the
following to my top boundary:
for (unsigned int face_number=0;
face_number::faces_per_cell; ++face_number)
if (cell->face(face_number)->at_boundary()
&&
(cell->face(face_number)->boundary_id() == 1))
{
fe_face_values.reinit (cell, face_number);
topstress.vector_value_list(fe_face_values.get_quadrature_points(),
topstress_values);
for (unsigned int q_point=0; q_point
>
> Jane,
>
> > Firstly, I decided not to use the normal vector for now. Since the
> > normal vector is 2D, i wasn't sure how to implement the rest so that it
> > is a 'double' since my g (neumann condition vector) has 3 components
> > when the normal vector will only have 2?
>
> I'm not sure I understand this -- in 2 space dimensions, the normal
> vector has 2 components; in 3 space dimensions, it has 3 components.
>
>
> > I went back to the beginning - step-22 tutorial and started from
> > absolute basics. Testing with 1D, i have started with exact solution
> > u=(0,-y^2)^T and p = y^3, ie purely vertical.
>
> I also don't understand this, I'm afraid: in d space dimensions, the
> velocity vector has d components. So if you're in 1d, how come your
> velocity vector has 2 components?
>
> Similarly...
>
> > Doing a manufactured
> > solution type, i get the right hand side vector f = (0, 4+3z^2, 2z)^T.
>
> the right hand side has d components. I suspect you are here showing the
> right hand side of the whole system, which should then have d+1
> components (which in 1d should be 2, not 3).
>
> Or do you want to suggest in all of this that you really have a
> two-dimensional domain but that the flow is one-dimensional?
>
>
> > with my solution, on a rectangular domain, i have no tangential stresses
> > anywhere, no flux on sides (unit normal being (/pm 1,0)), and plus and
> > minus (0, z^3 + 4z) at the top.
>
> Your domain is [0,1]^2, so I understand how there is no flux through the
> bottom because there u=[0,0]. But how do you arrive at the flux you show
> here? The velocity field is quadratic in z (y?), and so should the flux,
> but you show something that is cubic with the distance.
>
> If instead of 'flux' you really mean 'stress', then it should actually
> only be linear in z since it involves the *derivative* of the velocity.
>
> 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
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.
For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] Re: step-22 partial boundary conditions
Jane, Firstly, I decided not to use the normal vector for now. Since the normal vector is 2D, i wasn't sure how to implement the rest so that it is a 'double' since my g (neumann condition vector) has 3 components when the normal vector will only have 2? I'm not sure I understand this -- in 2 space dimensions, the normal vector has 2 components; in 3 space dimensions, it has 3 components. I went back to the beginning - step-22 tutorial and started from absolute basics. Testing with 1D, i have started with exact solution u=(0,-y^2)^T and p = y^3, ie purely vertical. I also don't understand this, I'm afraid: in d space dimensions, the velocity vector has d components. So if you're in 1d, how come your velocity vector has 2 components? Similarly... Doing a manufactured solution type, i get the right hand side vector f = (0, 4+3z^2, 2z)^T. the right hand side has d components. I suspect you are here showing the right hand side of the whole system, which should then have d+1 components (which in 1d should be 2, not 3). Or do you want to suggest in all of this that you really have a two-dimensional domain but that the flow is one-dimensional? with my solution, on a rectangular domain, i have no tangential stresses anywhere, no flux on sides (unit normal being (/pm 1,0)), and plus and minus (0, z^3 + 4z) at the top. Your domain is [0,1]^2, so I understand how there is no flux through the bottom because there u=[0,0]. But how do you arrive at the flux you show here? The velocity field is quadratic in z (y?), and so should the flux, but you show something that is cubic with the distance. If instead of 'flux' you really mean 'stress', then it should actually only be linear in z since it involves the *derivative* of the velocity. Best W. -- Wolfgang Bangerth email: bange...@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 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. For more options, visit https://groups.google.com/d/optout.
