Identifying dofs

Asked by Marcos Samudio

I am solving Navier-Stokes using both NSCoupled and StandardKE_Coupled, and I am interested in understanding the apply() method for wall functions in KE, so I am inspecting the matrix and vector modifications being performed by it (Wall.py file).

When solving coupled systems, how can I identify the dofs corresponding to one or another variable?(energy-dissipation in the KE case and velocity-pressure in NS). I suppose the unknowns vector must be ordered somehow, I just do not see how. Has this any relation with the weak formulation or something else in the problem formulation?

Thanks,

Marcos

Question information

Language:
English Edit question
Status:
Solved
For:
CBC.PDESys Edit question
Assignee:
No assignee Edit question
Solved by:
Mikael Mortensen
Solved:
Last query:
Last reply:
Revision history for this message
Mikael Mortensen (mikael-mortensen) said :
#1

Dear Marcos,

Den Mar 28, 2012 kl. 8:05 PM skrev Marcos Samudio:

> New question #191963 on CBC.PDESys:
> https://answers.launchpad.net/cbcpdesys/+question/191963
>
> I am solving Navier-Stokes using both NSCoupled and StandardKE_Coupled, and I am interested in understanding the apply() method for wall functions in KE, so I am inspecting the matrix and vector modifications being performed by it (Wall.py file).
>
> When solving coupled systems, how can I identify the dofs corresponding to one or another variable?(energy-dissipation in the KE case and velocity-pressure in NS). I suppose the unknowns vector must be ordered somehow, I just do not see how. Has this any relation with the weak formulation or something else in the problem formulation?

The ordering of dofs is a complicated issue and unfortunately I have not found a better way of implementing this boundary condition. You really need to understand how dofmaps work, look in the C++ definition of the dofmap class. For a single processor it is rather straight forward though. For a mixed k-epsilon space KE one has

K = FunctionSpace(mesh, 'CG', 1) # or something else
E = FunctionSpace(mesh, 'CG', 1) # has to be the same for routines in Wall.py to work
KE = K * E
ke = Function(KE)
kev = ke.vector()
N = K.dim()
kev[:N] are the dofs of k
kev[N:] are the dofs of epsilon
kev[0] and kev[N] are k and epsilon for the same node 0 in the dofmap.

The boundary condition reads

  epsilon = 2*nu*k/y**2

where y is the distance to the wall. So if node (dof) 1 is a node just inside the wall, then I want to set the equation

epsilon[1] = 2*nu*k[1]/y[1]**2 (1)

This is accomplished in the apply method of KEWall. The boundary condition needs to be set implicitly, which is the only stable way of imposing this boundary condition and a major reason for why k and epsilon needs to be coupled. The N+1'th row of the coefficient matrix A corresponding to (1) will look like

A[N+1, (0, -2*nu/y**2, 0, …. ,0, 1, 0, ...) ], where columns 1 and N+1 are nonzero

which is what you see is accomplished in apply.

Hope that helps

Mikael

>
> Thanks,
>
> Marcos
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.

Revision history for this message
Marcos Samudio (marcossamudio) said :
#2

Mikael,

Thanks for the fast response. I had already figured that out by reading the code, which is very clear and well written!

In apply(), the matrix is being manipulated, and in order to do so you have to know perfectly well how it is assembled. For instance, you 'move' to the epsilon dofs by doing self.N + dof_index, which is assuming epsilon dofs are in the end of the unknowns vector, as you said. This is what motivated my question in first place, as I did not know where the dofs where in the matrix! I will take a look at the dofmap class in C++ (although I am not a very experienced C++ programmer). Does it work for parallel processors?

May I ask further questions concerning the k-epsilon model? The Wall BCs are being applied both on k and epsilon. The additional KEWall BC is also modifying the values of epsilon at inner nodes (through the bnd_to_in mapping). Am I correct? Is that all?

Is there a way to apply the logarithmic velocity law close to the wall as well? Is there an implemented method for doing so? If not, I may try to write it and maybe even submit it. I am simulating a river bifurcation, so extreme mesh refinement is impossible!

Thanks again,

Marcos

Revision history for this message
Mikael Mortensen (mikael-mortensen) said :
#3

Den Mar 28, 2012 kl. 9:10 PM skrev Marcos Samudio:

> Question #191963 on CBC.PDESys changed:
> https://answers.launchpad.net/cbcpdesys/+question/191963
>
> Status: Answered => Open
>
> Marcos Samudio is still having a problem:
> Mikael,
>
> Thanks for the fast response. I had already figured that out by reading
> the code, which is very clear and well written!
>

Thanks:-)

> In apply(), the matrix is being manipulated, and in order to do so you
> have to know perfectly well how it is assembled. For instance, you
> 'move' to the epsilon dofs by doing self.N + dof_index, which is
> assuming epsilon dofs are in the end of the unknowns vector, as you
> said. This is what motivated my question in first place, as I did not
> know where the dofs where in the matrix! I will take a look at the
> dofmap class in C++ (although I am not a very experienced C++
> programmer). Does it work for parallel processors?
>

Unfortunately not. The dofmaps are more complicated in parallel and I haven't really had a chance to look into it.

> May I ask further questions concerning the k-epsilon model? The Wall BCs
> are being applied both on k and epsilon. The additional KEWall BC is
> also modifying the values of epsilon at inner nodes (through the
> bnd_to_in mapping). Am I correct? Is that all?

Yeah. It is possible to set only the inner dofs, the wall dofs or both. If you set both, then epsilon is constant in the entire element closest to the wall, which seems to be most stable. If you set only the inner dofs, then the wall dofs for epsilon don't enter the computations. Some CFD codes compute the epsilon value for the inner node and extrapolate a value on the wall. The most common, though, is to fix epsilon

>
> Is there a way to apply the logarithmic velocity law close to the wall
> as well? Is there an implemented method for doing so? If not, I may try
> to write it and maybe even submit it. I am simulating a river
> bifurcation, so extreme mesh refinement is impossible!

It has not been implemented, but you just need to approximate yplus and modify A more or less exactly as already done, just use the log law instead of the epsilon = 2nu*k/y**2 equation. If you add a new wall model I'm more than happy to add it to the repository.

Mikael

>
> Thanks again,
>
> Marcos
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.

Revision history for this message
Best Mikael Mortensen (mikael-mortensen) said :
#4

Den Mar 28, 2012 kl. 9:10 PM skrev Marcos Samudio:

> Question #191963 on CBC.PDESys changed:
> https://answers.launchpad.net/cbcpdesys/+question/191963
>
> Status: Answered => Open
>
> Marcos Samudio is still having a problem:
> Mikael,
>
> Thanks for the fast response. I had already figured that out by reading
> the code, which is very clear and well written!
>

Thanks:-)

> In apply(), the matrix is being manipulated, and in order to do so you
> have to know perfectly well how it is assembled. For instance, you
> 'move' to the epsilon dofs by doing self.N + dof_index, which is
> assuming epsilon dofs are in the end of the unknowns vector, as you
> said. This is what motivated my question in first place, as I did not
> know where the dofs where in the matrix! I will take a look at the
> dofmap class in C++ (although I am not a very experienced C++
> programmer). Does it work for parallel processors?
>

Unfortunately not. The dofmaps are more complicated in parallel and I haven't really had a chance to look into it.

> May I ask further questions concerning the k-epsilon model? The Wall BCs
> are being applied both on k and epsilon. The additional KEWall BC is
> also modifying the values of epsilon at inner nodes (through the
> bnd_to_in mapping). Am I correct? Is that all?

Yeah. It is possible to set only the inner dofs, the wall dofs or both. If you set both, then epsilon is constant in the entire element closest to the wall, which seems to be most stable. If you set only the inner dofs, then the wall dofs for epsilon don't enter the computations. Some CFD codes compute the epsilon value for the inner node and extrapolate a value on the wall. The most common, though, is to fix epsilon

>
> Is there a way to apply the logarithmic velocity law close to the wall
> as well? Is there an implemented method for doing so? If not, I may try
> to write it and maybe even submit it. I am simulating a river
> bifurcation, so extreme mesh refinement is impossible!

It has not been implemented, but you just need to approximate yplus and modify A more or less exactly as already done, just use the log law instead of the epsilon = 2nu*k/y**2 equation. If you add a new wall model I'm more than happy to add it to the repository.

Mikael

>
> Thanks again,
>
> Marcos
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.

Revision history for this message
Marcos Samudio (marcossamudio) said :
#5

Thanks Mikael Mortensen, that solved my question.