Improvement of the P1-model

[Jonathan Gerardin]

yes, I think I have used source from beta 2.0. I don't find sources from 2.0-rc1. 
(remember Alexandre, the difference of line we encountered in raydom.f90)

[Alexandre Douce]

Jonathan, here http://code-saturne.org/forum/viewtopic.php?f=9&t=706 you have confirm that you use the 2.0-rc1. Perhaps the best way to avoid conflict between versions is that you delete the 2.0-beta2 (it is radical I know :) )

If you use the python installer, the sources of 2.0-rc1 are in the "install/ncs-2.0.0-rc1" directory. If you don't find them, you have just to reload the archive from the forum. Could you please retry your computation with the raypun.f90 from the 2.0-rc1 sources in the SRC directory? If raypun.f90 is unmodified,  the results with or without raypun.f90 must be identical (you can check that with the command sdiff -s listin1 listing2 for example). Of course this procedure must be applied for all your sources, not only raypun.f 90! I believe that your troubles of optical thickness will solve themselves.

However, in the scope of your development, don't forget to comment the verification of the adequacy of the optical thickness for the P1 model (in usray3.f90 and in raydom.f90).

Thanks, and keep us informed.

[Jonathan Gerardin]

It works! Thank you! I have eradicate the 2.0-beta2 version now! I can make some test during my development now.
 
I have commented the verification of adequacy of optical thickness in usray3.f90 (but it is write in comment "Arret en fin de pas de temps si épaisseur optique trop grande", but it is stop if the optical thickness is too low no? Same thing in raydom.f90 "Alerte si épaisseur optique trop grande"). 
 
Another thing, I have stocked the diffusion coefficient and the linear anisotropic coefficient in two new variables. But, I realise that in propce, you stocked all the properties of center of cells, so I have to put the diffusion and the linear anisotropic coefficient in propce too. How can I change the size of propce to add 2 new columns? (I can let the two new variables for the moment, we can see that after.)

[Jonathan Gerardin]

Is there a subroutine to calculate the symmetry? I have one center of the element, one center of a face which is a symetry, and I want to find the coordinate of the wall which is a real wall. So I have to calculate a specular symmetry from the center of the element to the center of the wall( symmetry).

[Alexandre Douce]

Good news!

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)

2) create your own arrays locally with an "allocate"

[Jonathan Gerardin]

For the moment, I use the second solution. Less good news, I think the calculation time will increase in an exponantiel type when the number of element will increase.

[Alexandre Douce]

Why? Do you have loop like:

Code: Select all

do iel1 = 1, ncel
  do iel2 = 1, ncel
    ...
  enddo
enddo

[Jonathan Gerardin]

To be sure,
 
surfbo(1,nfabor) gives the norm in x of the boundary face , same thing in y and z with ndim = 2 or 3.
ra(diipb(1,nfabor)) gives the norm in x of the normal vector of the boundary face, same thing in y and z with ndim = 2 or 3. I have the components of the normal vector with ra(diipb(ndim,nfabor)). It's right?
 
(and I don't know how I can verify the value of intermediate results...)
 
thank you

[Jonathan Gerardin]

I've got a loop when I search if the boundary face could see the element. So for each element, I have to calculate if each boundary face can see it, and to do this, I have to calculate the intercection between the vector MP (M center of the cell of the medium and P the center of the boundary face) and each boundary face, and I have to see if a boundray face could be between the element and the boundary face considered.

I think my explanation is not very clear... I post raypun.f90, there is a lot of comment on this party (at the end of the file), and I join an explanation of the calculation (in french... :) )

Previously Alexandre Douce wrote:

Why? Do you have loop like:

Code: Select all

do iel1 = 1, ncel
  do iel2 = 1, ncel
    ...
  enddo
enddo

[David Monfort]

You have a do-loop over a cell x a do-loop over the boundary faces x a while-loop over the boundary faces... I'm not surprise it takes somes times! However I don't have a better solution yet.

As regards the diipb array, it it related to the distance between I and I' (the projection of I to the normal of the boundary face). You'll find more information in the theory manual.