# Urban Wind CFD with gloBFPr and OpenFOAM

    library(gloBFPr)
    library(sf)
    library(terra)
    library(ggplot2)

`gloBFPr` connects three steps in one R session: collecting urban spatial data, writing OpenFOAM configuration files, and visualising CFD results as maps.

# 1 Quick start

The package ships with a small Detroit-area dataset. No API keys are needed.

    data(globfp_example)
    buildings_list <- list(poly = globfp_example, binary = NULL, graduated = NULL)
    foam_case <- file.path(path.expand("~"), "openfoam_quickstart")

    foam_inputs <- prepare_openfoam_inputs(
      case_dir            = foam_case,
      buildings_list      = buildings_list,
      height_col          = "Height",
      include_buildings   = TRUE,
      include_fused_dsm   = FALSE,
      include_tree_canopy = FALSE,
      include_morphology  = FALSE,
      include_neighbors   = FALSE,
      include_greenspace  = FALSE,
      overwrite = TRUE, quiet = TRUE
    )

    case_files <- prepare_foam_case(
      case_dir = foam_inputs$case_dir,
      stl_file = foam_inputs$files$building_stl,
      domain   = foam_inputs$domain,
      inlet_velocity = c(5, 0, 0), z_ref = 10, z0 = 0.1,
      base_cell_size = 5, sim_hours = 0.5, overwrite = TRUE
    )

    run_openfoam_docker(case_dir = foam_inputs$case_dir,
                        image = "opencfd/openfoam-run:2506", wait = TRUE)

    maps <- read_foam_pedestrian_slice(case_dir = foam_inputs$case_dir,
                                       base_cell_size = 5, resolution = 2)

    plot_foam_map(maps, layer = "U_mag",
                  palette = "YlOrRd", legend_title = "Speed (m/s)")

# 2 Inputs

A case draws on two input sources: spatial layers from `prepare_openfoam_inputs()` and meteorological conditions from `get_era5_met()`. Set them up once, then pass the results to any simulation.

## 2.1 Spatial layers

`prepare_openfoam_inputs()` collects four categories of data and writes them into a structured case folder:

Layer | File | Used for
--- | --- | ---
Building footprints | `constant/triSurface/buildings.stl` | Solid geometry (`snappyHexMesh`)
Fused DSM | `rasters/fused_dsm.tif` | Terrain source — see Section 2.2
Building heights | `rasters/building_height.tif` | Terrain recovery
Canopy height (CHM) | `rasters/canopy_height.tif` | Canopy drag zone
Ground roughness (z0) | `rasters/ground_roughness_z0.tif` | Wall roughness

    foam_inputs <- prepare_openfoam_inputs(
      case_dir            = "~/openfoam_cases/detroit",
      bbox                = c(-83.065644, 42.333792, -83.045217, 42.346988),
      # target_crs is optional — the local UTM zone is chosen automatically
      include_buildings   = TRUE,
      include_fused_dsm   = TRUE,           # OpenTopography DEM + buildings + canopy
      include_tree_canopy = TRUE,           # Meta CHM → canopy drag
      canopy_source       = "metachm",
      min_tree_height     = 2,
      include_morphology  = TRUE,
      include_neighbors   = TRUE,
      include_greenspace  = TRUE,           # land cover → z0 raster
      mask_tree_cover     = TRUE,
      landcover_source    = "esri",         # "esa" (2020-2021) or "esri" (2017-2025)
      landcover_year      = 2022,
      opentopo_key        = Sys.getenv('OPENTOPO_API'),
      overwrite           = TRUE
    )

## 2.2 Terrain and canopy geometry

`prepare_openfoam_inputs()` writes terrain and canopy as rasters. `prepare_foam_geometry()` turns them into the STL geometry the mesher needs. This step is optional — skip it and you get a flat ground plane. If you pass terrain or canopy files to `prepare_foam_case()`, they must exist; explicitly supplied but unreadable inputs stop the case setup instead of being skipped.

    geo <- prepare_foam_geometry(
      case_dir        = foam_inputs$case_dir,
      fused_dsm       = foam_inputs$files$fused_dsm,
      building_height = foam_inputs$files$building_height_raster,
      canopy_height   = foam_inputs$files$canopy_height_raster,
      buildings       = readRDS(foam_inputs$files$buildings_rds),
      height_col      = "Height",
      domain          = foam_inputs$domain,
      base_cell_size  = 10
    )
    # geo$terrain_stl   → constant/triSurface/terrain.stl
    # geo$canopy_stl    → constant/triSurface/canopy.stl
    # geo$building_stl  → buildings.stl, re-based onto the terrain
    # geo$dem           → recovered bare-earth SpatRaster

Canopy becomes distributed drag, not solid blocks. A tree crown passes and drags air; it does not block it. Meshing crowns as solid over-blocks the flow. The canopy volume selects cells for a `leafAreaDensity` field, which drives `atmPlantCanopyUSource` and `atmPlantCanopyTurbSource` in `constant/fvOptions`

Drag is proportional to leaf area density, so cells outside the canopy are left at zero and feel nothing — which is why `selectionMode all` needs no cell set. Tune with `leaf_area_density` (default 0.4 m⁻¹; roughly 0.2 for sparse street trees, 0.8 for dense woodland) and `plant_cd` (default 0.2).

## 2.3 Ground roughness z0 from land-cover data

When `include_greenspace = TRUE`, land cover is reclassified to aerodynamic roughness length z0 and saved to `ground_roughness_z0.tif`. It sets `nutkRoughWallFunction` on the ground patch, so unresolved surface texture still slows the flow near the ground.

`"esa"` — ESA WorldCover (years 2020, 2021 only):

ESA class | Description | z0 (m)
--- | --- | ---
10 | Tree cover | `NA` — canopy drag zone
20 | Shrubland | 0.20
30 | Grassland | 0.05
40 | Cropland | 0.05
50 | Built-up | 0.03 — paved ground between resolved buildings
60 | Bare / sparse | 0.01
70 | Snow / ice | 0.001
80 | Permanent water | 0.0002
90 | Herbaceous wetland | 0.05
95 | Mangroves | 0.50
100 | Moss and lichen | 0.01

`"esri"` — Sentinel-2 10 m ESRI LULC Time Series (years 2017–2025). Use this when you need a year outside 2020–2021, or want consistent annual maps for before-and-after comparisons:

ESRI class | Name | z0 (m)
--- | --- | ---
2 | Trees | `NA` — canopy drag zone
11 | Rangeland — parks, lawns, pastures | 0.10
5 | Crops | 0.05
4 | Flooded vegetation | 0.05
7 | Built area | 0.03 — paved ground between resolved buildings
8 | Bare ground | 0.01
9 | Snow / ice | 0.001
1 | Water | 0.0002
10 | Clouds (no data) | 0.03 (fallback)

> Built-up z0 is 0.03 m, not the 0.5 m typical of unresolved urban terrain. Buildings here are explicit geometry, so the roughness applies only to the paved surface between them. Applying an urban-canopy z0 on top of resolved buildings would count the same drag twice. For the same reason tree cover is masked to `NA` — it is handled as canopy drag instead.

    z0_rast <- rast(foam_inputs$files$roughness_raster)
    z0_mean <- global(z0_rast, "mean", na.rm = TRUE)[[1]]
    cat(sprintf("Mean z0 = %.4f m\n", z0_mean))

## 2.4 ERA5 weather conditions

`get_era5_met()` fetches hourly ERA5 reanalysis from the Copernicus Climate Data Store. A wind case needs the 10-m wind components; `T_ref` is used as the uniform reference temperature that switches buoyancy off.

One-time setup: register at https://cds.climate.copernicus.eu, copy your personal access token, then store it:

    # Add to ~/.Renviron so it loads automatically every session
    CDS_API_KEY=your-token-here

    # Daytime — fetches 10-m wind components and 2-m temperature
    met <- get_era5_met(
      lon      = -83.05,
      lat      =  42.34,
      datetime = "2023-07-15 14:00",
      cds_key  = Sys.getenv("CDS_API_KEY")
    )
    # met$inlet_velocity  →  c(u10, v10, 0) m/s, ready for prepare_foam_case()
    # met$T_ref           →  2-m temperature in K

# 3 Wind simulation

## 3.1 Write the case

    case_files <- prepare_foam_case(
      case_dir            = foam_inputs$case_dir,
      # geo$building_stl when you ran Section 2.2, otherwise the flat-ground STL
      stl_file            = if (exists("geo") && !is.null(geo$building_stl)) {
                              geo$building_stl
                            } else {
                              foam_inputs$files$building_stl
                            },
      domain              = foam_inputs$domain,
      inlet_velocity      = met$inlet_velocity,   # ERA5 u10, v10 — any direction
      z_ref               = met$z_ref,            # 10 m
      T_ref               = met$T_ref,            # uniform: buoyancy inactive
      z0                  = z0_mean,
      terrain_stl         = if (exists("geo")) geo$terrain_stl else NULL,
      terrain_dem         = if (exists("geo")) geo$dem else NULL,
      canopy_stl          = if (exists("geo")) geo$canopy_stl else NULL,
      base_cell_size      = 15,
      building_refinement = 2L,
      sim_hours           = 1,
      overwrite           = TRUE
    )
    case_files$params$patch_roles
    #> xMin      xMax      yMin      yMax
    #> "inlet"   "outlet"  "lateral" "lateral"

### Any wind direction works

`blockMesh` emits four separately named lateral patches — `xMin`, `xMax`, `yMin`, `yMax` — and each is assigned an inlet, outlet or lateral role from the sign of `dot(flowDir, outward_normal)`. A south-westerly gives two inlets and two outlets:

    prepare_foam_case(..., inlet_velocity = c(3, 3, 0))
    #> Wind: 4.24 m/s at 10 m, dir (0.707 0.707)
    #> patches: xMin=inlet xMax=outlet yMin=inlet yMax=outlet

> Earlier versions pinned the inlet to the x-min face and expressed direction as a `flowDir` vector on it. For a northerly (`c(0, 5, 0)`) that injected momentum through a face whose normal is −x, so the inflow normal component was zero and effectively nothing entered — silently, with no error. If you have results from an earlier non-x-aligned legacy OpenFOAM run, re-run them.

## 3.2 Run via Docker

    run_openfoam_docker(
      case_dir = foam_inputs$case_dir,
      image    = "opencfd/openfoam-run:2506",
      wait     = TRUE
    )

## 3.3 Visualise

    maps <- read_foam_pedestrian_slice(
      case_dir       = foam_inputs$case_dir,
      T_ref          = case_files$params$T_ref,
      base_cell_size = case_files$params$base_cell_size,
      resolution     = 2
    )

    # Wind speed — buildings overlaid automatically
    plot_foam_map(maps, layer = "U_mag",
                  palette = "YlOrRd", legend_title = "Speed (m/s)",
                  title = "Pedestrian wind speed at 1.5 m AGL")

    # Speed with flow vectors
    p <- plot_foam_map(maps, layer = "U_mag",
                       palette = "Blues", reverse = TRUE,
                       legend_title = "Speed (m/s)")
    add_flow_vectors(p, maps, spacing = 20, colour = "grey20",
                     alpha = 0.7, linewidth = 0.25, arrow_size = 0.07)

    # Wind speed ratio (> 1.3 = potentially uncomfortable)
    plot_foam_map(maps, layer = "U_mag", palette = "RdYlGn", reverse = TRUE,
                  legend_title = "U / U_ref",
                  max_u_ref = sqrt(sum(met$inlet_velocity^2)))

> All `plot_foam_map()` calls return a plain `ggplot` object. Export with `ggsave("wind_map.png", width = 10, height = 8, dpi = 300)`.

The temperature field is uniform, so `T_cool` is ≈ 0 and only the velocity layers carry information.

> The map above is an instantaneous field, not a mean. The case is transient and `read_foam_pedestrian_slice()` takes the latest write, so in the wakes behind large buildings it shows whichever phase of the vortex shedding the run stopped on. Pedestrian comfort is assessed on a mean.

> The case therefore also runs a `fieldAverage` function object and writes a time-averaged slice of its own, starting one flow-through in (`endTime / 3`) so the startup transient is excluded:

>     postProcessing/pedestrianSliceMean/<time>/z1p5m_UMean.raw

> Same raw format as the instantaneous slice — x, y, z then the three velocity components. `read_foam_pedestrian_slice()` does not read it yet; for now parse it directly, or compare it against the instantaneous map to see how much the wakes move.

> `read_foam_pedestrian_slice()` reads the 1.5 m slice that the case writes automatically. To sample other heights or other fields — `p_rgh`, `k`, `nut` — use `sample_foam_slice()` instead; it runs `postProcess` in the container after the fact, so you can re-sample without re-running the simulation:

>     slice <- sample_foam_slice(case_dir = foam_inputs$case_dir,
>                                fields = c("U", "p_rgh"), z = 10,
>                                image = "opencfd/openfoam-run:2506", resolution = 5)

> Note the field is `p_rgh` (pressure minus hydrostatic), not `p`: the buoyant solver does not carry a separate kinematic pressure. Static pressure is `p_rgh + ρgh` if you need it.