Re: [deal.II] Problem with "constraints.distribute_local_to_global(local_rhs, local_dof_indices, system_rhs)"

[Magdalini Ntetsika]

I see. Thank you Daniel!

Best,
Magda

On Thursday, March 5, 2020 at 1:39:26 PM UTC-8, Daniel Arndt wrote:
>
> Magdalini,
>
> The problem only appears with inhomogeneous boundary conditions, but 
> periodic boundary conditions don't set any inhomogeneity.
>
> Best,
> Daniel
>
> Am Do., 5. März 2020 um 16:13 Uhr schrieb Magdalini Ntetsika <
> ntet...@berkeley.edu >:
>
>> Hi Wolfgang,
>>
>> thank you for the prompt reply. I see, I'll look up both. But I have a 
>> quick question: what if I have periodic boundary conditions all around my 
>> domain? It seems to be working in that case, even if the boundary values 
>> are nonzero, I get to have the correct system_rhs vector. I guess will need 
>> first to watch the topic on Dirichlet Boundary conditions to understand 
>> better, but it will be very helpful if I could get some feedback on that 
>> specific observation from you.
>>
>> Best,
>> Magda
>>
>> On Thursday, March 5, 2020 at 9:41:11 AM UTC-8, Wolfgang Bangerth wrote:
>>>
>>> On 3/5/20 9:24 AM, Magdalini Ntetsika wrote: 
>>> > 
>>> > but I don't seem to get the same *system_rhs *when 
>>> > *assemble_system=false* as when *assemble_system=true*. To be more 
>>> > specific, it seems that there is something wrong with how 
>>> > constraints.distribute_local_to_global(local_rhs, local_dof_indices, 
>>> > system_rhs) creates the system_rhs from local_rhs and 
>>> > local_dofs_indices. I don't change anything about constraints etc 
>>> > throughout the time steps since I set those things up at the very 
>>> > beginning and then I just run several times without updating anything. 
>>> > Any idea about what the problem could be? 
>>>
>>> Magdalini, 
>>> it's not in the function you point to, but in your understanding. In 
>>> order to compute the correct values of the right hand side when you have 
>>> non-zero constraints (e.g., when you have nonzer boundary values), the 
>>> function you call needs to have access to the matrix elements. 
>>>
>>> The issue is a bit complicated to understand, but you may want to look 
>>> up the video lectures on the topic of Dirichlet boundary conditions. 
>>>
>>> As for your case, take a look at step-26, which deals with exactly this 
>>> problem by building the matrix only once, but then using it for several 
>>> time steps. 
>>>
>>> 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 dea...@googlegroups.com .
>> To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/a320a562-6800-4650-871f-83a5eff9b79a%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/a320a562-6800-4650-871f-83a5eff9b79a%40googlegroups.com?utm_medium=email_source=footer>
>> .
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/1964a904-4c6e-4ca6-8ba7-8c7a2c10769d%40googlegroups.com.

Re: [deal.II] Problem with "constraints.distribute_local_to_global(local_rhs, local_dof_indices, system_rhs)"

[Magdalini Ntetsika]

Hi Wolfgang,

thank you for the prompt reply. I see, I'll look up both. But I have a 
quick question: what if I have periodic boundary conditions all around my 
domain? It seems to be working in that case, even if the boundary values 
are nonzero, I get to have the correct system_rhs vector. I guess will need 
first to watch the topic on Dirichlet Boundary conditions to understand 
better, but it will be very helpful if I could get some feedback on that 
specific observation from you.

Best,
Magda

On Thursday, March 5, 2020 at 9:41:11 AM UTC-8, Wolfgang Bangerth wrote:
>
> On 3/5/20 9:24 AM, Magdalini Ntetsika wrote: 
> > 
> > but I don't seem to get the same *system_rhs *when 
> > *assemble_system=false* as when *assemble_system=true*. To be more 
> > specific, it seems that there is something wrong with how 
> > constraints.distribute_local_to_global(local_rhs, local_dof_indices, 
> > system_rhs) creates the system_rhs from local_rhs and 
> > local_dofs_indices. I don't change anything about constraints etc 
> > throughout the time steps since I set those things up at the very 
> > beginning and then I just run several times without updating anything. 
> > Any idea about what the problem could be? 
>
> Magdalini, 
> it's not in the function you point to, but in your understanding. In 
> order to compute the correct values of the right hand side when you have 
> non-zero constraints (e.g., when you have nonzer boundary values), the 
> function you call needs to have access to the matrix elements. 
>
> The issue is a bit complicated to understand, but you may want to look 
> up the video lectures on the topic of Dirichlet boundary conditions. 
>
> As for your case, take a look at step-26, which deals with exactly this 
> problem by building the matrix only once, but then using it for several 
> time steps. 
>
> 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/a320a562-6800-4650-871f-83a5eff9b79a%40googlegroups.com.

[deal.II] Problem with "constraints.distribute_local_to_global(local_rhs, local_dof_indices, system_rhs)"

[Magdalini Ntetsika]

Hi,

I want to assemble my system matrix only once since it doesn't change 
throughout the time steps. For that my code is similar to step-57:

if (assemble_matrix){
 constraints.distribute_local_to_global(local_matrix, local_rhs, 
local_dof_indices, system_matrix, system_rhs);
}
else{
constraints.distribute_local_to_global(local_rhs, local_dof_indices, 
system_rhs);
}

but I don't seem to get the same *system_rhs *when *assemble_system=false* 
as when *assemble_system=true*. To be more specific, it seems that there is 
something wrong with how constraints.distribute_local_to_global(local_rhs, 
local_dof_indices, system_rhs) creates the system_rhs from local_rhs and 
local_dofs_indices. I don't change anything about constraints etc 
throughout the time steps since I set those things up at the very beginning 
and then I just run several times without updating anything. Any idea about 
what the problem could be?

Best,
Magda

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/174a9f9e-491a-47f2-9e36-d737bca51d13%40googlegroups.com.

Re: [deal.II] Re: Problem with Block Preconditioner and Fully Periodic Boundary Conditions

[ntetsika]

Hi Daniel and Wolfgang,

When I apply periodic boundary conditions only on velocities at all the 
four sides of my unit square domain, then the rank of B (mxn) becomes 
rank(B) = n-1, which I think is as you say due to the fact that the 
pressure is up to a constant . However, when I try to impose periodic 
boundary conditions for both velocities AND pressure, then the rank of B 
becomes much less than the number of its columns and there comes my problem 
with inverting the B^T (diag M)^{-1} B. 

Any suggestions how to overcome this problem? I've tried to release some 
nodes on the boundaries (i.e. have periodic b.c. on all the sides except 
from 4 nodes - 1 free node per side) but still not working.

Thank you in advance!

Best,
Magda

On Thursday, January 17, 2019 at 7:17:25 PM UTC-5, Daniel Arndt wrote:
>
> Magda,
>
> [...] 
>>
> B has rank much less than the number of columns which is due to the fully 
>> periodic boundary conditions and the corresponding sparsity pattern. So I 
>> don't know how I could solve this issue. Any ideas?
>>
> Can you identify the corresponding eigenvector of B^T B? I would still not 
> be surprised if the pressure is only unique up to a constant implying that 
> a vector representing a constant pressure is such an eigenvector.
>
> Best,
> Daniel
>

-- 
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: Problem with Block Preconditioner and Fully Periodic Boundary Conditions

[ntetsika]

Hi Wolfgang,

that's exactly the problem, B has rank much less than the number of columns 
which is due to the fully periodic boundary conditions and the 
corresponding sparsity pattern. So I don't know how I could solve this 
issue. Any ideas?

Best,
Magda

On Wednesday, January 16, 2019 at 7:31:50 PM UTC-5, Wolfgang Bangerth wrote:
>
> On 1/16/19 10:19 AM, ntet...@berkeley.edu  wrote: 
> > 
> > thanks for the reply. I'm working on Immersed Boundary Method and 
> flagella 
> > swimming, and, more specifically, I'm trying to model a unit square 
> domain, 
> > which is infinitely replicated in two dimensions, governed by the 
> > incompressible time-dependent stokes equations for low Reynolds number. 
> For 
> > this reason I need periodic boundary conditions for velocity in both 
> > directions, but also for pressure. Therefore, for my unit square domain, 
> I 
> > need to have *pure* periodicity, i.e., 
> > 
> >   \bold{u}_left = \bold{u}_left and p_left = p_right 
> > 
> >   \bold{u}_bottom = \bold{u}_top and p_bottom = p_top 
> > 
> > Does that make sense to you? 
>
> Hm, I see. You are thinking of an infinite lattice of flagella at periodic 
> positions, all of which are swimming in synchrony. (They should make that 
> an 
> Olympic sport!) You then just cut out one lattice square. I guess the 
> boundary 
> conditions make sense then. 
>
>
> > My main problem with that is that when I'm trying to use the Block 
> > Preconditioner for the time-dependent Stokes (without convection), i.e. 
> when I 
> > need to form the {\tilde{S}}^{-1}  \approx [B^T (diag M)^{-1} B]^{-1} + 
> > \Delta{t}(\nu)M_p^{-1}, I have problems to invert the B^T (diag M)^{-1} 
> B - it 
> > says not invertible. Any ideas? 
>
> The term 
>B^T (diag M)^{-1} B 
> is not invertible if one of the following three conditions is true: 
> * B is an  n x m  matrix where m>n, i.e., it is not "tall and skinny" 
> * B is an  n x m  matrix with n<=m but its column rank is less than n. 
> * (diag M) is not invertible. 
>
> You can now work your way down this list and check which condition is the 
> problem. For example, I would check that all diagonal entries of M are in 
> fact 
> nonzero (you most likely want them to actually be positive as well). 
>
> 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: Problem with Block Preconditioner and Fully Periodic Boundary Conditions

[ntetsika]

Hi Wolfgang,

thanks for the reply. I'm working on Immersed Boundary Method and flagella 
swimming, and, more specifically, I'm trying to model a unit square domain, 
which is infinitely replicated in two dimensions, governed by the 
incompressible time-dependent stokes equations for low Reynolds number. For 
this reason I need periodic boundary conditions for velocity in both 
directions, but also for pressure. Therefore, for my unit square domain, I 
need to have *pure* periodicity, i.e.,

 \bold{u}_left = \bold{u}_left and p_left = p_right

 \bold{u}_bottom = \bold{u}_top and p_bottom = p_top

Does that make sense to you?

My main problem with that is that when I'm trying to use the Block 
Preconditioner for the time-dependent Stokes (without convection), i.e. 
when I need to form the {\tilde{S}}^{-1}  \approx [B^T (diag M)^{-1} 
B]^{-1} + \Delta{t}(\nu)M_p^{-1}, I have problems to invert the B^T (diag 
M)^{-1} B - it says not invertible. Any ideas?

Best,
Magda

On Tuesday, January 15, 2019 at 6:22:12 PM UTC-5, Wolfgang Bangerth wrote:
>
> On 1/15/19 10:02 AM, ntet...@berkeley.edu  wrote: 
> > thank you for replying. I have looked thoroughly step 45 but it's not of 
> much 
> > help for my case. I actually want to impose periodic boundary conditions 
> for 
> > both velocities and pressure (I need to have a representative periodic 
> cell). 
>
> Is this really what you want to do? Can you elaborate on your setup? 
>
> I'm asking this because generally, mathematically speaking, you can only 
> impose the *stress* at the boundary, not the pressure alone. Second, if 
> you 
> think for example about an infinite pipe and you cut a piece out, then at 
> the 
> left and right end of that pipe segment, you can enforce that the velocity 
> is 
> periodic -- but the pressure will *not* be periodic, as there will be a 
> pressure drop along the pipe which counters the friction forces. 
>
> 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.

[deal.II] Re: Problem with Block Preconditioner and Fully Periodic Boundary Conditions

[ntetsika]

Hi Daniel,

thank you for replying. I have looked thoroughly step 45 but it's not of 
much help for my case. I actually want to impose periodic boundary 
conditions for both velocities and pressure (I need to have a 
representative periodic cell). I used the following lines to set pressure 
to zero at the first node, but still gives me the same problem when I try 
to invert [B^T (diag M)^{-1} B] with cg solver for the purpose of Block 
Preconditioning - it says [B^T (diag M)^{-1} B] not invertible.

// constraint first pressure dof to be 0

std::vector boundary_dofs (dof_handler.n_dofs(), false);

DoFTools::extract_boundary_dofs (dof_handler, fe.component_mask(pressure), 
boundary_dofs);

const int first_boundary_dof = std::distance (boundary_dofs.begin(), 
std::find (boundary_dofs.begin(), boundary_dofs.end(), true));

constraints.add_line(first_boundary_dof);


I was thinking maybe the time-dependent stokes problem I'm trying to solve 
with fully periodic boundary conditions (velocities+pressure) is 
over-constrained and maybe that's why [B^T (diag M)^{-1} B] is not 
invertible. What do you think? Any ideas?


Thanks again for the help!


Best,

Magda

On Tuesday, January 15, 2019 at 11:05:38 AM UTC-5, Daniel Arndt wrote:
>
> Magda,
>
> [...]
>> DoFTools::make_periodicity_constraints(dof_handler, 2, 3, 1, 
>> constraints);
>> DoFTools::make_periodicity_constraints(dof_handler, 4, 1, 0, 
>> constraints);
>>
>>
> These two lines imply that you want to impose periodic boundary conditions 
> for all the components. Assuming the dof_handler variable is of type 
> FESystem and describes both velocity and pressure (as is the case in the 
> code-gallery example),
> this in particular holds true for the pressure. Is that what you are 
> intending?
> If so, you probably need another constraint to make the pressure unique.
> Otherwise, you might want to have a look at step-45 
>  where a 
> Stokes problem is solved with periodic boundary conditions on part of the 
> boundary.
> There only periodicity constraints for the velocity are stored in the 
> ConstraintMatrix object.
>
> Best,
> Daniel
>

-- 
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.

[deal.II] Problem with Block Preconditioner and Fully Periodic Boundary Conditions

[ntetsika]

Hi,


any chance I could get your help with that? I'm trying to use the code of 
@Jie-Cheng  for Block Preconditioner 
(https://github.com/dealii/code-gallery/tree/893bdb5b0a770758fce98859ab4d137f95e62a50/time_dependent_navier_stokes)
 
for a case of fully periodic boundary conditions, but then it gives me the 
error that [B^T (diag M)^{-1} B]^{-1} not invertible. Any ideas why this 
could happen or how could I fix it? I have increased the number of 
iterations for the CG solver and I also tried to use gmres, but none of 
them ways worked.


The way I impose the Periodic Boundary conditions is the following. I have 
a unit square domain with boundary ids 1,2,3 and 4 for the four sides.


std::vector > periodicity_vector;
std::vector > periodicity_vector1;

FullMatrix rotation_matrix(dim);
rotation_matrix[0][0]=1.;
rotation_matrix[1][1]=1.;

GridTools::collect_periodic_faces(triangulation, 2, 3, 1,
  periodicity_vector, Tensor<1, dim>(),
  rotation_matrix);
GridTools::collect_periodic_faces(triangulation, 4, 1, 0,
  periodicity_vector1, Tensor<1, dim>(),
  rotation_matrix);

triangulation.add_periodicity(periodicity_vector);
triangulation.add_periodicity(periodicity_vector1);

DoFTools::make_periodicity_constraints(dof_handler, 2, 3, 1, constraints);
DoFTools::make_periodicity_constraints(dof_handler, 4, 1, 0, constraints);


Thank you for your help in advance, I've been stuck with that a long time 
now and I really don't understand what's the problem.


Best,
Magda

-- 
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.

[deal.II] Incompressible fluid with fully periodic boundary conditions

[ntetsika]

Hi,

I'm a kind of new user of dealii and I'm trying to set up the boundary 
conditions for a domain [0, 1] X [0, 1] of incompressible fluid with fully 
periodic boundary conditions. I know that for an incompressible fluid with 
fully periodic boundary conditions (velocity both directions and pressure) 
I need to release one boundary DOF, otherwise the problem will be over 
constrained. 

Therefore, the question is:
Is there any way, after I have set up my periodic boundary conditions, to 
release one DOF? I was thinking it should be something in the constraint 
matrix but I wasn't sure. 

I have spent a lot of time to solve that problem but unsuccessfully. I 
would be great if you could give me some help. Thanks in advance!

Best,
Magda

-- 
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.