Improvement of the P1-model

[Jonathan Gerardin]

I don't find another solution for the moment. But maybe we could recalculate a rough mesh. We don't need a fine mesh to do the resolution, so maybe by combining some element in bigger element and boundary condition in bigger boundary condition, we can gain time calculation.

For diipb, I have understand that it is the distance between I and I' . And as I' is the projection of I to the normal of the boundary face, this gives me the 3 norm of the normal vector of the boundary face. I need the coordinate of the normal vector, so I can have this with diipb. no?

[David Monfort]

Actually, the face normal coordinates are given by:

Code: Select all

nx = surfbo(1, ifac)
ny = surfbo(2, ifac)
nz = surfbo(3, ifac)
surf = sqrt(nx*nx + ny*ny + nz*nz)

So you don't need the ra(idiipb) array to compute it. It represents the II' vector, which is orthogonal to the face normal.

But maybe did I misunderstand something in your question ?

[Jonathan Gerardin]

Hello,
I think my english is bad. ;)
I want the coordinate of a vector which is perpendicular to the boundary face ("la normale à la paroi" in french), so ra(diipb(ndim,nfabor)) seems to give me what I want.
My understanding of surfbo is that gives me the coordinate of 3 vectors ( in each ndim) which caracterise the boundary face. The boundary face is caracterised by 1 pt (cdgfbo) and 3 vectors (surfbo) and thanks to that, you can know how the boundary face is (orientation, surface...)
If surfbo doesn't give me that, i've got errors in my calculation. (I define the boundary face with 1 point and 2 vectors u and v which are a composition of the 3 vectors nx,ny,nz : u = nx + 1/2 ny et v = 1/2 ny + nz if I take your nx,ny,nz).
 
Is there a way to see intermediate result etc...
 
Thank you for your answers!

[Jonathan Gerardin]

I have manage to see intermediate result by writing them in listing. So I can see that surfbo gives the coordinate of a vector which is perpendicular to the surface. So I will search a new way to define my face.
So you have understand me, my english is not so bad! ;)

[David Monfort]

:)

Indeed the "normale à la paroi" is given by surfbo, and a vector perpendicular to this normal can be given by ra(idiipb) provided that the II' distance is not zero!

What you can do is using one of the vertices defining the boundary face and define your second vector with thanks to this vertice and the face centre. You can access the face vertices like this:

Code: Select all

do ii = ipnfbr(ifac), ipnfbr(ifac+1)-1
  inod = nodfbr(ii)
  xx = xyznod(1,inod)
  ...
enddo

This way, you can use all the boundary face vertices.

However, perhaps someone can come up with a bettter solution !

[Jonathan Gerardin]

Thank you for this advice. The calculation which determine if there is a wall between the boundary and the element works now. I've just got some calculation error which gives me for example Delta = 0.9999999999989 instead of 1. If Delta = 1, there is no obstacle, if Delta <1, there is an obstacle. So I have to put my condition to Delta < 0.99999.
 
The time calculation grows a lot due to that. But this calculation could be done only once, even if it's a unsteady calculation (except if the mesh is moving). So I will add  a condition to switch this calculation if it has already be done before. I will stock an array of integer visible(ncel,nfabor) which contain the value 1 or 0 if the cell see the boundary or not.
 
Now I will implement the symmetry calculation. When a cell see a boundary which type is symmetry, I'll find the real boundary by symmetry. After, I will be able to do the calculation of IDA for all the cell.
 
I will continue to give you informations on my progress (and my difficulties too. sorry...)

[Jonathan Gerardin]

Hello,

I used a new approach of my problem. Before, for each cell, I go to each boundary to make the IDA calculation. But we see huge time calculation when the number of boundary face is high. The calculation was done in a lot of direction. So, I will try a new approach. For each cell, I go to a fixed number of direction, which reduce a lot the time calculation.

I'm close to finish the program of IDA before test it but I need a new information. I've got a point and I want to know the cell where this point is. (I need the intensity, the radiative heat flux, the absorption coefficient and the scattering coefficient of this cell, so I need the number of the cell) How can I do that? I could with the cell vertices but as I understand, I can find the boundary face vertices, or maybe the internal face vertices, but I don't know if I can find the cell vertices. And search that by scanning all the internal face vertices seems to be long (I haven't try it, but if it's is the unique solution, I will try that)

Another question,  about propce :

Alexandre Douce :

You have two possibilities for your new arrays:

1) modified the number of variables contained in PROPCE: in oder to do that you have to modifed the varpos.f90 (line 2066) subroutine, and the include file radiat.h (lines 76 and 87)

So I add in varpos.f90 idif(iphas) and iani(iphas) between iabs(iphas) and iemi(iphas) to add the scattering coefficient and the anisotropic linear coefficient for each cell.

In radiat.h, i add idif(1+nclcpm) and iani(1+nclcpm) (and idif and iani in the common).

But for paraview for example, I need modify cs_gui_radiative_transfer.c too? In order to put the name of the scattering coefficient and the anisotropic linear coefficient, and make the good reorganisation of the cell properties.

After, I also want you to show me how to modify the python code to modify the gui. But we will see that after the implementation of the IDA.

thank you for your help!

[Jonathan Gerardin]

Hello,
 
I've got a new question. I have a memory problem. When I put some write(nfecra,*) in my code, my results are good, but when I comment those write(nfecra,*), my result are wrong. It's a memory problem. So I have to compile my .f90 with some option, like -fbound-check.
How can I add some option of compilation?
 
Thank you!

[Yvan Fournier]

Hello,
Both -fbounds-check and running under Valgrind regularly are useful. If you are able to reproduce the bug when compiled in debug mode, things should be easier: you'll get line numbers under Valgrind, and you'll automatically have -fbounds-check enabled using GCC.
To build in debug mode, simply add "--enable-debug" to the configure options for Code_Saturne, and re-run configure && make && make install.
Otherwise, to change compiler options, edit the "bin/cs_auto_flags.sh" file, and re-run configure && make && make install.
Best regards,

[Jonathan Gerardin]

Hello,
 
Thank you for your help, I manage to see where the problem was. (I use a table in a loop but this table was allocated before the loop, I allocate/deallocate it in the loop and it works)
 
Now, I've got a new problem with the intensity. I need the intensity to make a calculation for the IDA (or MDA). The intensity is in propce(iel,ipproc(ilumin)). But when I call it, all the values are equal to 0.0d0 . I call it at the end of raypun.f90 . The intensity is not calculated yet in raypun.f90?
Other question, to use the MDA, I need to make a P1 calculation with the wall at 0K instead of the real value of the wall. The temperature of the wall is stocked in tparoi(nfabor)? If I modify it in raydom.f90 before the call of raypun, it will works?
 
Thank you for your answers.
 
Jonathan