Re: [deal.II] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

On 4/15/19 9:51 AM, jane@jandj-ltd.com wrote:
> 
> I can get a table if it would be useful.

In the end, *you* need to convince yourself that something works or doesn't 
work. A table is a useful tool :-)


> I see what you mean in terms of convergence. I guess I was looking for the 
> accuracy pointwise on a boundary where the Dirichlet condition for the 
> pressure is imposed weakly. In my case, the value of the output on the 
> boundary was important, which is why I was looking at the value of the 
> pressure on the boundary.

I understand why it is important to have an accurate solution. But before 
accuracy comes convergence.


> I will try using Q2-Q2 elements, though I realise that this doesn't change 
> the 
> application of the BC itself, so the issues may be the same, but definitely 
> worth a try.

Correct. But it also has a higher convergence rate than the low-order RT 
element -- which is known to work, but converge really slowly. The lowest 
order RT element is really not a particularly good choice for anything at all.


> For a neumann boundary condition though (my problem has that on another 
> boundary), am I correct in thinking that project_boundary_values doesn't 
> allow 
> for component masking? the condition is only on u = f(grad p) and not on p 
> itself so I would need something similar to component masking the velocities.

I think this is true. The function in the library doesn't allow you to choose 
a component mask -- there is a discussion on this on the mailing list sometime 
in the last couple of weeks, if you want to read up on it.

You might have to write something yourself that does the projection.

Best
  WB


-- 

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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

Hi Wolfgang,

I can get a table if it would be useful. 

I see what you mean in terms of convergence. I guess I was looking for the 
accuracy pointwise on a boundary where the Dirichlet condition for the 
pressure is imposed weakly. In my case, the value of the output on the 
boundary was important, which is why I was looking at the value of the 
pressure on the boundary. 

I will try using Q2-Q2 elements, though I realise that this doesn't change 
the application of the BC itself, so the issues may be the same, but 
definitely worth a try. 

For a neumann boundary condition though (my problem has that on another 
boundary), am I correct in thinking that project_boundary_values doesn't 
allow for component masking? the condition is only on u = f(grad p) and not 
on p itself so I would need something similar to component masking the 
velocities. 

Many thanks
 

On Monday, April 15, 2019 at 3:23:28 PM UTC+1, Wolfgang Bangerth wrote:
>
>
> Jane, 
>
> > I continued to find out why I wasn't getting the correct applied 
> Dirichlet 
> > values on the boundary for a code very similar to step-20, where the 
> Dirichlet 
> > condition is applied weakly using 
> > 
> > for (unsigned int face_no=0; 
> > face_no::faces_per_cell; 
> > ++face_no) 
> > if (cell->at_boundary(face_no)) 
> > { 
> > fe_face_values.reinit 
> > <
> https://www.dealii.org/8.4.1/doxygen/deal.II/classFEFaceValues.html#af6e079ca7429d54433343d50bd334c3c>
>  
>
> > (cell, face_no); 
> > pressure_boundary_values 
> > .value_list (fe_face_values.get_quadrature_points 
> > <
> https://www.dealii.org/8.4.1/doxygen/deal.II/classFEValuesBase.html#a5f8732ebe2d3c6746f6de26a79cb1e45>(),
>  
>
> > boundary_values); 
> > for (unsigned int q=0; q > for (unsigned int i=0; i > local_rhs(i) += -(fe_face_values[velocities].value (i, q) * 
> > fe_face_values.normal_vector 
> > <
> https://www.dealii.org/8.4.1/doxygen/deal.II/classFEValuesBase.html#a130eea0fa89263d93b20521addc830c7>(q)
>  
>
> > * 
> > boundary_values[q] * 
> > fe_face_values.JxW 
> > <
> https://www.dealii.org/8.4.1/doxygen/deal.II/classFEValuesBase.html#ad097580a2f71878695096cc73b271b9d>(q));
>  
>
> > } 
> > 
> > 
> > I then looked at step-20 - I used the exact code but solved directly 
> instead, 
> > giving me the same results as in the tutorial (for the errors etc). 
> > 
> > Even in step-20, the boundary values aren't correct for most of the 
> tests, or 
> > rather, have a lot of error in itself. 
> > For example, at the point (1,1) which is on the boundary, p = -1.1 with 
> the 
> > given test problem. 
> > 
> > These are the values I obtained having run the code: 
> > at lowest degree (0), refinement level (RL) 3, p=-0.941992, which is 
> rather 
> > far off the -1.1 value for the Dirichlet condition applied. even at RL6, 
> p=1.07976 
> > I tried the next degree up (1), at RL3, p= -1.09984. eventually at RL6 
> for 
> > this degree, we get -1.1. 
> > Similarly, at degree 2, at RL3. p is still not accurate at p = -1.10004 
>
> Individual values are usually not particularly meaningful. What matters is 
> the 
> convergence towards the correct value. Can you produce a table for each 
> order 
> in which you show the value as a function of refinement level? 
>
> I will add that it's not clear to me that the finite element method 
> applied to 
> this problem guarantees that the error *at a single point* converges to 
> zero. 
> We know that the error converges to zero when measured in integral 
> quantities 
> such as the L2 norm, but that mean not imply pointwise convergence. 
>
>
> > How else can we imposed the Dirichlet condition (without having to use 
> very 
> > high refinement levels and degrees) on the boundary for this problem? 
>
> I don't know of other ways. The question I would ask first is whether what 
> you 
> have *converges*. That's all you can really hope for. The second question 
> is 
> whether the error is *small*. If the method converges, then the error 
> *will* 
> be small if only the mesh is fine enough; but "fine enough" may be smaller 
> than you can computationally afford, and in that case one can start to ask 
> about alternatives (such as using the Q2-Q1 element used in step-43 
> instead of 
> the RT element you have here), but the first step is to unambiguously 
> establish convergence. 
>
> Best 
>   WB 
>
> -- 
>  
> 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 

Re: [deal.II] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

Jane,

> I continued to find out why I wasn't getting the correct applied Dirichlet 
> values on the boundary for a code very similar to step-20, where the 
> Dirichlet 
> condition is applied weakly using
> 
> for (unsigned int face_no=0;
> face_no::faces_per_cell;
> ++face_no)
> if (cell->at_boundary(face_no))
> {
> fe_face_values.reinit 
> 
>  
> (cell, face_no);
> pressure_boundary_values
> .value_list (fe_face_values.get_quadrature_points 
> (),
> boundary_values);
> for (unsigned int q=0; q for (unsigned int i=0; i local_rhs(i) += -(fe_face_values[velocities].value (i, q) *
> fe_face_values.normal_vector 
> (q)
>  
> *
> boundary_values[q] *
> fe_face_values.JxW 
> (q));
> }
> 
> 
> I then looked at step-20 - I used the exact code but solved directly instead, 
> giving me the same results as in the tutorial (for the errors etc).
> 
> Even in step-20, the boundary values aren't correct for most of the tests, or 
> rather, have a lot of error in itself.
> For example, at the point (1,1) which is on the boundary, p = -1.1 with the 
> given test problem.
> 
> These are the values I obtained having run the code:
> at lowest degree (0), refinement level (RL) 3, p=-0.941992, which is rather 
> far off the -1.1 value for the Dirichlet condition applied. even at RL6, 
> p=1.07976
> I tried the next degree up (1), at RL3, p= -1.09984. eventually at RL6 for 
> this degree, we get -1.1.
> Similarly, at degree 2, at RL3. p is still not accurate at p = -1.10004

Individual values are usually not particularly meaningful. What matters is the 
convergence towards the correct value. Can you produce a table for each order 
in which you show the value as a function of refinement level?

I will add that it's not clear to me that the finite element method applied to 
this problem guarantees that the error *at a single point* converges to zero. 
We know that the error converges to zero when measured in integral quantities 
such as the L2 norm, but that mean not imply pointwise convergence.


> How else can we imposed the Dirichlet condition (without having to use very 
> high refinement levels and degrees) on the boundary for this problem?

I don't know of other ways. The question I would ask first is whether what you 
have *converges*. That's all you can really hope for. The second question is 
whether the error is *small*. If the method converges, then the error *will* 
be small if only the mesh is fine enough; but "fine enough" may be smaller 
than you can computationally afford, and in that case one can start to ask 
about alternatives (such as using the Q2-Q1 element used in step-43 instead of 
the RT element you have here), but the first step is to unambiguously 
establish convergence.

Best
  WB

-- 

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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

As an update to this thread (please let me know if you think i should start 
a new one):
I continued to find out why I wasn't getting the correct applied Dirichlet 
values on the boundary for a code very similar to step-20, where the 
Dirichlet condition is applied weakly using

for (unsigned int face_no=0;
face_no::faces_per_cell;
++face_no)
if (cell->at_boundary(face_no))
{
fe_face_values.reinit 

 
(cell, face_no);
pressure_boundary_values
.value_list (fe_face_values.get_quadrature_points 

(),
boundary_values);
for (unsigned int q=0; qhttps://www.dealii.org/8.4.1/doxygen/deal.II/classFEValuesBase.html#a130eea0fa89263d93b20521addc830c7>(q)
 
*
boundary_values[q] *
fe_face_values.JxW 

(q));
}


I then looked at step-20 - I used the exact code but solved directly 
instead, giving me the same results as in the tutorial (for the errors etc).

Even in step-20, the boundary values aren't correct for most of the tests, 
or rather, have a lot of error in itself. 
For example, at the point (1,1) which is on the boundary, p = -1.1 with the 
given test problem. 

These are the values I obtained having run the code: 
at lowest degree (0), refinement level (RL) 3, p=-0.941992, which is rather 
far off the -1.1 value for the Dirichlet condition applied. even at RL6, 
p=1.07976
I tried the next degree up (1), at RL3, p= -1.09984. eventually at RL6 for 
this degree, we get -1.1.
Similarly, at degree 2, at RL3. p is still not accurate at p = -1.10004

How else can we imposed the Dirichlet condition (without having to use very 
high refinement levels and degrees) on the boundary for this problem? 
Earlier in the thread, I mention that the boundary value I have is 
important as it is used in the next set of equations I solve after this 
darcy-like system. 

I'd appreciate any ideas or suggestions

-- 
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] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

On 3/18/19 4:30 PM, jane@jandj-ltd.com wrote:
> 
> To impose strongly - would you just 
> useVectorTools::compute_nonzero_tangential_flux_constraints with the 
> ZeroFunction?
> or is there a function similar to compute_no_normal_flux_constraints?

Yes, this will compute the constraints that correspond to u.t=0.

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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

Thanks Jean-Paul, I set n_components to dim and it ran. 

However, no difference whatsoever in the solution to when I was 
equivalently imposing in the weak form (where the tangential term 
disappears due to being zero), so I do continue to wonder whether the two 
are equivalent. Thank you

On Tuesday, March 19, 2019 at 11:35:45 AM UTC, Jean-Paul Pelteret wrote:
>
> Dear Jane,
>
> const Functions::ZeroFunction no_tang_bcs;
>
>
> This should be const Functions::ZeroFunction 
> no_tang_bcs(n_components);
>
> Best,
> Jean-Paul
>
>

-- 
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] Accuracy of Dirichlet condition for p in step-20

[Jean-Paul Pelteret]

Dear Jane,

> const Functions::ZeroFunction no_tang_bcs;

This should be const Functions::ZeroFunction no_tang_bcs(n_components);

Best,
Jean-Paul

-- 
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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

In addition to above, I thought I'd try VectorTools::compute_
nonzero_tangential_flux_constraints:

I did:
std::set no_tang_flux_boundaries;
no_tang_flux_boundaries.insert(1);
const Functions::ZeroFunction no_tang_bcs;
typename FunctionMap::type no_tang_map;
no_tang_map[1] = _tang_bcs;
VectorTools::compute_nonzero_tangential_flux_constraints (dof_handler, 0, 
no_tang_flux_boundaries, no_tang_map, constraints);


And this gives me a dimension miss match error - I'm guessing coming from 
the ZeroFunction, somehow, I thought I might need to use a 
fe.component_mask(velocities), but this isn't one of the arguments. 
Where is the dimension miss match coming from?? 



On Monday, March 18, 2019 at 10:30:32 PM UTC, jane...@jandj-ltd.com wrote:
>
> No that's fair enough. 
>
> I had thought the way I was doing it would be the equivalent of setting no 
> tangential stresses. I actually also did it this way as I wasn't sure how 
> you impose it strongly. 
>
> To impose strongly - would you just use 
> VectorTools::compute_nonzero_tangential_flux_constraints 
> with the ZeroFunction?
> or is there a function similar to compute_no_normal_flux_constraints?
>
> I'll try that and see if there's a difference. 
>
> Thanks
>
> On Monday, March 18, 2019 at 8:22:31 PM UTC, Wolfgang Bangerth wrote:
>>
>> On 3/18/19 12:37 PM, jane...@jandj-ltd.com wrote: 
>> > 
>> > For the u_t=0 condition, I had been imposing weak. So basically, I have 
>> > separated a Neumann boundary condition into: 
>> > n.[pI-2e] = (n.[pI-2e]n)n + (n.[pI-2e]t)t 
>> > and saying that the second term on the rhs is 0 so disappears, and the 
>> first 
>> > term is known weakly imposed into 
>> > topstress_values[q_point]*fe_face_values.normal_vector[q_point] etc. 
>> > Is this not a valid way? 
>>
>> I don't know. I don't, because I don't know how you impose the u.t=0 
>> boundary 
>> condition weakly. The way this is typically done is via the Nitsche 
>> method 
>> (which is essentially what DG approaches use for all interior faces as 
>> well). 
>> But I'd have to spend substantially more time than I have thinking 
>> through the 
>> implications for the case you pose here... 
>>
>> 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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

No that's fair enough. 

I had thought the way I was doing it would be the equivalent of setting no 
tangential stresses. I actually also did it this way as I wasn't sure how 
you impose it strongly. 

To impose strongly - would you just use 
VectorTools::compute_nonzero_tangential_flux_constraints 
with the ZeroFunction?
or is there a function similar to compute_no_normal_flux_constraints?

I'll try that and see if there's a difference. 

Thanks

On Monday, March 18, 2019 at 8:22:31 PM UTC, Wolfgang Bangerth wrote:
>
> On 3/18/19 12:37 PM, jane...@jandj-ltd.com  wrote: 
> > 
> > For the u_t=0 condition, I had been imposing weak. So basically, I have 
> > separated a Neumann boundary condition into: 
> > n.[pI-2e] = (n.[pI-2e]n)n + (n.[pI-2e]t)t 
> > and saying that the second term on the rhs is 0 so disappears, and the 
> first 
> > term is known weakly imposed into 
> > topstress_values[q_point]*fe_face_values.normal_vector[q_point] etc. 
> > Is this not a valid way? 
>
> I don't know. I don't, because I don't know how you impose the u.t=0 
> boundary 
> condition weakly. The way this is typically done is via the Nitsche method 
> (which is essentially what DG approaches use for all interior faces as 
> well). 
> But I'd have to spend substantially more time than I have thinking through 
> the 
> implications for the case you pose here... 
>
> 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] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

On 3/18/19 12:37 PM, jane@jandj-ltd.com wrote:
> 
> For the u_t=0 condition, I had been imposing weak. So basically, I have 
> separated a Neumann boundary condition into:
> n.[pI-2e] = (n.[pI-2e]n)n + (n.[pI-2e]t)t
> and saying that the second term on the rhs is 0 so disappears, and the first 
> term is known weakly imposed into 
> topstress_values[q_point]*fe_face_values.normal_vector[q_point] etc.
> Is this not a valid way?

I don't know. I don't, because I don't know how you impose the u.t=0 boundary 
condition weakly. The way this is typically done is via the Nitsche method 
(which is essentially what DG approaches use for all interior faces as well). 
But I'd have to spend substantially more time than I have thinking through the 
implications for the case you pose here...

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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

Hi Wolfgang,

step-20: Yes indeed I do agree that that is what I am doing. I guess I'm 
trying now to find out what else it could be that is producing: the correct 
boundary points as in the Dirichlet condition when local_rhs=0 is done 
again (overwriting), but the wrong boundary points when it isn't. 


For the u_t=0 condition, I had been imposing weak. So basically, I have 
separated a Neumann boundary condition into:
n.[pI-2e] = (n.[pI-2e]n)n + (n.[pI-2e]t)t
and saying that the second term on the rhs is 0 so disappears, and the 
first term is known weakly imposed into 
topstress_values[q_point]*fe_face_values.normal_vector[q_point] etc. 
Is this not a valid way?
(PS apologies - earlier, i had said my boundary condition is partial but 
i've imposed it in a Neumann sense. sorry for any confusion)


On Monday, March 18, 2019 at 5:17:08 PM UTC, Wolfgang Bangerth wrote:
>
>
> Jane, 
>
> > Ok, re the first step-20 bc issue, I'll have another think, but I still 
> am not 
> > sure why then it isn't giving me the exact figure, whilst my suggestion 
> is (I 
> > understand your point here and I would have said I agreed with you, but 
> my 
> > implementation does work). I've verified this part of the code using 
> exact 
> > solutions so I really have no idea what the issue is, but I guess i'll 
> have to 
> > have another think. 
>
> In your exact solution, you may simply have had zero boundary conditions 
> on 
> parts of the boundary, so removing these contributions again made no 
> difference. 
>
> I don't know if that rings any kind of bell, but if you're in 2d and you 
> have 
> a rectangular geometry, you first visit the left and right boundaries of 
> cells, and then the bottom and top boundaries. So if you do local_rhs=0 
> before 
> every boundary face, in essence you are overwriting everything you've done 
> for 
> the left/right boundaries if you are on a cell that also has a bottom or 
> top 
> boundary. 
>
>
> > On the second, step-22 issue, i don't have u.n=0 conditions. I have 
> partial 
> > boundary conditions in the form of that stated in the tutorial (1st in 
> the 
> > fourth set of Bc explanations) of: 
> > u_t = 0 
> > n.(n.(pI-2*epsilon)) = nonzero value. 
> > So this is imposed weakly, and no strong conditions are set on that 
> boundary 
> > at all, and that nonzero value equals the topstress_value I had before 
> in the 
> > weak form. 
>
> But u_t=0 is a strong boundary condition. Are you not calling a function 
> that 
> imposes it? 
>
> 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] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

Jane,

> Ok, re the first step-20 bc issue, I'll have another think, but I still am 
> not 
> sure why then it isn't giving me the exact figure, whilst my suggestion is (I 
> understand your point here and I would have said I agreed with you, but my 
> implementation does work). I've verified this part of the code using exact 
> solutions so I really have no idea what the issue is, but I guess i'll have 
> to 
> have another think.

In your exact solution, you may simply have had zero boundary conditions on 
parts of the boundary, so removing these contributions again made no difference.

I don't know if that rings any kind of bell, but if you're in 2d and you have 
a rectangular geometry, you first visit the left and right boundaries of 
cells, and then the bottom and top boundaries. So if you do local_rhs=0 before 
every boundary face, in essence you are overwriting everything you've done for 
the left/right boundaries if you are on a cell that also has a bottom or top 
boundary.


> On the second, step-22 issue, i don't have u.n=0 conditions. I have partial 
> boundary conditions in the form of that stated in the tutorial (1st in the 
> fourth set of Bc explanations) of:
> u_t = 0
> n.(n.(pI-2*epsilon)) = nonzero value.
> So this is imposed weakly, and no strong conditions are set on that boundary 
> at all, and that nonzero value equals the topstress_value I had before in the 
> weak form.

But u_t=0 is a strong boundary condition. Are you not calling a function that 
imposes it?

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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

Hi Wolfang, thanks for your time on this. 

Ok, re the first step-20 bc issue, I'll have another think, but I still am 
not sure why then it isn't giving me the exact figure, whilst my suggestion 
is (I understand your point here and I would have said I agreed with you, 
but my implementation does work). I've verified this part of the code using 
exact solutions so I really have no idea what the issue is, but I guess 
i'll have to have another think. 


On the second, step-22 issue, i don't have u.n=0 conditions. I have partial 
boundary conditions in the form of that stated in the tutorial (1st in the 
fourth set of Bc explanations) of:
u_t = 0
n.(n.(pI-2*epsilon)) = nonzero value. 
So this is imposed weakly, and no strong conditions are set on that 
boundary at all, and that nonzero value equals the topstress_value I had 
before in the weak form. 




On Saturday, March 16, 2019 at 3:53:10 PM UTC, Wolfgang Bangerth wrote:
>
> On 3/16/19 7:16 AM, jane...@jandj-ltd.com  wrote: 
> > So that's what I meant by having extra contributions. I thought you 
> needed a 
> > local_rhs = 0 after/within: 
> >   for (unsigned int face_n = 0; 
> >face_n < GeometryInfo::faces_per_cell; 
> >++face_n) 
> > if (cell->at_boundary(face_n)) 
> >   { 
> > fe_face_values.reinit(cell, face_n); 
> > 
> > In the same way that if you enforce them strongly, it basically ignores 
> the 
> > contributions to the cell vertices that are on the boundary. Because 
> this is a 
> > Dirichlet condition, should we not reset the local_rhs contribution that 
> have 
> > accumulated just prior to the boundary condition? putting local_rhs = 0 
> again 
> > within the boundary contribution terms seems to actually do the job. 
>
> You need to distinguish between Dirichlet/Neumann boundary conditions, and 
> strong/weak boundary conditions. Weak boundary conditions are the ones 
> that go 
> into the weak formulation; for the typical Laplace equation, Neumann 
> boundary 
> conditions are weak and Dirichlet boundary conditions are strong. But for 
> the 
> mixed form you are using in step-20, it is the other way around. 
>
> As for whether you should reset local_rhs=0: Think about this from a data 
> flow 
> perspective. In the block you quote above, you are computing something for 
> each face. You are *writing* it into local_rhs. But what are you *doing* 
> with 
> what you have computed? If you set local_rhs=0 between each face you 
> visit, 
> then you are simply throwing away what you have computed on previous 
> faces, 
> before you use it for something (namely, putting it into the global 
> matrix). 
>
> So what the code is doing is to first accumulate the contributions of all 
> (boundary) faces of a cell before it puts it into the global matrix. Then 
> it 
> moves to the next cell, sets local_rhs=0, accumulates contributions of all 
> (boundary) faces of that cell, puts it into the global matrix, and so on. 
>
> If you reset local_rhs=0 between each two faces, then you are computing 
> the 
> contribution of each (boundary) face alright, but because you overwrite it 
> with zeros when you move to the next face, the only thing that ends up in 
> the 
> global matrix are only the contributions of the last (boundary) face you 
> visit 
> on each cell. 
>
> If this makes a difference, you may ask why that is the case and when it 
> actually matters. Specifically, if you put local_rhs=0 right before or 
> after 
> the call to fe_face_values.reinit(...) above, then you will be zeroing out 
> any 
> computation you may have done on the previous boundary faces of the 
> current 
> cell -- but this is only going to make a difference on cells that actually 
> have two or more boundary faces, i.e., the ones in the corners of your 
> domain. 
> I don't know why that would matter in your case, but it's an observation 
> that 
> may help you think about what could be wrong in your program. 
>
>
> > For Stokes - normal component of the normal stress: 
> > Ok, that's interesting. I seems to be getting a non-zero contribution 
> when I do: 
> > 
> > local_rhs(i) = topstress_values[q] * fe_face_values.normal_vector(q)* 
> >  
> fe_face_values[velocities].value(i,q_point)* 
> > fe_face_values_rock.JxW(q); 
> > 
> > where topstress_values is the value of the normal component of the 
> normal stress. 
>
> Yes, this may be nonzero, but when you call 
> constraints.distribute_local_to_global(), it zeros out rows and columns 
> that 
> correspond to degrees of freedom with strong boundary conditions. So, for 
> example, if topstress_values*normal_vector has a nonzero normal component, 
> then you will get something nonzero here, but you have prescribed strong 
> boundary conditions of the form u.n=0, then the nonzeroness here will be 
> zeroed out by the constraints. 
>
> Best 
>   W. 
>
>
> -- 
> 

Re: [deal.II] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

On 3/16/19 7:16 AM, jane@jandj-ltd.com wrote:
> So that's what I meant by having extra contributions. I thought you needed a 
> local_rhs = 0 after/within:
>           for (unsigned int face_n = 0;
>                face_n < GeometryInfo::faces_per_cell;
>                ++face_n)
>             if (cell->at_boundary(face_n))
>               {
>                 fe_face_values.reinit(cell, face_n);
> 
> In the same way that if you enforce them strongly, it basically ignores the 
> contributions to the cell vertices that are on the boundary. Because this is 
> a 
> Dirichlet condition, should we not reset the local_rhs contribution that have 
> accumulated just prior to the boundary condition? putting local_rhs = 0 again 
> within the boundary contribution terms seems to actually do the job.

You need to distinguish between Dirichlet/Neumann boundary conditions, and 
strong/weak boundary conditions. Weak boundary conditions are the ones that go 
into the weak formulation; for the typical Laplace equation, Neumann boundary 
conditions are weak and Dirichlet boundary conditions are strong. But for the 
mixed form you are using in step-20, it is the other way around.

As for whether you should reset local_rhs=0: Think about this from a data flow 
perspective. In the block you quote above, you are computing something for 
each face. You are *writing* it into local_rhs. But what are you *doing* with 
what you have computed? If you set local_rhs=0 between each face you visit, 
then you are simply throwing away what you have computed on previous faces, 
before you use it for something (namely, putting it into the global matrix).

So what the code is doing is to first accumulate the contributions of all 
(boundary) faces of a cell before it puts it into the global matrix. Then it 
moves to the next cell, sets local_rhs=0, accumulates contributions of all 
(boundary) faces of that cell, puts it into the global matrix, and so on.

If you reset local_rhs=0 between each two faces, then you are computing the 
contribution of each (boundary) face alright, but because you overwrite it 
with zeros when you move to the next face, the only thing that ends up in the 
global matrix are only the contributions of the last (boundary) face you visit 
on each cell.

If this makes a difference, you may ask why that is the case and when it 
actually matters. Specifically, if you put local_rhs=0 right before or after 
the call to fe_face_values.reinit(...) above, then you will be zeroing out any 
computation you may have done on the previous boundary faces of the current 
cell -- but this is only going to make a difference on cells that actually 
have two or more boundary faces, i.e., the ones in the corners of your domain. 
I don't know why that would matter in your case, but it's an observation that 
may help you think about what could be wrong in your program.


> For Stokes - normal component of the normal stress:
> Ok, that's interesting. I seems to be getting a non-zero contribution when I 
> do:
> 
> local_rhs(i) = topstress_values[q] * fe_face_values.normal_vector(q)*
>                              fe_face_values[velocities].value(i,q_point)*
>                             fe_face_values_rock.JxW(q);
> 
> where topstress_values is the value of the normal component of the normal 
> stress.

Yes, this may be nonzero, but when you call 
constraints.distribute_local_to_global(), it zeros out rows and columns that 
correspond to degrees of freedom with strong boundary conditions. So, for 
example, if topstress_values*normal_vector has a nonzero normal component, 
then you will get something nonzero here, but you have prescribed strong 
boundary conditions of the form u.n=0, then the nonzeroness here will be 
zeroed out by the constraints.

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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

So that's what I meant by having extra contributions. I thought you needed 
a local_rhs = 0 after/within:
 for (unsigned int face_n = 0; 
  face_n < GeometryInfo::faces_per_cell; 
  ++face_n) 
   if (cell->at_boundary(face_n)) 
 { 
   fe_face_values.reinit(cell, face_n); 

In the same way that if you enforce them strongly, it basically ignores the 
contributions to the cell vertices that are on the boundary. Because this 
is a Dirichlet condition, should we not reset the local_rhs contribution 
that have accumulated just prior to the boundary condition? putting 
local_rhs = 0 again within the boundary contribution terms seems to 
actually do the job. 


For Stokes - normal component of the normal stress:
Ok, that's interesting. I seems to be getting a non-zero contribution when 
I do:

local_rhs(i) = topstress_values[q] * fe_face_values.normal_vector(q)*
fe_face_values[velocities].value(i,q_point)*
   fe_face_values_rock.JxW(q);

where topstress_values is the value of the normal component of the normal 
stress. 

Is this incorrect??




On Saturday, March 16, 2019 at 12:35:40 AM UTC, Wolfgang Bangerth wrote:
>
>
> Jane, 
>
> > Always setting local_rhs = 0 immediately before the below implementation 
> > would take into account all cases so that would be the best: 
> > 
> > for (unsigned int q=0; q > for (unsigned int i=0; i > local_rhs(i) += -(fe_face_values[velocities].value (i, q) * 
> > fe_face_values.normal_vector 
> > <
> https://www.dealii.org/current/doxygen/deal.II/classFEValuesBase.html#a130eea0fa89263d93b20521addc830c7>(q)
>  
>
> > * 
> > boundary_values[q] * 
> > fe_face_values.JxW 
> > <
> https://www.dealii.org/current/doxygen/deal.II/classFEValuesBase.html#ad097580a2f71878695096cc73b271b9d>(q));
>  
>
>
> But it's already there, in the marked line (lkine 494 in the current 
> sources): 
>
>  for (const auto  : dof_handler.active_cell_iterators()) 
>{ 
>  fe_values.reinit(cell); 
>  local_matrix = 0; 
>  local_rhs= 0; // *** 
>
>  right_hand_side.value_list(fe_values.get_quadrature_points(), 
> rhs_values); 
>  k_inverse.value_list(fe_values.get_quadrature_points(), 
>   k_inverse_values); 
>
>  for (unsigned int q = 0; q < n_q_points; ++q) 
>for (unsigned int i = 0; i < dofs_per_cell; ++i) 
>  ...cell integration... 
>
>
>  for (unsigned int face_n = 0; 
>   face_n < GeometryInfo::faces_per_cell; 
>   ++face_n) 
>if (cell->at_boundary(face_n)) 
>  { 
>fe_face_values.reinit(cell, face_n); 
>
>pressure_boundary_values.value_list( 
>  fe_face_values.get_quadrature_points(), 
>  boundary_values); 
>
>for (unsigned int q = 0; q < n_face_q_points; ++q) 
>  for (unsigned int i = 0; i < dofs_per_cell; ++i) 
>local_rhs(i) += ... 
>
> So all that is happening is that we accumulate contributions from all 
> faces of the cell before we put it all into the global matrix. This is 
> the right thing to do. 
>
>
> > As for the Stokes case, I am making the analogy for the INHOMOGENEOUS 
> > normal component of normal stress boundary condition. \ 
> > You are right in that the no tangential stress condition requires 
> > nothing, but then i have a nonzero value of the normal component of the 
> > normal stress which needs to be in the weak form. 
> > I can't quite figure out if this is a "dirichlet' condition so I also 
> > set local_rhs = 0 before the boundary implementation, or whether it is 
> > just an extra addition for the cells on that boundary. 
>
> It is correct that the normal component of the normal stress is nonzero. 
> But it is multiplied by a test function whose normal component is zero 
> (because the normal component of the solution is prescribed). So the 
> product of these two terms is zero, and there is no additional 
> contribution to the weak form. 
>
> 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] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

Jane,

> Always setting local_rhs = 0 immediately before the below implementation 
> would take into account all cases so that would be the best:
> 
> for (unsigned int q=0; q for (unsigned int i=0; i local_rhs(i) += -(fe_face_values[velocities].value (i, q) *
> fe_face_values.normal_vector 
> (q)
>  
> *
> boundary_values[q] *
> fe_face_values.JxW 
> (q));

But it's already there, in the marked line (lkine 494 in the current 
sources):

 for (const auto  : dof_handler.active_cell_iterators())
   {
 fe_values.reinit(cell);
 local_matrix = 0;
 local_rhs= 0; // ***

 right_hand_side.value_list(fe_values.get_quadrature_points(),
rhs_values);
 k_inverse.value_list(fe_values.get_quadrature_points(),
  k_inverse_values);

 for (unsigned int q = 0; q < n_q_points; ++q)
   for (unsigned int i = 0; i < dofs_per_cell; ++i)
 ...cell integration...


 for (unsigned int face_n = 0;
  face_n < GeometryInfo::faces_per_cell;
  ++face_n)
   if (cell->at_boundary(face_n))
 {
   fe_face_values.reinit(cell, face_n);

   pressure_boundary_values.value_list(
 fe_face_values.get_quadrature_points(),
 boundary_values);

   for (unsigned int q = 0; q < n_face_q_points; ++q)
 for (unsigned int i = 0; i < dofs_per_cell; ++i)
   local_rhs(i) += ...

So all that is happening is that we accumulate contributions from all 
faces of the cell before we put it all into the global matrix. This is 
the right thing to do.


> As for the Stokes case, I am making the analogy for the INHOMOGENEOUS 
> normal component of normal stress boundary condition. \
> You are right in that the no tangential stress condition requires 
> nothing, but then i have a nonzero value of the normal component of the 
> normal stress which needs to be in the weak form.
> I can't quite figure out if this is a "dirichlet' condition so I also 
> set local_rhs = 0 before the boundary implementation, or whether it is 
> just an extra addition for the cells on that boundary.

It is correct that the normal component of the normal stress is nonzero. 
But it is multiplied by a test function whose normal component is zero 
(because the normal component of the solution is prescribed). So the 
product of these two terms is zero, and there is no additional 
contribution to the weak form.

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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

>What line of the code would that be? Do you think it would be wrong to 
just 
>*always* set local_rhs=0? 

Always setting local_rhs = 0 immediately before the below implementation 
would take into account all cases so that would be the best:

for (unsigned int q=0; qhttps://www.dealii.org/current/doxygen/deal.II/classFEValuesBase.html#a130eea0fa89263d93b20521addc830c7>(q)
 
*
boundary_values[q] *
fe_face_values.JxW 

(q));

As for the Stokes case, I am making the analogy for the INHOMOGENEOUS 
normal component of normal stress boundary condition. \
You are right in that the no tangential stress condition requires nothing, 
but then i have a nonzero value of the normal component of the normal 
stress which needs to be in the weak form. 
I can't quite figure out if this is a "dirichlet' condition so I also set 
local_rhs = 0 before the boundary implementation, or whether it is just an 
extra addition for the cells on that boundary. 

On Friday, March 15, 2019 at 5:00:01 AM UTC, Wolfgang Bangerth wrote:
>
> On 3/14/19 5:02 PM, jane...@jandj-ltd.com  wrote: 
> > 
> > I think just a comment in the bit on how to implement the dirichlet bc 
> in the 
> > weak form would be sufficient - something to say 'In the case of an 
> > inhomogeneous boundary condition, you would need to set local_rhs = 0 
> before 
> > adding the cell contributions for the boundary condition'. 
>
> What line of the code would that be? Do you think it would be wrong to 
> just 
> *always* set local_rhs=0? 
>
>
> > I'm still unsure about the step-22 condition, for the normal component 
> for the 
> > normal stress. Is this equivalent to a dirichlet condition on the 
> pressure 
> > only? I'm a little confused on this one and any thoughts would be 
> helpful. 
>
> The normal stress shows up in the weak form after integrating by parts 
> both 
> the div(2*eta*eps(u)) and the grad p terms. The term is going to be 
> something like 
>(v, (2*eta*eps(u) - pI).n)_{Gamma} 
>
> So if you prescribe 
> * no normal flux 
> * no tangential stress 
> (i.e., "free slip"), then you will have: 
> * the normal component of v is zero 
> * the tangential component of (2*eta*eps(u) - pI).n is zero 
> As a consequence, the entire boundary term disappears, as you can see if 
> you 
> write the product of test function and normal stress as 
>(v, s.n) = (v_n, (s.n).n) + (v_t, (s.n) x n) 
> (normal component plus tangential component of vectors), where the first 
> of 
> the two terms on the right disappears because v_n=0 and the second because 
> (s.n) x n = 0. 
>
> 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] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

On 3/14/19 5:02 PM, jane@jandj-ltd.com wrote:
> 
> I think just a comment in the bit on how to implement the dirichlet bc in the 
> weak form would be sufficient - something to say 'In the case of an 
> inhomogeneous boundary condition, you would need to set local_rhs = 0 before 
> adding the cell contributions for the boundary condition'.

What line of the code would that be? Do you think it would be wrong to just 
*always* set local_rhs=0?


> I'm still unsure about the step-22 condition, for the normal component for 
> the 
> normal stress. Is this equivalent to a dirichlet condition on the pressure 
> only? I'm a little confused on this one and any thoughts would be helpful.

The normal stress shows up in the weak form after integrating by parts both 
the div(2*eta*eps(u)) and the grad p terms. The term is going to be something 
like
   (v, (2*eta*eps(u) - pI).n)_{Gamma}

So if you prescribe
* no normal flux
* no tangential stress
(i.e., "free slip"), then you will have:
* the normal component of v is zero
* the tangential component of (2*eta*eps(u) - pI).n is zero
As a consequence, the entire boundary term disappears, as you can see if you 
write the product of test function and normal stress as
   (v, s.n) = (v_n, (s.n).n) + (v_t, (s.n) x n)
(normal component plus tangential component of vectors), where the first of 
the two terms on the right disappears because v_n=0 and the second because 
(s.n) x n = 0.

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] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

Hi Jane,

> I would think that it might be useful in the tutorial!!

The same applies here: I deal with too much email to be clear about what 
"it" replies to :-) If you think "it" would be useful, can you make a 
concrete suggestion?


> Please do let me know if you have any suggestions on a perhaps analogous 
> question for step-22 in my other reply to this thread. I'm making the 
> analogy between the dirichlet boundary condition implementation and the 
> normal component of the normal stress condition for the partial boundary 
> condition in step-22.

The same would apply there as well, if it isn't already done there. So 
if you think that something ought to be changed in step-22, then I'll be 
happy to get *concrete* suggestions for that as well!

Thanks
  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] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

On 3/14/19 11:42 AM, jane@jandj-ltd.com wrote:
> 
> as an addition to this post, For Dirichlet condition that are in the 
> weak form, would this have to be done?
> 
> And for example, for an inhomogeneous normal component of the normal 
> stress condition in step-22 for Stokes (1st of the partial bc discussion 
> in the tutorial), do you also have to initialise local_rhs? I have a 
> similar condition for the top boundary with the Stokes equation but 
> can't quite figure out if the pressure should equal the normal component 
> of the normal stress or not (if the stress tensor/velocity is non-zero, 
> then it shouldn't equal the inhomogeneous condition input?)

I'm not entirely sure I know what you mean with "it" or "this" -- I 
suspect you are asking whether local_rhs should be set to zero every 
time we move to a new cell.

The answer is yes. We think of these as each cell's contribution to the 
global matrix. We want to compute it separately on each cell, rather 
than carry over some state from the previous cell.

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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

Hi Wolfgang,

I would think that it might be useful in the tutorial!!

Please do let me know if you have any suggestions on a perhaps analogous 
question for step-22 in my other reply to this thread. I'm making the 
analogy between the dirichlet boundary condition implementation and the 
normal component of the normal stress condition for the partial boundary 
condition in step-22. 

Thank you!


On Thursday, March 14, 2019 at 5:12:50 PM UTC, Wolfgang Bangerth wrote:
>
> On 3/14/19 9:13 AM, jane...@jandj-ltd.com  wrote: 
> > Note that I think this was the issue. 
> > 
> > You would need to put local_rhs = 0 within the part below to apply the 
> > Dirichlet condition. Please do let me know if you think this is wrong 
> > but I am getting the correct output. Thanks 
>
> That seems reasonable and correct. 
>
> Are you suggesting that step-20 could be made clearer by adding a line 
> that in that program doesn't make a difference but would if it were used 
> in a more general context? If so, could you suggest what to put and where? 
>
> Best & glad to hear that you solved the problem! 
>   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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

Hi all,

as an addition to this post, For Dirichlet condition that are in the weak 
form, would this have to be done?

And for example, for an inhomogeneous normal component of the normal stress 
condition in step-22 for Stokes (1st of the partial bc discussion in the 
tutorial), do you also have to initialise local_rhs? I have a similar 
condition for the top boundary with the Stokes equation but can't quite 
figure out if the pressure should equal the normal component of the normal 
stress or not (if the stress tensor/velocity is non-zero, then it shouldn't 
equal the inhomogeneous condition input?)

Thanks


On Thursday, March 14, 2019 at 3:13:59 PM UTC, jane...@jandj-ltd.com wrote:
>
> Note that I think this was the issue. 
>
> You would need to put local_rhs = 0 within the part below to apply the 
> Dirichlet condition. Please do let me know if you think this is wrong but I 
> am getting the correct output. Thanks
>
> for (unsigned int face_no=0;
> face_no::faces_per_cell;
> ++face_no)
> if (cell->face(face_no)->at_boundary()
> &&
> (cell->face(face_no)->boundary_id() == 1)) // top has boundary id 1 
> On Thursday, March 14, 2019 at 2:24:55 PM UTC, jane...@jandj-ltd.com 
> wrote:
>>
>> Hi Wolfgang, Sorry for the reply - for some reason i didn't get the 
>> notification. 
>>
>> What i mean by correct is that I do have two uncoupled equations, but 
>> that it is correct in that I had already verified it using what you 
>> suggested (by using exact solutions - MMS).
>>
>> By it fails, I mean that the top pf value output is not equal to pr as it 
>> should be, which I set in the weak form as in my first message in my 
>> thread: 
>>
>> local_rhs(i) += -( fe_face_values_pf[velocities].value (i, q) *
>> fe_face_values_pf.normal_vector(q) *
>> pr_boundary_values[q] *
>> fe_face_values_pf.JxW(q));
>>
>>  this is in the same way as step-20.
>>
>> I have now done further tests by just letting pr_boundary_values[q] to 
>> 1.0 in the blurb copied here - I do not get 1.0 as my value for the 
>> solution for pf  on boundary with id 1, which I believe I should get. (i 
>> get something close but not 1.0. eg 0.96..)
>>
>> I am now wondering whether I am doubling up on the local_rhs in the 
>> assembling of cells on the boundary. To clarify, in step-20, f = 0 so that 
>> in the rhs_values, there is nothing added to the cells on the boundary 
>> other than that which is put in the boundary condition. for my problem, f 
>> is not equal to zero, so I am wondering whether I have extra contributions 
>> to the cells on the boundary and that's where the difference comes from. 
>>
>> Thank you for your help.
>>
>>
>>
>>
>>
>>
>>
>>
>> On Friday, February 22, 2019 at 10:51:11 PM UTC, Wolfgang Bangerth wrote:
>>>
>>> On 2/22/19 4:32 AM, jane...@jandj-ltd.com wrote: 
>>> > 
>>> > Wolfgang, so that's exactly what I had done with MMS, and that was 
>>> > verified, so i assumed the way i was imposing it was correct. 
>>>
>>> When you say "correct", you mean that you have two uncoupled equations? 
>>>
>>>
>>> > It then 
>>> > fails to get the correct value at the top. Initially, it's correct up 
>>> to 
>>> > about 4 decimal places, then to 2, then to the unit value. 
>>>
>>> Try to be precise. When you say "It then fails...", what has changed? 
>>> Above you say that things are correct, and now they are no longer. What 
>>> changes have you made in between? And what is "then to 4,... then to 2" 
>>> -- you are suggesting that something changes form one step to the next, 
>>> but what is it? 
>>>
>>> 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] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

On 3/14/19 9:13 AM, jane@jandj-ltd.com wrote:
> Note that I think this was the issue.
> 
> You would need to put local_rhs = 0 within the part below to apply the 
> Dirichlet condition. Please do let me know if you think this is wrong 
> but I am getting the correct output. Thanks

That seems reasonable and correct.

Are you suggesting that step-20 could be made clearer by adding a line 
that in that program doesn't make a difference but would if it were used 
in a more general context? If so, could you suggest what to put and where?

Best & glad to hear that you solved the problem!
  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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

Note that I think this was the issue. 

You would need to put local_rhs = 0 within the part below to apply the 
Dirichlet condition. Please do let me know if you think this is wrong but I 
am getting the correct output. Thanks

for (unsigned int face_no=0;
face_no::faces_per_cell;
++face_no)
if (cell->face(face_no)->at_boundary()
&&
(cell->face(face_no)->boundary_id() == 1)) // top has boundary id 1 
On Thursday, March 14, 2019 at 2:24:55 PM UTC, jane...@jandj-ltd.com wrote:
>
> Hi Wolfgang, Sorry for the reply - for some reason i didn't get the 
> notification. 
>
> What i mean by correct is that I do have two uncoupled equations, but that 
> it is correct in that I had already verified it using what you suggested 
> (by using exact solutions - MMS).
>
> By it fails, I mean that the top pf value output is not equal to pr as it 
> should be, which I set in the weak form as in my first message in my 
> thread: 
>
> local_rhs(i) += -( fe_face_values_pf[velocities].value (i, q) *
> fe_face_values_pf.normal_vector(q) *
> pr_boundary_values[q] *
> fe_face_values_pf.JxW(q));
>
>  this is in the same way as step-20.
>
> I have now done further tests by just letting pr_boundary_values[q] to 1.0 
> in the blurb copied here - I do not get 1.0 as my value for the solution 
> for pf  on boundary with id 1, which I believe I should get. (i get 
> something close but not 1.0. eg 0.96..)
>
> I am now wondering whether I am doubling up on the local_rhs in the 
> assembling of cells on the boundary. To clarify, in step-20, f = 0 so that 
> in the rhs_values, there is nothing added to the cells on the boundary 
> other than that which is put in the boundary condition. for my problem, f 
> is not equal to zero, so I am wondering whether I have extra contributions 
> to the cells on the boundary and that's where the difference comes from. 
>
> Thank you for your help.
>
>
>
>
>
>
>
>
> On Friday, February 22, 2019 at 10:51:11 PM UTC, Wolfgang Bangerth wrote:
>>
>> On 2/22/19 4:32 AM, jane...@jandj-ltd.com wrote: 
>> > 
>> > Wolfgang, so that's exactly what I had done with MMS, and that was 
>> > verified, so i assumed the way i was imposing it was correct. 
>>
>> When you say "correct", you mean that you have two uncoupled equations? 
>>
>>
>> > It then 
>> > fails to get the correct value at the top. Initially, it's correct up 
>> to 
>> > about 4 decimal places, then to 2, then to the unit value. 
>>
>> Try to be precise. When you say "It then fails...", what has changed? 
>> Above you say that things are correct, and now they are no longer. What 
>> changes have you made in between? And what is "then to 4,... then to 2" 
>> -- you are suggesting that something changes form one step to the next, 
>> but what is it? 
>>
>> 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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

Hi Wolfgang, Sorry for the reply - for some reason i didn't get the 
notification. 

What i mean by correct is that I do have two uncoupled equations, but that 
it is correct in that I had already verified it using what you suggested 
(by using exact solutions - MMS).

By it fails, I mean that the top pf value output is not equal to pr as it 
should be, which I set in the weak form as in my first message in my 
thread: 

local_rhs(i) += -( fe_face_values_pf[velocities].value (i, q) *
fe_face_values_pf.normal_vector(q) *
pr_boundary_values[q] *
fe_face_values_pf.JxW(q));

 this is in the same way as step-20.

I have now done further tests by just letting pr_boundary_values[q] to 1.0 
in the blurb copied here - I do not get 1.0 as my value for the solution 
for pf  on boundary with id 1, which I believe I should get. (i get 
something close but not 1.0. eg 0.96..)

I am now wondering whether I am doubling up on the local_rhs in the 
assembling of cells on the boundary. To clarify, in step-20, f = 0 so that 
in the rhs_values, there is nothing added to the cells on the boundary 
other than that which is put in the boundary condition. for my problem, f 
is not equal to zero, so I am wondering whether I have extra contributions 
to the cells on the boundary and that's where the difference comes from. 

Thank you for your help.








On Friday, February 22, 2019 at 10:51:11 PM UTC, Wolfgang Bangerth wrote:
>
> On 2/22/19 4:32 AM, jane...@jandj-ltd.com  wrote: 
> > 
> > Wolfgang, so that's exactly what I had done with MMS, and that was 
> > verified, so i assumed the way i was imposing it was correct. 
>
> When you say "correct", you mean that you have two uncoupled equations? 
>
>
> > It then 
> > fails to get the correct value at the top. Initially, it's correct up to 
> > about 4 decimal places, then to 2, then to the unit value. 
>
> Try to be precise. When you say "It then fails...", what has changed? 
> Above you say that things are correct, and now they are no longer. What 
> changes have you made in between? And what is "then to 4,... then to 2" 
> -- you are suggesting that something changes form one step to the next, 
> but what is it? 
>
> 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] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

On 2/22/19 4:32 AM, jane@jandj-ltd.com wrote:
> 
> Wolfgang, so that's exactly what I had done with MMS, and that was 
> verified, so i assumed the way i was imposing it was correct.

When you say "correct", you mean that you have two uncoupled equations?


> It then 
> fails to get the correct value at the top. Initially, it's correct up to 
> about 4 decimal places, then to 2, then to the unit value.

Try to be precise. When you say "It then fails...", what has changed? 
Above you say that things are correct, and now they are no longer. What 
changes have you made in between? And what is "then to 4,... then to 2" 
-- you are suggesting that something changes form one step to the next, 
but what is it?

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] Accuracy of Dirichlet condition for p in step-20

[jane . lee]

Hi Both,

Thanks for your replies. 

Daniel, the weak forms are as in the steps, but with Neumann ocnditions on 
boundaries other than the top. With the chosen spaces, I believe the 
problem is well-posed. 

Wolfgang, so that's exactly what I had done with MMS, and that was 
verified, so i assumed the way i was imposing it was correct. It then fails 
to get the correct value at the top. Initially, it's correct up to about 4 
decimal places, then to 2, then to the unit value. I'm not sure why this 
would be happening when I am setting system_rhs and system_matrix to 0 each 
time and the only top boundary condition is as written above into the weak 
form. Doing std::cout of the value, it is inputting in the correct value, 
so not sure where the rest is coming from.

I'd be grateful for any further suggestions...!!


On Thursday, February 21, 2019 at 5:55:44 PM UTC+3, Wolfgang Bangerth wrote:
>
>
> > I am trying to solve a system of equations that do this: 
> > Stokes to solve for v_r and p_r for one fluid (viscous rock), I use 
> these 
> > solutions on the RHS of a Darcy type equation solved like step-20 for 
> the 
> > pressure p_f in the fluid in the domain. Using the 3 solutions, I update 
> > another component (porosity - using a simple time dependent equation 
> like 
> > dphi/dt = rhs(pr-pf)) and put it back in the Stokes equation and so 
> on 
> > 
> > One of my conditions for the Darcy eqn is that for the top boundary, the 
> > pressure there equals pr found previously. 
> > 
> > Looking at the output, my values for p_f at the top boundary does not 
> equal 
> > the p_r just found. This is causing a massive issue as my porosity 
> equation 
> > needs pr-pf to be 0 (or at least very small) at the top boundary, and 
> this 
> > seems to be causing a real mess. I'm guessing because it's weakly 
> imposed, 
> > there might be issues. 
> > 
> > Does anyone have any suggestions on what to do for this? 
> > 
> > My base profile for the equations separately all work and have been 
> verified 
> > with MMS. With test values of the porosity, I have also verified that 
> the 
> > profiles separately give me what I expect (which then gives me pr-pf 
> monotonic 
> > at least, but with what is happening at the top, this bit is messed up). 
>
> Jane -- the problem as you state it is too difficult to debug. Here's ho 
> to 
> approach this: 
>
> * Pick a problem for which you know the exact solution. 
> * In the Stokes equation, in all of the places where you currently use the 
> Darcy solution, use the exact solution instead 
> * In the Darcy equation, in all of the places where you currently use the 
> Stokes solution, use the exact solution instead 
>
> With this scheme, your iteration needs to converge in one step. If it 
> doesn't, 
> then you know where to look. If it does, revert the steps above one by one 
> until the problem stops working. 
>
> I hope this helps! 
> 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] Accuracy of Dirichlet condition for p in step-20

[Wolfgang Bangerth]

> I am trying to solve a system of equations that do this:
> Stokes to solve for v_r and p_r for one fluid (viscous rock), I use these 
> solutions on the RHS of a Darcy type equation solved like step-20 for the 
> pressure p_f in the fluid in the domain. Using the 3 solutions, I update 
> another component (porosity - using a simple time dependent equation like 
> dphi/dt = rhs(pr-pf)) and put it back in the Stokes equation and so on
> 
> One of my conditions for the Darcy eqn is that for the top boundary, the 
> pressure there equals pr found previously.
> 
> Looking at the output, my values for p_f at the top boundary does not equal 
> the p_r just found. This is causing a massive issue as my porosity equation 
> needs pr-pf to be 0 (or at least very small) at the top boundary, and this 
> seems to be causing a real mess. I'm guessing because it's weakly imposed, 
> there might be issues.
> 
> Does anyone have any suggestions on what to do for this?
> 
> My base profile for the equations separately all work and have been verified 
> with MMS. With test values of the porosity, I have also verified that the 
> profiles separately give me what I expect (which then gives me pr-pf 
> monotonic 
> at least, but with what is happening at the top, this bit is messed up).

Jane -- the problem as you state it is too difficult to debug. Here's ho to 
approach this:

* Pick a problem for which you know the exact solution.
* In the Stokes equation, in all of the places where you currently use the 
Darcy solution, use the exact solution instead
* In the Darcy equation, in all of the places where you currently use the 
Stokes solution, use the exact solution instead

With this scheme, your iteration needs to converge in one step. If it doesn't, 
then you know where to look. If it does, revert the steps above one by one 
until the problem stops working.

I hope this helps!
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.