Design: new geometry and material input format for SEB

[mvreeuwijk]

New Geometry and Material Format for uDALES SEB

Problem

Current SEB inputs are spread across ~10 text files with rigid positional columns. No surface grouping (faces → surfaces → materials), no standard format readable by other tools, and all materials are forced to share the same number of thermal layers (nfaclyrs is global).

High-Level Architecture: Three Files

1. Urban model — OBJ + MTL

OBJ natively supports everything we need:

• Vertices and face connectivity (polygons, not just triangles)
• Groups (g) — faces grouped into named surfaces (e.g. g roof_bldg1, g wall_north_bldg1)
• Materials (usemtl) — material assignment per face/group
• MTL companion file — material visual properties (albedo as Kd, roughness, etc.)
• Text-based, ~100-200 MB for a few million faces — manageable
• Opens directly in Blender, SketchUp, ParaView, MeshLab — zero conversion
• Python: trimesh reads and writes OBJ with groups and materials

What OBJ/MTL cannot carry: thermal layer stacks, conductivities, emissivities, view factors. These are uDALES-specific.

2. Material/thermal properties — TOML extension

A TOML file extending MTL with uDALES-specific properties per material name. Each material defines its own layer stack as parallel arrays:

[concrete]
albedo = 0.5
emissivity = 0.85
z0m = 0.05
z0h = 0.00035
lGR = false
d      = [0.12, 0.12, 0.12]
C      = [2500000, 2500000, 2500000]
lambda = [1.0, 1.0, 1.0]

[green_roof]
albedo = 0.25
emissivity = 0.95
z0m = 0.05
z0h = 0.00035
d      = [0.05, 0.10, 0.15, 0.10, 0.20]
C      = [2000000, 3000000, 5000000, 2500000, 2500000]
lambda = [0.5, 1.2, 2.0, 1.0, 1.0]

Material names match the MTL material names used in the OBJ. Easy to read and edit, arrays are natural in TOML. Python reads via tomllib (stdlib since 3.11).

3. Solver setup — NetCDF (udales.nc.inp.<expnr>)

Everything the solver needs at startup: facet geometry, surface properties (z0, z0h, emissivity, albedo, lGR), IBM staggered-grid mapping (facet sections, solid points, boundary points), and initial profiles (currently prof.inp). Needed for all IBM simulations, with or without SEB. Fortran never parses OBJ — it reads this file only.

Naming follows existing Fortran arrays (facnorm, faca, facz0, facz0h, facem, faclGR, etc.).

4. SEB input — NetCDF (seb.nc.inp.<expnr>)

SEB-specific data, only needed when lEB=.true.: thermal layer properties, view factors, sky view factors, net shortwave, initial facet temperatures. References the same facet ordering as udales.nc.inp.

Solver NetCDF schemas

udales.nc.inp.<expnr> — solver setup (geometry + surface properties + profiles)

dimensions:
    nfcts                       # number of facets
    ndim = 3
    ktot                        # number of vertical grid levels
    nsolpts_u, nsolpts_v, nsolpts_w, nsolpts_c
    nbndpts_u, nbndpts_v, nbndpts_w, nbndpts_c
    nfctsecs_u, nfctsecs_v, nfctsecs_w, nfctsecs_c

variables:
    # --- Facet properties (naming follows Fortran) ---
    double facnorm(nfcts, ndim)             # outward normal
    double faca(nfcts)                      # facet area [m^2]
    int    facets(nfcts)                    # wall/material ID (identifies which material this facet has)
    double facz0(nfcts)                     # momentum roughness [m] (needed for wall functions, even without SEB)
    double facz0h(nfcts)                    # heat roughness [m]

    # --- Solid points per staggered grid ---
    int    solid_u(nsolpts_u, ndim)         # (i,j,k) on u-grid
    int    solid_v(nsolpts_v, ndim)
    int    solid_w(nsolpts_w, ndim)
    int    solid_c(nsolpts_c, ndim)

    # --- Boundary points per staggered grid ---
    int    bndpts_u(nbndpts_u, ndim)        # (i,j,k) on u-grid
    int    bndpts_v(nbndpts_v, ndim)
    int    bndpts_w(nbndpts_w, ndim)
    int    bndpts_c(nbndpts_c, ndim)

    # --- Facet sections per staggered grid ---
    int    secfacids_u(nfctsecs_u)          # facet ID for each section
    int    secbndptids_u(nfctsecs_u)        # boundary point ID
    double secareas_u(nfctsecs_u)           # section area [m^2]
    double bnddst_u(nfctsecs_u)            # distance to boundary point [m]
    double recpts_u(nfctsecs_u, ndim)      # reconstruction point coordinates [m]
    int    recids_u_u(nfctsecs_u, ndim)    # cell indices of rec point on u-grid
    int    recids_u_v(nfctsecs_u, ndim)    # cell indices of rec point on v-grid
    int    recids_u_w(nfctsecs_u, ndim)    # cell indices of rec point on w-grid
    int    recids_u_c(nfctsecs_u, ndim)    # cell indices of rec point on c-grid
    # ... repeated for _v, _w, _c

    # --- Initial profiles (currently prof.inp) ---
    double zf(ktot)                         # full-level heights [m]
    double thl(ktot)                        # initial liquid water potential temperature [K]
    double qt(ktot)                         # initial total water specific humidity [kg/kg]
    double u(ktot)                          # initial u-velocity [m/s]
    double v(ktot)                          # initial v-velocity [m/s]
    double tke(ktot)                        # initial TKE [m^2/s^2]

global attributes:
    :Conventions = "uDALES-1.0"
    :source_obj = "urban.obj"
    :source_materials = "materials.toml"

Note: Time-dependent forcings (lscale.inp, time-dependent shortwave) are out of scope for this file — they remain separate (potentially moving to NetCDF in future).

seb.nc.inp.<expnr> — SEB properties + radiation

dimensions:
    nfcts                       # must match udales.nc.inp
    ntotal_layers               # sum of nfaclyrs across all facets
    ntotal_nodes                # sum of (nfaclyrs+1) across all facets
    nnz                         # non-zero view factor entries

variables:
    # --- Per-facet layer counts (ragged) ---
    int    nfaclyrs(nfcts)                   # layer count per facet
    int    faclyr_offset(nfcts + 1)          # index into layer arrays
    int    facnode_offset(nfcts + 1)         # index into node arrays

    # --- Layer properties (ragged 1D) ---
    double facd(ntotal_layers)               # layer thickness [m]
    double faccp(ntotal_layers)              # volumetric heat capacity [J/(K m^3)]
    double faclam(ntotal_nodes)              # thermal conductivity at nodes [W/(m K)]

    # --- Initial temperature (ragged 1D) ---
    double facT(ntotal_nodes)                # temperature at nodes [K]

    # --- Radiation ---
    double svf(nfcts)                        # sky view factor [-]
    int    vf_row(nnz)                       # sparse VF row (0-based)
    int    vf_col(nnz)                       # sparse VF col (0-based)
    double vf_val(nnz)                       # sparse VF value

global attributes:
    :Conventions = "uDALES-seb-1.0"

Facet n has:

• nfaclyrs(n) layers, nfaclyrs(n)+1 nodes
• Layer data at facd(faclyr_offset(n) : faclyr_offset(n+1)-1)
• Node data at facT(facnode_offset(n) : facnode_offset(n+1)-1)

Data Model

User-facing files (OBJ + TOML):

• A face is a triangle (or polygon) in the OBJ mesh.
• Faces are grouped into surfaces via OBJ groups (g wall_north_bldg1).
• Each surface has a material assigned via OBJ usemtl.
• Each material has visual properties in MTL and thermal/physical properties in TOML (albedo, emissivity, z0m, z0h, lGR, and layer stacks of d, C, lambda — variable count per material).

Solver files (NetCDF):

• udales.nc.inp contains per-facet properties (normal, area, z0, emissivity, albedo), the IBM staggered-grid mapping (solid points, boundary points, facet sections for u/v/w/c grids), and initial profiles (currently prof.inp). Each facet section maps a portion of a facet to one boundary point on one staggered grid.
• seb.nc.inp contains thermal layer properties (ragged per facet), radiation (svf, netsw, view factors), and initial temperatures. Only needed when lEB=.true..
• Time-dependent shortwave (when ltimedepsw=.true.) lives in a separate NetCDF file — not in seb.nc.inp.

Key Requirement: Variable Layers Per Material

Each material defines its own layer count. In the solver:

• Replace global nfaclyrs with per-facet nfaclyrs(n) (read from seb.nc.inp) in the EB loop
• Matrix construction loop: do j = 1, nfaclyrs(n) instead of do j = 1, nfaclyrs
• Matrix inversion: gaussji(HM, IDM, nfaclyrs(n)+1) — already takes size as argument
• Remove redundant matinv4 fast-path
• Allocate work matrices to max(nfaclyrs)+1, reused per face
• Layer properties via ragged offset: facd(faclyr_offset(n) + j) instead of facd(n, j)
• Temperature nodes via offset: facT(facnode_offset(n) + j)

Variable Layers: Ragged Storage

Each facet has nfaclyrs(n) thermal layers → nfaclyrs(n) + 1 temperature nodes. Expected range: typically ~10, up to ~20. With millions of faces on a parallel code, padding to max wastes significant memory across ranks. Use ragged 1D arrays with offsets — the same pattern as face-node connectivity.

Layout

nfaclyrs(nfcts)               # layer count per facet
faclyr_offset(nfcts + 1)      # into layer arrays (facd, faccp): cumulative sum of nfaclyrs
facnode_offset(nfcts + 1)     # into node arrays (facT, faclam, facTdash): cumulative sum of nfaclyrs+1

facd(ntotal_layers)            # concatenated thicknesses
faccp(ntotal_layers)           # concatenated heat capacities
faclam(ntotal_nodes)           # concatenated conductivities (nfaclyrs+1 per facet)
facT(ntotal_nodes)             # concatenated temperatures (nfaclyrs+1 per facet)
facTdash(ntotal_nodes)         # concatenated temperature gradients

For facet n: layers at facd(faclyr_offset(n) : faclyr_offset(n+1)-1), nodes at facT(facnode_offset(n) : facnode_offset(n+1)-1).

Where this applies

• Solver NetCDF input (seb.nc.inp): ragged 1D arrays + offset arrays
• Solver Fortran arrays: same layout. EB loop indexes into flat arrays via offsets.
• Solver NetCDF output (facT.nc): same ragged format
• EB loop: work matrices allocated to max(nfaclyrs)+1, reused per facet. gaussji takes size as argument.
• Temperature initialisation: linear interpolation between outer and inner T, writing into the correct slice

Remaining Questions

1. Material extension format — resolved: TOML.

2. Surface assignment workflow — resolved: three-level approach:

   • split_surfaces — automatic: identify all walls, roofs, floors, ground per building (extends existing split_buildings)
   • merge_surfaces — combine multiple surfaces into one (e.g. merge all east-facing walls of a building)
   • Blender — for advanced editing: import OBJ, manually assign/rename groups, re-export. OBJ groups round-trip cleanly.
   • Use case example: select faces for a green wall → assign to a new surface with a different material. All remaining facets become surface 2 with the default material.

3. Migration — resolved: no dual code path in Fortran. The solver only reads NetCDF. A conversion tool (convert_legacy_to_nc.py) detects legacy text files (facets.inp, factypes.inp, etc.) and automatically produces the NetCDF. Old cases just need one conversion run. Small user base means we can make this clean break.

Implementation Phases

1. Define the formats — OBJ group/material conventions, TOML schema, NetCDF schema
2. Python writer — new preprocessing section: OBJ + TOML + solver NetCDF
3. Conversion tool — convert_legacy_to_nc.py: detect legacy text files → produce NetCDF
4. Fortran reader — replace text-file reader in initfac.f90 with NetCDF-only reader
5. Per-face variable layers — face_nlayers array in modEB.f90, remove matinv4
6. Postprocessing — teach udbase.py to read udales.nc.inp

Critical Files

• src/initfac.f90 — Fortran facet reader
• src/modEB.f90 — energy balance solver (EB loop, matrix construction)
• src/modstartup.f90 — namelist parameters
• tools/python/udprep/udprep_ibm.py — current facet file writer
• tools/python/udprep/udprep_seb.py — SEB preprocessing
• tools/python/udprep/udprep_radiation.py — radiation preprocessing
• tools/python/udgeom/udgeom.py — geometry class (trimesh-based)
• tools/python/udbase.py — postprocessing