Speed up UCB urban model
Opened this issue · 19 comments
I tried with this input and it works, but much slower than using Hamada.
https://drive.google.com/file/d/119VjwXsv0s4xQr8MhSc7dtfCBecbMNPD/view?usp=sharing
Thanks Dwi. Are you using this building_fuel_models.csv file?
Yes, I am
One thing I noticed right away - in elmfire_spread_rate.f90 line 407 in the subroutine UMD_UCB_BLDG_SPREAD, this line causes division by zero if V_S is zero:
C%LOW = AMIN1((C%VELOCITY_DMS+C%VBACK)/2/V_S,10.0)
I'll change it to something like:
IF (V_S .GT. 1E-4) THEN
C%LOW = AMIN1((C%VELOCITY_DMS+C%VBACK)/2/V_S,10.0)
ELSE
C%LOW = 1.0
ENDIF
OK, I created a new branch urban-spread-optimization and you can see changes here:
main...urban-spread-optimization
These are simple optimizations - A*A instead of A**2, multiplying by reciprocal cellsize instead of dividing by cellsize, etc. I doubt it's going to bring down the run time significantly though.
The computational challenge is that the linked list LB contains thousands of nodes and we're looping over all of those nodes multiple times. Do you think it's possible to limit the number of LB nodes that are looped over using a distance criterion or similar?
The range of LB is only 100m radius, so its 4-ish cells (16 cells area in total). What I know is we can do active cells, I never try it before, but its like only a limited number of cells within the range is going to be calculated.
OK, it seems the biggest problem is probably with the LB list. I did some basic profiling yesterday and saw that the LB started off with a few nodes but quickly grew to 10,000+ nodes. So we're looping over 10,000+ LB nodes for every cell with fuel model 91.
Are you running in Visual Studio for debugging? If not, I highly recommend it. If you get stuck I'm happy to help you get up and running. It will save a lot of time in the long run.
Yes, I use Visual Studio now. Maria helped me set it up.
Excellent. Are you seeing any improvement in run time with the urban-spread-optimization branch:
main...urban-spread-optimization
If the changes look good to you I'd like to merge it with main before going further.
The changes work fine. The time improvement, however, is only around 2%. This is from 10 repetitions of running. I use one of the tutorial (01 - case)
That's about what I expected.
To get a significant speed improvement, a change in the algorithm is going to be needed. Rather than a single global building list that may contain tens of thousands of nodes, I think we're going to have to create a local building list for each "tagged" cell that contains at most ~10 ish building nodes. Right now, lots of cpu time is spent analyzing nodes in the LB list where TARGET_R_METERS is > 1 km.
BTW, I'm going to merge the urban-spread-optimization branch into main
Yes, I agree with that.
OK - are you familiar enough with the code at this point to create a new branch and make an initial attempt at it?
Yes, I will do it
I have created a pull request for the first attempt of speeding up the simulation. It makes it 5 times faster. For now, it is just using IF command in which if the burning cell is within the predefined neighborhood, then the subsequent commands are executed, otherwise, they are skipped. It means it still goes through all the burning cells, its just that only a limited number of these cells will be thoroughly calculated.
I think this can be made faster by limiting the number of elements in the loop, so the time to execute the IF command will be taken back. Still thinking of a way of doing this
This is the pseudocode in MATLAB for the next speed up. What is the best wayto implement it in FORTRAN.
allIX = [LIST_BURNED.IX]; %make an array of all IX coordinate of burning cells and its corresponding index in LIST_BURNED
IX_INCLUDED = (allIX >= IX_LOW_BORDER && allIX <= IX_HIGH_BORDER); % Make all the values outside the limit to be 0 and 1 otherwise
ind = find(IX_INCLUDED); % Find the index number that have 1 value for IX_INCLUDED
for i = 1:numel(ind) % iterate only for the number element in the ind variable
DEL_X = C.IX - LIST_BURNED{ind(i)}.IX % only the burning cell that has index stored in ind are considered
and continue the codes......
end
I assume you do the same thing for IY as well?
This seems similar (but not identical) to how the narrow band level set method has been implemented in ELMFIRE. Basically LIST_TAGGED tracks all cells that are within BANDTHICKNESS cells of the fire front. That way you calculate gradients and spend cpu time only in a narrow band adjacent to the fire front. The LIST_TAGGED list is generated by calling the subroutine TAG_BAND (in elmfire_level_set.f90) with inputs IXLOC and IYLOC (x and y locations to around which to tag a band), here's the subroutine:
! *****************************************************************************
SUBROUTINE TAG_BAND(NX, NY, IXLOC, IYLOC, T)
! *****************************************************************************
INTEGER, INTENT(IN) :: NX, NY, IXLOC, IYLOC
REAL, INTENT(IN) :: T
INTEGER :: IXTAGSTART, IXTAGSTOP, IYTAGSTART, IYTAGSTOP, IX, IY
IXTAGSTART = MAX(3, IXLOC - BANDTHICKNESS)
IXTAGSTOP = MIN(NX-2, IXLOC + BANDTHICKNESS)
IYTAGSTART = MAX(3, IYLOC - BANDTHICKNESS)
IYTAGSTOP = MIN(NY-2, IYLOC + BANDTHICKNESS)
DO IY = IYTAGSTART, IYTAGSTOP
DO IX = IXTAGSTART, IXTAGSTOP
IF (ISNONBURNABLE(IX,IY)) CYCLE
IF (.NOT. TAGGED(IX,IY) .AND. (.NOT. EVERTAGGED(IX,IY)) ) THEN
TAGGED (IX,IY) = .TRUE.
EVERTAGGED(IX,IY) = .TRUE.
CALL APPEND(LIST_TAGGED, IX, IY, .FALSE., T)
NUM_EVERTAGGED = NUM_EVERTAGGED + 1
EVERTAGGED_IX(NUM_EVERTAGGED) = IX
EVERTAGGED_IY(NUM_EVERTAGGED) = IY
ENDIF
ENDDO
ENDDO
! *****************************************************************************
END SUBROUTINE TAG_BAND
! *****************************************************************************
The 2D matrices TAGGED(IX,IY) and EVERTAGGED(IX,IY) are used to avoid adding a tagged cell to the tagged list more than once.
Yes, its for IY as well.
So basically I need to declare a new variable that is the same type as LIST_TAGGED?
I use CYCLE. It tremendously speed up the simulation. But if I use it, the building model is not executed.
This is my snippet
IF (.....) THEN
BUILDING_MODEL
ELSE
CYCLE
ENDIF
This input deck was posted at Issue #35 which has subsequently been closed, re-posting it here:
https://drive.google.com/file/d/1QZNJwqDyKahPp-am9lqUfx1LAvFSs_Qr/view?usp=sharing