(EVAC) Question about code

[springcb]

Following is 'line 279' of read.f90.

           EVACUATION_Z_OFFSET(NM) = EVAC_Z_OFFSET

I think this line has no meaning since EVAC_Z_OFFSET is arbitrary.

At subroutine READ_VENT,
FDS makes two more VENTs for EVAC.

             IF (EVACUATION_ONLY(NM))           N_VENT = N_VENT + 2

However, later
These two VENTs are rejected by following "Evacuation criteria"

             IF (.NOT.EVACUATION .AND. EVACUATION_ONLY(NM)) REJECT_VENT=.TRUE.

As a result, adding the VENTs has no meaning.
How about my thought ?

at line 5892 of evac.f90

            EVAC_Z_MIN =  HUGE(EVAC_Z_MIN)
            EVAC_Z_MAX = -HUGE(EVAC_Z_MIN)
            DO nm = 1, NMESHES
                        MFF=>MESHES(nm)
                        EVAC_Z_MIN = MIN(EVAC_Z_MIN,REAL(MFF%ZS,FB))
                        EVAC_Z_MAX = MAX(EVAC_Z_MAX,REAL(MFF%ZF,FB))
            END DO

=>

            EVAC_Z_MIN =  HUGE(EVAC_Z_MIN)
            EVAC_Z_MAX = -HUGE(EVAC_Z_MIN)
            DO nm = 1, NMESHES
                        MFF=>MESHES(nm)
                        IF (.NOT.EVACUATION_GRID(nm)) CYCLE
                        EVAC_Z_MIN = MIN(EVAC_Z_MIN,REAL(MFF%ZS,FB))
                        EVAC_Z_MAX = MAX(EVAC_Z_MAX,REAL(MFF%ZF,FB))
            END DO

At line 6014~6016 In evac.f90

                        END SELECT
            END IF
            EXIT SS_Loop

=>

                        END SELECT
                        EXIT SS_Loop
            END IF

How about my thought ?

Thanks in advance.

[springcb]

in SUBROUTINE CHECK_EVAC_NODES

"N_MESH_TMP" is not deallocated.

[springcb]

399 line in evac.f90

I_HERDING_TYPE is not READ(EB) but INTEGER.

[springcb]

6401 line in evac.f90

    IF (.NOT.Is_Solid) THEN
        II = FLOOR( CELLSI(FLOOR((x_tmp(1)-XS)*RDXINT)) + 1.0_EB )
        JJ = FLOOR( CELLSJ(FLOOR((y_tmp(1)-YS)*RDYINT)) + 1.0_EB )
        I_OBST = OBST_INDEX_C(CELL_INDEX(II,JJ,KK))
        Is_Solid = (Is_Solid .OR. (SOLID(CELL_INDEX(II,JJ,KK)) .AND. .NOT. OBSTRUCTION(I_OBST)%HIDDEN))
        II = FLOOR( CELLSI(FLOOR((x_tmp(3)-XS)*RDXINT)) + 1.0_EB )
        JJ = FLOOR( CELLSJ(FLOOR((y_tmp(3)-YS)*RDYINT)) + 1.0_EB )
        Is_Solid = (Is_Solid .OR. (SOLID(CELL_INDEX(II,JJ,KK)) .AND. .NOT. OBSTRUCTION(I_OBST)%HIDDEN))
        II = FLOOR( CELLSI(FLOOR((x_tmp(2)-XS)*RDXINT)) + 1.0_EB )
        JJ = FLOOR( CELLSJ(FLOOR((y_tmp(2)-YS)*RDYINT)) + 1.0_EB )
        Is_Solid = (Is_Solid .OR. (SOLID(CELL_INDEX(II,JJ,KK)) .AND. .NOT. OBSTRUCTION(I_OBST)%HIDDEN))
    END IF

=>

    IF (.NOT.Is_Solid) THEN
        II = FLOOR( CELLSI(FLOOR((x_tmp(1)-XS)*RDXINT)) + 1.0_EB )
        JJ = FLOOR( CELLSJ(FLOOR((y_tmp(1)-YS)*RDYINT)) + 1.0_EB )
        I_OBST = OBST_INDEX_C(CELL_INDEX(II,JJ,KK))
        Is_Solid = (Is_Solid .OR. (SOLID(CELL_INDEX(II,JJ,KK)) .AND. .NOT. OBSTRUCTION(I_OBST)%HIDDEN))
        II = FLOOR( CELLSI(FLOOR((x_tmp(3)-XS)*RDXINT)) + 1.0_EB )
        JJ = FLOOR( CELLSJ(FLOOR((y_tmp(3)-YS)*RDYINT)) + 1.0_EB )
        I_OBST = OBST_INDEX_C(CELL_INDEX(II,JJ,KK))
        Is_Solid = (Is_Solid .OR. (SOLID(CELL_INDEX(II,JJ,KK)) .AND. .NOT. OBSTRUCTION(I_OBST)%HIDDEN))
        II = FLOOR( CELLSI(FLOOR((x_tmp(2)-XS)*RDXINT)) + 1.0_EB )
        JJ = FLOOR( CELLSJ(FLOOR((y_tmp(2)-YS)*RDYINT)) + 1.0_EB )
        I_OBST = OBST_INDEX_C(CELL_INDEX(II,JJ,KK))
        Is_Solid = (Is_Solid .OR. (SOLID(CELL_INDEX(II,JJ,KK)) .AND. .NOT. OBSTRUCTION(I_OBST)%HIDDEN))
    END IF

[rmcdermo]

Timo,
Please address issues on GitHub. We want to encourage use of this issue tracker. Thanks!

[tkorhon1]

Hi Randy,
Yes, I'll start using this now. But I filter "assigned to me", other issues I do not read,
not even the subject lines. Too many to read...

Timo

[tkorhon1]

And now the SpringTiger comments:

1. This is already modified in my local copy. Need to "commit" it sometime. This
   does not affect the results, it is just a line of code, that is not needed anywere,
   i.e., EVACUATION_Z_OFFSET(NM) is not used for these kind of evacuation
   meshes. Just the main evacuation meshes (and automatic STRS main
   evacuation meshes) use this information and it is set correctly for those
   meshes.

2. This I commented also at the discussion forum (google). The 2 additional
   vents are the (automatic) -z and +z boundary conditions on the evacuation
   meshes (mirror boundaries).

3. This is also commented in the discussion forum, no need to chage the code.
   This has just the effect how agents are shown in the z direction, does not affect
   the evacuation results. It just affects how agents are shown for visualization.
   And my choise is that agents that would go "out of bounds in Smokeview"
   are shown at the highest (or lowest) z coordinate that Smokeview shows
   (the max/min of the z of the all meshes).

4. This was a bug and it is corrected in my local copy. I should "commit"
   the changes. But I have 1.5 months updates to source code to be checked.

And then: "N_MESH_TMP" is not deallocated.: Thanks, now this change
is also in my local copy.

399 line in evac.f90: Thanks, now it is a couple of lines earlier, where INTEGERS are.

6401 line in evac.f90: I do not get this. I do not see any difference between the
two pieces of the code.

[springcb]

6401 line in evac.f90

1. "I_OBST" is updated just once.
2. "I_OBST" is updated two more times.

Thanks for your reply.

[tkorhon1]

Thanks!

Now I have already commited these corrections to the GitHub.

[springcb]

(User's Guide) page 74

........ NO_EVACUATION is set to .FALSE. on the MISC namelist.
=>
........ NO_EVACUATION is set to .TRUE. on the MISC namelist.

[springcb]

(User's Guide) page 80

The value of DIA_MEAN or if not given (not needed for all distributions)
then the distribution mean is used to scale the other body dimensions like
D_TORSO_MEAN.

=> Is this a complete sentence ?

[springcb]

(User's Guide) page 81
default value of FED_DOOR_CRIT : 0.000001 ?

That's not the case with source code.

[springcb]

(User's Guide) page 91

=>
The default value of FAC_SPEED is 1.0 but 0.6.

[tkorhon1]

Thanks!

Now I have started to do the above things:

1. NO_EVACUATION is set to .TRUE. now.
2. Yes, the sentence is not correct, something is missing.

The value of DIA_MEAN is used to scale the other body
dimensions like D_TORSO_MEAN}. If DIA_MEAN is not
given (not needed for all distributions) then the distribution mean
is used to scale the other body dimensions.

1.   !Timo: FED_DOOR_CRIT = 0.000001_EB ! Which doors are 'smoke free' Evac <= 2.2.1
   

   FED_DOOR_CRIT = -100.0_EB ! Which doors are 'smoke free' Evac >= 2.2.2
   Now this is updated (was already on my todo list...)

2. Thanks, that is now also corrected. They are now in my own LaTeX2e source code. I need to recompile and put the pdf at the correct place. I'm am not doing that right now, because I need to update the V&V results also. I have rerun my V&V cases (not all "monte carlo" type data), I just need to replot the figures. So, I should update the gide to FDS 6.4.0 version (and evac 2.5.1).

Timo

[springcb]

If one sets the area of &EVHO as same as that of &EVAC by mistake, FDS+EVAC falls into an infinite loop.

I suggest that 6346~6357 lines move to 6381 line position in evac.f90.

[tkorhon1]

Thanks, it is now moved a couple of lines later. I am not going to put the new source
code to Github right now. This correction will come to Github next time, I need to
change some other parts of the code. This correction is not crucial one, because
it can be fixed by changing user inputs.

TimoK

[springcb]

evac.f90

 IF (TRIM(EVAC_EXITS(N)%ID) /= TRIM(EMESH_EXITS(JJ_NOW)%ID)) THEN
      WRITE(LU_ERR,*)'*** ERROR IN FINE_PREFERRED_DIRECTION'
      WRITE(LU_ERR,*)'*** ',TRIM(EVAC_DOORS(N)%ID), ' is not ',TRIM(EMESH_EXITS(JJ_NOW)%ID)

->

 IF (TRIM(EVAC_EXITS(N)%ID) /= TRIM(EMESH_EXITS(JJ_NOW)%ID)) THEN
     WRITE(LU_ERR,*)'*** ERROR IN FIND_PREFERRED_DIRECTION'
     WRITE(LU_ERR,*)'*** ',TRIM(EVAC_EXITS(N)%ID), ' is not ',TRIM(EMESH_EXITS(JJ_NOW)%ID)

FINE -> FIND
EVAC_DOORS -> EVAC_EXITS

[tkorhon1]

Thanks, I made the change today. I was on a parental leave 3 months, now back at the office.

Wbr,
Timo

[godisreal]

The speed-dependent social force is implemented as below (evac.f90). However, I do not quite get the idea. It seems that A_WALL is decreased for some reason. Does anyone know how the force is adjusted? So if the actual speed is less than desired speed, the social force and wall force will decrease. Is that right?

  ! =======================================================
  ! Speed dependent social force
  ! =======================================================
  IF (C_HAWK >= 0.0_EB .AND. HR%I_DoorAlgo==1) THEN
     HR_A = HR%DoseCrit1
  END IF
  HR_A =  HR_A*MAX(0.5_EB,(SQRT(HR%U**2+HR%V**2)/HR%SPEED))
  A_WALL = MIN(A_WALL, FAC_A_WALL*HR_A)

[godisreal]

Excuse me, I have a question.
The units of the forces are set to be 'N/m' as below in data.f90, but I do not get why it is 'N/m', not 'N'? I mean the unit of force should be Newton, isn't it?

OUTPUT_QUANTITY(245)%NAME = 'HUMAN_CONTACT_LINEFORCE'
OUTPUT_QUANTITY(245)%UNITS = 'N/m'
OUTPUT_QUANTITY(245)%SHORT_NAME = 'force_c'

OUTPUT_QUANTITY(246)%NAME = 'HUMAN_TOTAL_LINEFORCE'
OUTPUT_QUANTITY(246)%UNITS = 'N/m'
OUTPUT_QUANTITY(246)%SHORT_NAME = 'force_tot'

[tkorhon1]

Well, the outputs (labels, units, etc. defined in the data.f90) are more like pressure outputs. Think that the agent is a two dimensional circle. If the agent is in a really dense human crowd, the other persons are more or less making the agent to feel large pressure. Units of pressure is N/m2. But now we have just a two dimensional circle ==> N/m is a better unit.

Speed dependent social force:
Forget next three lines:
IF (C_HAWK >= 0.0_EB .AND. HR%I_DoorAlgo==1) THEN
HR_A = HR%DoseCrit1
END IF
This was some simple hawk-dove game that on person (Simo Heliövaara) did. So, he chekced some things using concepts coming from system analysis branch of applied mathematics (including game theory).

FAC_A_WALL = 1.0 by default. And:
A_WALL = FAC_A_WALLHR%A
So, this user input can be used to change the wall - agent social force (FAC_B_WALL also). But, by default, the human-human and human-wall social forces are the same.
The line:
A_WALL = MIN(A_WALL, FAC_A_WALLHR_A)
When the agent is walking faster than its preferred walking speed v0_i ==> The social force strength is increased for human-human intercation ==> agents try to be futher away from each others, because they are going fast (and they need more space for their feets, e.g., longer step lenght). But for the human-wall interaction the social force is not increased. So, the agents are walking as close to the walls as normally (sideways, motive force tries to keep the agents from colliding "head-on" the walls).

If the agent is walking slower than the desired walking speed v0_i, the human-human and human-wall social forces are both made weaker. So, the agent can go closer to the (side) walls.

Why not to change the A_WALL larger, when going fast? Well, you could try. I did not find this very usuful, because usually the egress routes are crowded and, typically, the walking speeds are less than the preferred (unimpeded) walking speeds. So, if the agents are going faster than v0_i ==> there are just few persons ==> things do not matter too much how close these are to walls. "do not matter" = the total evacuation time "sense". But things could be better looking in Smokeview, e.g., more natural looking.
Note: A disabled person might have quite low v0_i ==> can have larger speed in a human crowd (due to the social force push by the other persons). But this disabled person is not going to walk too fast, mainly the person just goes with the "flow" ==> no need to change the disabled person - wall interaction.

If the agent speed is low (or zero, i.e., standing) => social force is 50% of its "normal speed value". Why? Well, you need some social force so that in a crowd the agents do not touch each others. What should be its strength? I just played with the values and this choice seemed to be good enough. See the FDS+Users guide section "5.3 Human Parameter Sensitivity".

TimoK

[godisreal]

Hi, Timo. Thanks. The speed-dependent social force is a good thing, and it adjusts the force in different occasions.

Also, I checked the 'HUMAN_CONTACT_LINEFORCE' and 'HUMAN_TOTAL_LINEFORCE' in evac.f90
'HUMAN_CONTACT_LINEFORCE' -> HR%Sumforces -> CONTACT_F
'HUMAN_TOTAL_LINEFORCE' -> HR%Sumforces2 -> SOCIAL_F + CONTACT_F

          ! Save the forces for the loop evac_move_loop (next time step)
          HR%F_X = P2P_U
          HR%F_Y = P2P_V
          HR%SUMFORCES  = CONTACT_F
          HR%SUMFORCES2 = SOCIAL_F + CONTACT_F
          HR%TORQUE = P2P_TORQUE

             CASE(245)  ! CONTACT_LINEFORCE, Human Pressure: contact forces
                QP(NPP,NN) = REAL(HR%SumForces ,FB)
             CASE(246)  ! TOTAL_LINEFORCE, Human Pressure2: contact + social
                QP(NPP,NN) = REAL(HR%SumForces2 ,FB)

Please double check. Thanks.

What do fed_max_door and k_ave_door represent in the code (evac.f90)? They seem to be important indicators related to hazard conditions, do they?

For example, the following code is for visible doors.

IF (FED_DOOR_CRIT >= TWO_EPSILON_EB) THEN
    L2_tmp = FED_max_Door(i) * SQRT((x1_old-x_o)**2 + (y1_old-y_o)**2)/Speed
ELSE
    L2_tmp = K_ave_Door(i)
END IF

Also, I do not quite get rad_flux in evac.f90. Is it the heat radiation represented by HUMAN_GRID(i,j)%RADFLUX? However, I thought radiation should be a vector???

OK. Here is some update about RADFLUX in evac.f90:
The direction of radiation flux is determined by the gradient of HUMAN_GRID(i,j)%TMP_G in a 2D space. The magtitude is given by HUMAN_GRID(i,j)%RADFLUX. So it is feasible to reconstruct the vector in a 2D evacuation mesh by using both together.

BTW, where is Timo? He does not like Github?

[tkorhon1]

Yes, the output CASE(245) is the contact forces (per metre of the "body circle") and CASE(246) is the total force, i.e, it includes the social force also.

Now the default in FDS+Evac is the K_ave_door criterion, i.e., the visibility (or smoke density). So, by, default, FED_DOOR_CRIT < 0.0. Some lines copied from evac.f90 (the present one in GitHub):

Line 15829 about:
K_ave_door = 0.0000001_EB*ABS(FED_DOOR_CRIT)
FED_max_Door = 0.0_EB

Make these zero before a sum loop. The K_ave_door is not put zero, because there
is later a line like "something = xxx/K_ave_door" (visibility calculation). If extinction
coefficient is zero => visibility is infinite. So, this is just to avoid division by zero.
So, the K_ave_door is the extinction factor, see, e.g., line 15484:
ave_K = EVAC_MASS_EXTINCTION_COEFF1.0E-6_EBM%HUMAN_GRID(i_r1,j_r1)%SOOT_DENS

So, for the visibility criterion, the average K is calculated/estimated for the different doors. For the FED criterion (FED_DOOR_CRIT >0) the FED dose is estimated.

The radflux is the integrated radiant intensity or something like that. It is not directly related to the radiation coming from some direction. But it tells something about average radiation over all angles. If it is large => probably too hot for a human. Well, might be that gas temperature would tell this information also.

TimoK

[tkorhon1]

Well, accidentally closed this one (I pushed a wrong button).

TimoK

[godisreal]

Hi, Timo. How to use the falling down model in evac module?
There is a falling-down model in evac.f90, but there is no description about it in the user guide. Is the model usable? I hope to use it in a simulation of doorway scenario. Thanks.

[tkorhon1]

Well, the falling down model was something that a student was checking. But that work was not finished. So, the model might not work anymore. It was a simple model. If I figure out correctly, an agent falls down, if the contact force was too high. There was some probability for this falling down also, if I remember correcty. And if you walked over a fallen down agent you could fall down. All of these features might be in the code that is inside evac.f90 or they are not. Some features might have just been something that we thought that we could implement.

FDS+Evac is treating the agent dynamics strictly in two dimensions, so falling down is not easy. One can not make any physical model, just some statistical model or some rule based model.

[godisreal]

Hi, Timo, I have tested the falling down model in a doorway example and it works fine generally.
An important issue refers to reducing the strength of social force so that contact force becomes effective when agents interact with each other at the bottleneck. If I use the default value A=2000 and B=0.08 the falling down scenario does not occur. So I reduce the value of A in the test case, and someone exactly falls down at the narrow doorway.

I read part of the code, and as you said, the model implementation is not complete.
One question about the code is

          L_FALLDOWN_THIS_TIME_STEP = .FALSE.
          IF (.NOT. L_FALLEN_DOWN .AND. TAU_FALL_DOWN >= 0.0_EB .AND. T-HR%T_CheckFallDown > 0.1_EB) THEN
             ! Check later if this agent falls down, do not do it at every time step
             HR%T_CheckFallDown = T ! Do the check at this time step, save the check time
             CALL RANDOM_NUMBER(RN_REAL)
             RN = REAL(RN_REAL,EB)
             ! P[falls down during Delta_t] = 1 - exp(-lambda*Delta_t), now Delta_t = 0.1 s
             IF (1.0_EB - EXP(-PROB_FALL_DOWN*0.1_EB) > RN) L_FALLDOWN_THIS_TIME_STEP = .TRUE.
          END IF

I think it should be HR%T_CheckFallDown = T + TAU_FALL_DOWN
What do you think? It is just like TAU_CHANGE_V0 or TAU_CHANGE_DOOR
Your comment is much appreciated.

OK. So the probablity is computed only for 0.1_EB time interval, and TAU_FALL_DOWN is not used there. Please omit my previous comment.
The probablity is used to deal with D_OVERLAP_FALL, and thus it is the probability of falling down when an agent walks over an existing fallen-down agent.

IF (TAU_FALL_DOWN >= 0.0_EB .AND. TIM_DIST <= R_TMP(III)+R_TMP(JJJ)-D_OVERLAP_FALL) THEN
                         ! Circles are touching each others more than D_OVERLAP_FALL => Fall down
                         IF (L_HRE_FALLEN_DOWN) COSPHIFAC = 0.0_EB ! On top of a fallen agent, no social force anymore
                         IF (T > HR%TDET .AND. L_FALLDOWN_THIS_TIME_STEP .AND. L_HRE_FALLEN_DOWN .AND. .NOT.L_FALLEN_DOWN) THEN
                            L_FALLEN_DOWN = .TRUE.
                            HR%T_FallenDown = T
                            HR%Angle_FallenDown = HR%angle
                            HR%SizeFac_FallenDown = 0.0_EB
                            N_DEAD = N_DEAD+1
                         END IF
                      END IF 

[tkorhon1]

I think it should be HR%T_CheckFallDown = T + TAU_FALL_DOWN
What do you think? It is just like TAU_CHANGE_V0 or TAU_CHANGE_DOOR
Your comment is much appreciated.

No, it should be as it is. HR%T_CheckFallDown is the time, when the agent (last) checked the falling down. See: " .AND. T-HR%T_CheckFallDown > 0.1_EB)". So, the IF statement checks if there has been more than 0.1s from the last check => recheck.

Yes, the D_OVERLAP_FALL is a measure when an agent falls down.

And note, that the A in the social force are speed dependent:

      ! =======================================================
      ! Speed dependent social force
      ! =======================================================
      HR_A =  HR_A*MAX(0.5_EB,(SQRT(HR%U**2+HR%V**2)/HR%SPEED))
      A_WALL = MIN(A_WALL, FAC_A_WALL*HR_A)

You could also make the B speed dependent and change the A and B speed reduction factors so
that the social force is smaller in bottlenecks. But you should check how these affect the specific
flows though doors and corridors (see the validation in FDS+Evac Guide).

TimoK