Skip to content

param_to_var on PDESystem parameters doesn't add initial condition BCs #200

@ctessum

Description

@ctessum

Description

When param_to_var promotes a PDESystem parameter to a dependent variable during coupling, no initial condition (t=0 BC) is added for the new variable. This causes discretize to fail with:

Please provide a default (`u0`), initialization equation, or guess
for the following variables:
(E_wind(t))[2, 2], (R_0(t))[2, 2], (u_y(t))[2, 2], (B_wind(t))[2, 2], (C_wind(t))[2, 2]

Root cause

The ODE coupling path doesn't have this problem because icbc(di, statevars) generates ICs for all state variables when promoting via sys + domain. But in the PDE coupling path:

  1. couple2 calls param_to_var(pdesys, :R_0) to promote a parameter to a DV
  2. The connector equation wires it: pdesys.R_0 ~ other_sys.R0
  3. merge_pdesystems collects DVs and BCs from all subsystems
  4. The new DV R_0 has no t=0 BC because it was a parameter (not a DV) in the original PDESystem

Affected variables

Any PDESystem parameter promoted via param_to_var during coupling. In the wildfire model, this includes: R_0, C_wind, B_wind, E_wind, u_x, u_y, dzdx, dzdy, φs_coeff, β_ratio, and others.

Suggested fix

Either:

  1. In param_to_var(sys::PDESystem, ...): When promoting a parameter, also add a dv(t0, spatial_ivs...) ~ default BC to the PDESystem's boundary conditions.

  2. In merge_pdesystems: After merging, detect DVs that lack a t=0 BC and add dv(t0, spatial_ivs...) ~ 0.0 for algebraic variables.

Workaround

Add missing ICs manually after convert(PDESystem, cs):

let t0 = DomainSets.infimum(pde.domain[1].domain)
    has_ic = Set{Symbol}()
    for bc in pde.bcs, dv in pde.dvs
        dvname = Symbol(Symbolics.tosymbol(dv, escape=false))
        if occursin(string(dvname, "(", t0), string(bc.lhs))
            push!(has_ic, dvname)
        end
    end
    new_bcs = copy(pde.bcs)
    for dv in pde.dvs
        dvname = Symbol(Symbolics.tosymbol(dv, escape=false))
        dvname  has_ic && continue
        args = Symbolics.arguments(Symbolics.unwrap(dv))
        op = Symbolics.operation(Symbolics.unwrap(dv))
        push!(new_bcs, op(t0, args[2:end]...) ~ 0.0)
    end
    if length(new_bcs) > length(pde.bcs)
        pde = PDESystem(equations(pde), new_bcs, pde.domain, pde.ivs, pde.dvs, pde.ps;
            name = nameof(pde), initial_conditions = pde.initial_conditions, checks = false)
    end
end

Related

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions