programmer's documentation
Basic examples

Local variables to be added

The following local variables need to be defined for the examples in this section:

integer iel, ielpdc, ikpdc
integer ilelt, nlelt
integer izone
double precision alpha, cosalp, sinalp, vit, ck1, ck2
double precision, dimension(:,:), pointer :: cvara_vel
integer, allocatable, dimension(:) :: lstelt

Initialization and finalization

The following initialization block needs to be added for the following examples:

! Allocate a temporary array for cells selection
allocate(lstelt(ncel))

At the end of the subroutine, it is recommended to deallocate the work array:

! Deallocate the temporary array
deallocate(lstelt)

In theory Fortran 95 deallocates locally-allocated arrays automatically, but deallocating arrays in a symetric manner to their allocation is good pratice, and avoids using a different logic for C and Fortran.

Map field array

! Map field arrays
call field_get_val_prev_v(ivarfl(iu), cvara_vel)

Body

Calculation and identification of the number of cells with imposed head loss term

if (iappel.eq.1.or.iappel.eq.2) then

2 calls:

Note
  • Do not use ckupdc in this section (it is defined with iappel = 3)
  • Use icepdc in this section only with (iappel = 2)

To be completed by the user: cell selection

Head loss examples

Example 1: No head loss (default)

ielpdc = 0

Example 2: Head losses defined by coordinates for zone

(4 <= x; 2 <= y <= 8)
No head losses else.

if (1.eq.0) then ! This test allows deactivating the example
izone = 0
ielpdc = 0
call getcel('X <= 6.0 and X >= 4.0 and Y >= 2.0 and'// &
'Y <= 8.0',nlelt,lstelt)

Generic subsection

For iappel = 1 , define ncepdp, the number of cells with head losses. This is valid for both examples above.

izone = izone + 1
do ilelt = 1, nlelt
iel = lstelt(ilelt)
izcpdc(iel) = izone
ielpdc = ielpdc + 1
if (iappel.eq.2) icepdc(ielpdc) = iel
enddo
endif

Defining the number of cells with head losses

if (iappel.eq.1) then
ncepdp = ielpdc
endif

Computing the head loss coefficient values

Third call, at each time step

Note:

Example 1: head losses in direction x

Diagonal tensor : Example of head losses in direction x

if (.true.) return ! (replace .true. with .false. or remove test to activate)
do ielpdc = 1, ncepdp
iel=icepdc(ielpdc)
vit = sqrt(cvara_vel(1,iel)**2 + cvara_vel(2,iel)**2 + cvara_vel(3,iel)**2)
ckupdc(ielpdc,1) = 10.d0*vit
ckupdc(ielpdc,2) = 0.d0*vit
ckupdc(ielpdc,3) = 0.d0*vit
enddo

Example 2: alpha = 45 degres

3x3 tensor: Example of head losses at alpha = 45 degres x,y direction x resists by ck1 and y by ck2
ck2 = 0 represents vanes as follows: in coordinate system x, y

orthogonal_reference_frame_sketch.gif
Orthogonal reference frame sketch
if (.true.) return ! (replace .true. with .false. or remove test to activate)
alpha = pi/4.d0
cosalp = cos(alpha)
sinalp = sin(alpha)
ck1 = 10.d0
ck2 = 0.d0

Filling the cells of the head loss matrix

do ielpdc = 1, ncepdp
iel = icepdc(ielpdc)
vit = sqrt(cvara_vel(1,iel)**2 + cvara_vel(2,iel)**2 + cvara_vel(3,iel)**2)
ckupdc(ielpdc,1) = (cosalp**2*ck1 + sinalp**2*ck2)*vit
ckupdc(ielpdc,2) = (sinalp**2*ck1 + cosalp**2*ck2)*vit
ckupdc(ielpdc,3) = 0.d0
ckupdc(ielpdc,4) = cosalp*sinalp*(-ck1+ck2)*vit
ckupdc(ielpdc,5) = 0.d0
ckupdc(ielpdc,6) = 0.d0
enddo
!-----------------------------------------------------------------------------
endif