fix build_parabolic_dz_array()

[zacharycope0]

Feature Description The LANL team gave me some MATLAB code for building a dynamic dz-array depending on the height of the QU-grid (wind grid) that we use. We want the grid to be greater(3x topo or (tallest fuel height + 100m). I tried to write the function as build_parabolic_dz_array() in print_inp_files.py and I think I created an infinite loop.

Getting Started

1. Fix build_parabolic_dz_array() (see screenshots). Feel free to start from scratch if it's easier.
2. Change lines printing dz array in print_QUIC_fire_inp(dom). Should look something like this. dz_array = build_parabolic_dz_array(nz=22, Lz=350, n_surf=4, dz_surf=2.5) input_file.write("{} ! Number of vertical layers of fire grid cells (integer)\n".format()) input_file.write("1 ! x - fire grid ratio = (QUIC-URB cell size)/(fire cell size), integer, can only be >= 1\n") input_file.write("1 ! y - fire grid ratio = (QUIC-URB cell size)/(fire cell size), integer, can only be >= 1\n") input_file.write("0 ! Vertical stretching flag: 0 = uniform dz, 1 = custom\n") for dz in dz_array: input_file.write("{}\n".format(dz))
3. Change the build_parabolic_dz_array() call so that: nz = 22 Lz = greater(3x topo or (tallest fuel height + 100m) see issue Implement Topography n_surf = 5 dz_surf = 1

Screenshots LANL description of function: MatLab Code: dzmax_high=zmax-dzmin*nucells;

    dzmax_low=0;

    zmax_temp=ncells*dzmin;

    dz=dzmin*ones(1,ncells);

    while abs(1-zmax_temp/zmax) > 0.001;

        dzmax=0.5*(dzmax_low+dzmax_high);

        C1=(dzmax-dzmin)/((ncells-nucells)^2);

        C2=-2*C1*nucells;

        C3=dzmin+C1*nucells^2;

        dz=zeros(1,ncells);

        dz(1:nucells)=dzmin;

        for k=nucells+1:ncells;

            dz(k)=(C1*k^2)+(C2*k)+C3;

        end;

        zmax_temp=sum(dz);

        if zmax_temp>zmax;

            dzmax_high=dzmax;

        elseif zmax_temp<zmax;

            dzmax_low=dzmax;

        else

            break;

        end;

    end;

    diff=zmax-sum(dz);

    dz(end)=dz(end)+diff;

    handles.domain.dz_array=dz;

    zmax=sum(handles.domain.dz_array);

[cbonesteel]

Couple questions.

Firstly, for the default function values, is Lz = greater(3x topo or (tallest fuel height + 100m) an option of which to implement?

Secondly, is dz_array = build_parabolic_dz_array(nz=22, Lz=350, n_surf=4, dz_surf=2.5) correct for the input file? It seems to contradict the idea of the default values or the idea of using the domain for some of these like nz (unless this nz is different).

[zacharycope0]

First question: No I don't want it to be a user option. I want the code to automatically build a QU_grid (wind grid) that will handle the domain that the user is specifying. If they use topo the wind grid needs to be 3x the change in elevation. However, if there is very little topo (like a relief of 10m across the domain) I'd rather use fuel height & 100 m.

Secondly: You are correct. The values that I left in there are just to make sure the function is working as expected (see dx[m] in screen shot). Once the function is fixed we don't want the values hard coded. We will want Lz = greater(3x topo or (tallest fuel height + 100m), n_surf=5, dz_surf=1, and probably nz=22. We might want to make nz bigger (up to 30 vert cells) if there is a lot of topo, but we can worry about that later.

[cbonesteel]

I implemented the matlab script and put the dz_array = build_parabolic_dz_array(nz=22, Lz=350, n_surf=4, dz_surf=2.5) function call in to test and the data is off. Instead of cell 5 being 2.639853 it is 2.5. Everything else after this is off as well. I can't seem to find what the issue is and attempts at debugging it against a matlab have been unsuccessful. I might be missing something obvious but another pair of eyes would be appreciated.

Once this is fixed I will implement the different options for Lz but for now it has one default value.

[cbonesteel]

I implemented the matlab script and put the dz_array = build_parabolic_dz_array(nz=22, Lz=350, n_surf=4, dz_surf=2.5) function call in to test and the data is off. Instead of cell 5 being 2.639853 it is 2.5. Everything else after this is off as well. I can't seem to find what the issue is and attempts at debugging it against a matlab have been unsuccessful. I might be missing something obvious but another pair of eyes would be appreciated.

This was discussed and decided that there would be some error compared to the Matlab code and the current implementation is fine.

The solution for Lz is below. It finds the terrarion relief if topography is being used and compares the relief * 3 to 100. If it is larger, the fuel_height + relief * 3 formula is used and in all other cases it is just fuel_height + 100. I assume I found the fuel_height correctly but a double check would be appreciated.

https://github.com/QUIC-Fire-TT/ttrs_quicfire/blob/1b0e7d5561410245bb33a83848e08deb242f5590/ttrs_quicfire/print_inp_files.py#L302-L310

[zacharycope0]

I added my reply to the topo issue: https://github.com/QUIC-Fire-TT/ttrs_quicfire/issues/2#issuecomment-1144071388

[zacharycope0]

I approved the pull request after changing how the fuel height was calculated: https://github.com/QUIC-Fire-TT/ttrs_quicfire/blob/dc19fbcaa8cfe8606155c0318385fb9b9ebfaceb/ttrs_quicfire/print_inp_files.py#L302-L305 Your solution worked but was harder to follow. I also wanted the fuel height to be the top of the fuel and with PY indexing it was 1m off. Great work!

[zacharycope0]

Feel free to close after reading this.

[cbonesteel]

Your solution is definitely more elegant! I will go ahead and setup a new release so pip install works with these latest changes. Should be a minor release so 1.1.0 would be the new version.