Adding new surface types

New ideas for code development

Adding new surface types

Postby Jaakko Leppänen » Thu Mar 25, 2010 5:03 am

Unlike MCNP and other Monte Carlo codes, the geometry model in Serpent doesn't allow the use of the union operator to construct cells. To compensate for this deficiency, there are a number of derived geometry types, consisting of multiple elementary surfaces, that can be used to construct most fuel and reactor geometries.

Adding new surface types is relatively easy as well. If your geometry needs a surface that is not supported by Serpent (see Sec. 3.2.1 in the User's Manual for available surface types), the easiest and fastest way for me to add it in the code is that you provide the description of two routines either in "pseudo-code" or, more preferably, in C.

The first routine determines whether an arbitrary point (x,y,z) is inside the boundary or not. For example, for a cube centered at (x0,y0,z0) with half-width r, this routine could be written as:

Code: Select all
        if ((x - x0 > r) || (x - x0 < -r) || (y - y0 > r) || (y - y0 < -r)|| (z - z0 > r) || (z - z0 < -r))
          return 0;
        else
          return 1;

Return value 0 indicates that the point is out and return value 1 that the point is in.

The second routine is a bit more complicated, and it determines the shortest possible distance to the boundary in an arbitrary direction (u,v,w). For the same cube this routine could be written as:

Code: Select all
       min = INFTY;

       if (u != 0.0)
          {
            if (((d = -(x - x0 - r)/u) > 0.0) && (d < min))
              min = d;
            if (((d = -(x - x0 + r)/u) > 0.0) && (d < min))
              min = d;
          }

        if (v != 0.0)
          {
            if (((d = -(y - y0 - r)/v) > 0.0) && (d < min))
              min = d;
            if (((d = -(y - y0 + r)/v) > 0.0) && (d < min))
              min = d;
          }

        if (w != 0.0)
          {
            if (((d = -(z - z0 - r)/w) > 0.0) && (d < min))
              min = d;
            if (((d = -(z - z0 + r)/w) > 0.0) && (d < min))
              min = d;
          }

        return min;

Notice that the returned distance doesn't have to be the actual shortest distance. It is sufficient that it is not longer than the shortest distance. In the cube example, if point (x,y,z) is outside the cube, the routine may actually find an intersection on one of the six planes that make up the cube, but also extend over the actual boundaries. This is not a problem, as long as none of the surfaces is ever missed.
- Jaakko
User avatar
Jaakko Leppänen
Site Admin
 
Posts: 1982
Joined: Thu Mar 18, 2010 10:43 pm
Location: Espoo, Finland

Planes parallel to an axis

Postby Andrei Fokau » Wed Aug 18, 2010 1:55 pm

Sometimes, it is not very desirable to calculate tangents for defining planes angle, and it is probably common enough to have planes parallel to an axis, so it would be useful to define special surface types for them:
Code: Select all
PPX  Y  Z  angle_between_surface_and_Y
PPY  X  Z  angle_between_surface_and_X
PPZ  X  Y  angle_between_surface_and_X
where the last one will certainly be most popular.
KTH Reactor Physics (Stockholm, Sweden) neutron.kth.se
Andrei Fokau
 
Posts: 77
Joined: Thu Mar 25, 2010 12:25 am
Location: KTH, Stockholm, Sweden

CYLPRISM, PADPRISM

Postby Andrei Fokau » Wed Sep 01, 2010 8:09 pm

Despite that pins (nests) are great for defining fuel elements, sometimes one needs to model the elements with more details. For such cases, it would be nice to have CYLPRISM surface type, which is a cylinder between z1 and z2. Additionally, it would also be great to have tube as PAD with no angles provided and PADPRISM:
Code: Select all
surf 1 cylprism  x0 y0 r z1 z2
surf 2 pad       x0 y0  r1 r2 # tube
surf 3 padprism  x0 y0  r1 r2  t1 t2  z1 z2

PADPRISM can work as tube prism and CYLPRISM:
Code: Select all
surf 4 padprism  x0 y0  r1 r2  0 360  z1 z2 # tube prism
surf 5 padprism  x0 y0  0  r   0 360  z1 z2 # cyl prism
KTH Reactor Physics (Stockholm, Sweden) neutron.kth.se
Andrei Fokau
 
Posts: 77
Joined: Thu Mar 25, 2010 12:25 am
Location: KTH, Stockholm, Sweden

Re: Adding new surface types

Postby Andrei Fokau » Thu Sep 02, 2010 8:53 pm

Jaakko, could you instruct me what I need to do in order to add a custom surface type to Serpent? I have tried to find out it myself, but got lost quickly due to not sufficient experience in C.

It may be more convenient to separate all surface-type-dependent logic and parameters to surface-type routines, so adding a new one would be just going through such "check-list".
KTH Reactor Physics (Stockholm, Sweden) neutron.kth.se
Andrei Fokau
 
Posts: 77
Joined: Thu Mar 25, 2010 12:25 am
Location: KTH, Stockholm, Sweden

Re: CYLPRISM, PADPRISM

Postby Andrei Fokau » Sat Sep 18, 2010 9:38 pm

Andrei Fokau wrote:Despite that pins (nests) are great for defining fuel elements, sometimes one needs to model the elements with more details. For such cases, it would be nice to have CYLPRISM surface type, which is a cylinder between z1 and z2.
Code: Select all
surf 1 cylprism  x0 y0 r z1 z2

I have tried to add CYLPRISM by myself and to my surprise discovered that the cylindric prism surface is available since a long time as CYL type:
Code: Select all
surf 1 cyl  x0 y0 r z1 z2

In analogy with CYL, it would be simpler to define hexagonal prisms as HEXXC and HEXYC types instead of having special types for them (namely, HEXXPRISM and HEXYPRISM):
Code: Select all
surf 1 hexxc  x0 y0 r r0 z1 z2
surf 2 hexyc  x0 y0 r r0 z1 z2

As I see in testsurface.c, it is extremely easy to implement, and HEXXPRISM and HEXYPRISM can actually be kept for version compatibility reason.

Jaakko, what do you think about it?
KTH Reactor Physics (Stockholm, Sweden) neutron.kth.se
Andrei Fokau
 
Posts: 77
Joined: Thu Mar 25, 2010 12:25 am
Location: KTH, Stockholm, Sweden

Re: Adding new surface types

Postby Jaakko Leppänen » Tue Sep 21, 2010 10:25 am

The hexagonal prism and infinite cylinder are not just for compatibility. The way the outer boundary is defined (single prism vs. cylinder and two axial surfaces) affects the implementation of the boundary conditions. In a prismatic surface, all boundaries, including the top and the bottom are reflective / periodic, while with the cylinder the repeated condition affects only that surface. (See Sec. 3.7 in the Manual).

It is true, though, that the same functionality could be implemented with a single surface type.
- Jaakko
User avatar
Jaakko Leppänen
Site Admin
 
Posts: 1982
Joined: Thu Mar 18, 2010 10:43 pm
Location: Espoo, Finland

Re: Adding new surface types

Postby Jaakko Leppänen » Tue Sep 21, 2010 10:40 am

In order to add new surface types on your own, you need to complete the following steps:

1) In locations.h, add an index for the new surface type (the list is found after #define SURFACE_TYPES). Increase the number of available surface types by one.

2) In processgeometry.c, add the surface name in vector surftp found at the beginning of the subroutine and the number of parameters in vector nparam. The first and last numbers are the minimum and maximum number of input parameters, respectively.

3) In testsurface.c, add source code that determines whether a point (x, y, z) is inside the surface or not. The existing definitions should be pretty instructive (see the spherical surface SURF_SPH, for example). The surface parameters are stored in DATA[ptr], DATA[ptr + 1], DATA[ptr + 2], etc...

4) In surfacedistance.c, add source code that calculates the minimum possible distance to the surface from point (x, y, z) in direction (u, v, w). The subroutine is structured very similar to testsurface.c.

I think that should do it. Let me know what happens if you decide to try this out.
- Jaakko
User avatar
Jaakko Leppänen
Site Admin
 
Posts: 1982
Joined: Thu Mar 18, 2010 10:43 pm
Location: Espoo, Finland


Return to Development

Who is online

Users browsing this forum: No registered users and 1 guest