In this paper we outline a strategy for rapid construction of prototype computational schemes for complex physical problems. We work with an algorithmic framework based on Cartesian grids and semi-implicit iterative schemes: this common and well studied framework is robust and simple to use, while it allows to focus on application-specific modeling challenges of prototype models for complicated physical applications, for which body fitting and unstructured grids are frequently replaced with immersed boundary and fictitious domain techniques [19,13]. These features have been successfully used for multiphysics applications and frameworks, see e.g., in [3,1].
Our focus in this paper is on systems which involve some form of mixed-size domain and couplings for models involving some discretizations of PDEs. Recently, there is abundance of theoretical and applications oriented work on mixed dimensional setting for a variety of applications: blood flow, hydro-mechanical coupling in fracture networks, and more. The theoretical basis: detailed delicate functional analytic setting of the PDEs and of the underlying finite element approximations, is delicate since it involves distributed line sources or hybrid coupled continua-network models; see, e.g., [29,12,16] following earlier modeling and applications oriented work including [7,9,11].
Our approach here is to start from images which define the geometry of (a portion of) the computational domain Ω ⊂ R d which we project on a Cartesian grid; we represent the lower dimensional domains by unions of the grid/voxel cells embedded within the original Ω; see Figure 1 for an illustration. Next we take advantage of the concept of overlapping continua [28,6], asymptotically correct in domains separated by low-conducting interface, or of the homogenized models [15,10] to handle the mixed dimensional coupling, while working with the original continuum models. For discretization, we use the well known practical computational framework of lowest order mixed finite element methods implemented as cell centered finite differences (CCFD), and we develop some model approximations and enhancements suitable for the selected phenomena.
This strategy allows to build prototype models, propose their extensions and enhancements, and investigate their robustness, flexibility, and sensitivity with case studies which guide how to extend, modify, and interpret the models, while paying attention to numerical discretization (h, τ ), and implementation (solver). This paper illustrates the simplicity and potential of this methodology which we recommend as a stepping stone towards refinement strategies beyond the applications we discuss.
We select three disparate applications dubbed H, RS, and T which involve coupled variables defined on model-and variable-specific domains which may overlap, with mixed dimensional setting; see Figure 1. The computational models solved on these domains are discretizations of PDEs which communicate on the intersection of the domains, or are "aware" of the variables defined on these other domains. The models H, RS and T are on modeling hypothermia problems in tissue, water uptake in large root-soil systems, and mixed-use traffic flow on a campus domain. From solver point of view, case H is semilinear, and case RS and T require care in the treatment of advective terms. The outline of this paper is as follows. In Section 2 we recall the basic algorithm for a scalar parabolic-hyperbolic PDE implemented as CCFD. In Section 3 we introduce the H (hypothermia) model, an extension of the linear bioheat flow model [20,10,15] for which we exploit fictitious domain approach as well as introduce the nonlinear feature of vasoconstriction. In Section 4 we employ an extension of d = 1 Richards-Darcy models for RS (root-soil) systems from [4,27] to mixed dimensional setting in d ≥ 2. In Section 5 we develop a T model for mixed-species traffic flow on a campus network involving several species whose trajectories may or may not be confined to pathways, and which may or may not be aware of other traffic participants.
While due to the scope we cannot cover all the details, the models H, RS, and T can help in further construction of more refined models, also for other complex applications. Such refined models can take advantage of sensitivity studies carried out with the prototype models.
In the paper we apply the following common notation. First, χ S denotes the characteristic function of a set S, |S| its measure, and S its closure. For a function f on S, ⟨f ⟩ S denotes its average on S. We work in spatial domains which are open bounded sets Ω ⊂ R d , d = 1, 2, 3, possibly partitioned into some material subdomains, open sets Ω m . The case d ≤ 2 is natural for Section 5, and the problems in Sections 3 and 4 work with d ≤ 3. For the boundary ∂Ω or that of any of the subdomains ∂Ω m , η is a unit normal vector pointing outward, and we assume that these boundaries are sufficiently smooth. We will consider Dirichlet conditions on ∂Ω D and flux conditions on ∂Ω N .
We will denote by 0 ≤ t ≤ T the time variable, with T the final time of the simulation. In time-discrete models, t n is the time step, and t 0 = 0 denotes the initial time. With uniform time stepping we have t n = nτ .
We consider first a generic scalar parabolic quasi-linear homogeneous PDE
which is supplemented by some boundary and initial conditions. Assume D is a symmetric uniformly positive definite tensor, c is uniformly nonnegative, f is sufficiently smooth, and C(u) = c b (x, u, t)u is the Helmholtz term. (See (H) model below). For nonlinear implicit parabolic problem, D = D(u), and c∂ t u derives from some ∂ t C(u) (See (RS) model below). When c = 0, C = 0, (1) is a steady flow problem, an elliptic PDE (See (T) model below).
For spatial discretization, we assume that Ω = Ω h = (ij)∈T h ω ij on some underlying background Cartesian grid of rectangular cells ω ij of center x ij whose edges E h align with ∂Ω and with the material interfaces, with max ij |ω ij | = h 2 . The cells are identified with voxels (image pixels); some ω ij ̸ ∈ Ω are "key-outs" outside the union T h of indices as in Figure 1. The unknown u is approximated by piecewise constants u ij ≈ u(x ij ); we collect u h = (u ij ) (ij)∈T h . The flux components q h defined over the edges E h are from the mixed finite element RT [0] space. This discretization is well known [21,18] for its locally conservative properties and easiness for modeling nonlinear multi-physics applications. Since Dirichlet conditions are satisfied weakly [21], fictitious domain strategies [13] are their natural extension.
Applying implicit-explicit time stepping, we seek u n h which solves
Here M h is the mass matrix, and A h is the stiffness matrix incorporating the boundary conditions; C applies pointwise to each degree of freedom of u n h . Further, F h is the numerical flux defined based on f , and ∇ h • is the discrete counterpart of ∇• on the grid T h . For scalar problems we use Godunov flux; when u is a vector, we require a Riemann solver [8,17]. For nontrivial F h , we apply operator splitting: we solve the advection portion c∂ t u + ∇ • f (u) = 0 explicitly, followed by an implicit diffusion-reaction step solved directly or by iteration, time-lagging the nonlinearity and/or the coupling. After u n h is found, we retrieve the fluxes q n h . In what follows we drop the reference to h, n.
We consider next a system based on (1)
The model (3) when Ω 1 = Ω 2 is known as the Barenblatt model for multiscale flow in porous media [6,28] postulated for the flow dynamics in composite multiscale media with vastly different storage c j , D j . In this paper Ω 1 is a continuum, and Ω 2 is "almost" one-dimensional, Ω 2 ⊂ Ω 1 with 0 < |Ω 2 | << |Ω 1 | in the same R d measure, and c = c(x) = Cχ Ω1∩Ω2 (x), C ≥ 0, only active on Ω 1 ∩ Ω 2 . Our scheme for (3) extends directly those in Section 2.1, and the coupling can be resolved by iteration or directly. This alternative to mixed dimensional setting allows an easy proof-of-concept complex nonlinear dynamics. We consider three models which we call in shorthand by H, RS, and T. In H model of hypothermia we have constant coefficients except possibly nonlinear c b to model vasoconstriction. For the root-soil problem RS, c 2 = 0, and the model for u 1 is the nonlinear Richards equation, where D 1 = D 1 (u 1 ) and c 1 = c 1 (u 1 ) feature degenerate behavior. We also consider a model T of overlapping continua involving mixed traffic participants utilizing different trajectories. In that model f (x, u) = v(x, u) f (u), and v(x, u) is computed solving an auxiliary elliptic problem similar to (1) with c = 0, C = 0.
We aim to build a qualitatively realistic model for the temperature in the complex tissue domain Ω such as in Figure 1 (a) connected to the rest of the body Ω body (not shown). When exposed to low external temperatures, the body thermoregulation system attempts various strategies to prevent hypothermia and tissue damage. While these mechanisms are complex and not completely understood by physiologists, our prototype model based on overlapping continua and fictitious domain concepts could eventually aid, e.g., in developing patientspecific hypothermia therapy or cryo-preservation strategies which require rapid turnaround time from an image to simulation.
The body is composed of Ω vessel , Ω capillaries , Ω tissue . Heat transport is by convection in Ω vessel , Ω capillaries and by conduction in Ω tissue . Metabolic sources and products and energy are exchanged between Ω capillaries and Ω tissue ; the latter process, called perfusion, makes heat conduction in Ω capillaries and Ω tissue similar to flow in porous medium; see, e.g., recent developments in [24].
Literature offers various models for thermal perfusion that simplify or enhance the overlapping continua model ( 3) on Ω = Ω tissue , with the blood temperature approximated as θ * = const, the arterial temperature. In particular, the so-called Pennes model from [20] postulates
Here c and k are the volumetric heat capacity and thermal conductivity, and
, where c bh and v b are blood volumetric heat capacity and perfusion rate. Here and below we ignore heat sources. In [10], c b (θ-θ * ) is derived by homogenization techniques in Ω capillaries , Ω tissue made of ϵ-size unit cells Y with a Robin boundary condition imposed on the interface Γ = Y ∩∂Ω capillaries ∩ ∂Ω tissue , and v b is shown to be proportional to |Γ |. In [7], Ω vessel is treated as a 1D domain, with convection in Ω vessel coupled to that in Ω capillaries ∪Ω tissue ⊂ R 3 , and several multi-equation models use a similar approach, e.g. [14,26].
Now we aim to apply (4) in a realistic setting involving Ω vessel and a complicated external boundary ∂Ω. For the former, we include a penalty term c vessel (θ -θ * ) on Ω vessel to enforce θ| Ω vessel ≈ θ * ; for the latter, we proceed similarly on some Ω \ Ω with Ω of simpler shape than Ω. These treatments resemble the Immersed Boundary [19] and fictitious domain [13] approaches, respectively. We define now
Similarly, we set θ * = θ * (x) equal θ body , θ vessel , θ D , respectively, and postulate an extension of (4) to Ω as follows
The coefficients c, k can be extended to Ω in any convenient way as long as ( 5) is well-posed. The model ( 5) is next approximated by the mixed finite elements/CCFD setting as described in Section 2.1.
We illustrate the effect of the blood perfusion term c b (θ -θ * ) and fictitious domain term c D (θ -θ * ) in a hypothermia example based on geometry in Figure 1 (a) aiming for a physically motivated scenario to simulate the onset of frostbite in the human hand exposed to cold air over t ≤ T = 1200 s. We extend Ω = Ω hand ⊂ Ω by Ω air so that Ω is a large enough rectangular domain. The tissue in Ω hand has heterogeneous properties including large blood vessels. We identify ∂Ω wrist = ∂Ω ∩{x :
We start with a d = 1 slice of the 2d domain where the effect of c b , c D is easily quantified. Convergence studies (not shown) show expected order O(τ + h 2 ). We consider next the 2d example, where the hand is protected by mitten material in Ω mitten , case (M ), or not, case (N ). Other notation is adapted easily. The model exhibits qualitatively intuitive behavior, as shown in Figure 3. Next we post-process the results to analyze the extent of hypothermia and possible frostbite, i.e. θ(x, t) < 0 • C. Frostbite is avoided in case (M ) at least until t = 1200, but occurs in case (N ), affecting about 11% of cells in Ω hand .
Next we address the quality of model adaptation. The term C(θ) simulates perfusion in Ω hand \ Ω vessels and acts as a penalty term in Ω vessels and Ω \ Ω.
While c b has physiological meaning, c vessel and c D are chosen somewhat ad-hoc to ensure, respectively, θ| Ω vessel ≈ θ * | Ω vessel and θ| ∂Ωnonwrist ≈ θ * | Ωair . We test the quality of the adapted model by checking if this first condition holds. We record θ max = max ij θ n ij at t n = T . When c vessel = 10 5 , case (N ) has θ max = 35.7. However, when c vessel = 1, case (N ) has θ max = 9.8. Case (M ) is similar.
With model sensitivity to c b , c D , c vessel reasonably understood, we extend (4) to account for defense against hypothermia with vasoconstriction, wherein body temperature information captured by thermoreceptor neurons induces arterial smooth muscle constriction, reducing blood flow in extremeties and retaining heat near the core body. Such action can prevent core hypothermia, but may cause permanent morphological changes, e.g. frostbite, in extremeties.
To model vasoconstriction, we recall c b depends on Γ . Thus a natural step is to include in (4) some nonlinear dependence of c b on θ| Ω . We hypothesize that c b decreases when a decrease in θ is detected through sensory data S(θ), e.g. S(θ) = ⟨θ⟩ Ω . Our simulations suggest that when c b decreases, so will S(θ), affecting the venous blood temperature ≈ S(θ) returning to the body. Therefore, without increased metabolism, the body core temperature θ body may decrease, further decreasing S(θ), c b , and θ ∂Ωwrist . We postulate the following simple model
with µ ≥ 0 representing metabolism and λ a Lagrange multiplier ensuring θ body ≤ 37, mimicking thermoregulation. We solve (6) coupled to (4) semiimplicitly, requiring some c B ≥ 0, µ ≥ 0, and some model for c b = c b (θ body ). 37 . I assume u B is θ body ?. For dramatic effects, we consider µ = 0. We study the value θ body (t) and retrieve the position x * ∈ Ω : θ(x * , t) = 0 which indicates the extent of frostbite.
Simulations for Example 3 confirm intuition. Upon vasoconstriction (b-c), frostbite is more extensive than in (a), but not dramatically. There is not much difference in θ between (b) and (c), but there is difference in θ body between (b-c).
There is little dependence on c b if it remains below 10 4 , perhaps because c vessel is fixed. As expected, higher temperatures and less frostbite occur when vasoconstriction is not active; see also Figure 3 for comparison of θ vaso and θ no vaso .
Summary: Based on the computational experiments, we believe the CCFD implementation of the H model is robust to model variants, including the use of fictitious domain. The results agree qualitatively with the intuition; the most crucial parameters are c b and c vessel .
Consider a complex domain shown in Figure 1 (b) representing a plant root embedded in soil, adapted from [23]. Our long-term goal is to simulate water flow in this root-soil system combined with other coupled phenomena including surface and above-surface models plus energy equation. Therefore, even though roots have a much smaller volume than reasonable soil volumes, we aim to build a general flexible physically meaningful RS model in d ≥ 1.
Our RS model extends the d = 1 overlapping continua root-soil models from [27,4], a nonlinear Richards-Darcy generalization of (3). We consider Ω r ⊂ Ω s ⊂ cD c b vaso # frostbite # off θmin θmax θave
10 6 10 -3 yes 42 165 -16.03 36.99 27.28 10 6 10 -3 no 42 165 -16.03 36.99 27.28 10 6 10 2 yes 42 163 -16.04 36.99 27.29 10 6 10 2 no 42 163 -15.98 36.99 27.29 10 6 10 4 yes 32 131 -15.37 36.99 28.18 10 6 10 4 no 22 124 -11.81 36.99 28.28 10 6 10 5 yes 1 3 -0.02 36.98 32.76 10 6 10 5 no 0 2 8.38 36.98 32.78 10 7 10 4 yes 6 40 -9.03 36.99 32.66 10 7 10 4 no 5 39 -6.65 36.99 32.70 10 7 10 5 yes 0 2 0.87 37 34.46 10 7 10 5 no 0 2 8.70 37 34.47 Table 2. Sensitivity of the solutions to the H model with vasoconstriction in Example 4 (2D) at t = T = 120 s to c b , cD and to the choice of vasoconstriction model (vaso or off). We define # frostbite as the number of the cells ωij where θij|t=T < 0, and # off those with θij|t=T < 10. Richards equation in unsaturated soil domain Ω s models (incompressible) water flow in Ω s coupled to saturated flow in Ω r as follows
∇ • q r = c(P s -P r ), x ∈ Ω r , t > 0, (7b)
where we denote the soil/root domain by subscripts s / r , the saturation (volume fraction) of water in the soil by S s , the pressure of the fluid in Ω m by P m . In soil we have P s = -P c (S s ), where P c is capillary pressure. In addition, q is flux of the fluid, ϕ and K are porosity and permeability of the medium, k(S) is the relative permeability, µ and ρ are viscosity and density of the fluid, g is the gravitational acceleration, D is the depth under the soil's surface, and c ∼ χ Ωr is the coefficient of the exchange term between root and soil. Of main interest is the nonnegative exchange coefficient c = c(x) with supp c = Ω r . Suppose now that the initial and general boundary conditions are
When d = 1, the model ( 7) is derived from first principles in [4] and independently in [27], making several modeling assumptions, and sharing principles.
Model in [4] emphasizes the geometrical complexity of the root-soil exchange, and works in potentials Ps g as variables. The idea is to set up the fluid exchange between root and soil in a simple manner and avoid discretization at the scale of individual xylems through which the actual water transport takes place. To this end, the radial character of the flow to the roots, the low conductivity of the inner portion of the roots (endodermis), along with high permeability of epidermis and medium permeability of cortex, are assumed. The models are of overlapping continua type (3) but nonlinear and without the overall periodic assumption. In [27] this model is extended to allow subroots and root network with stochastic updates to this d = 1 model, and predict that dry and wet zones will develop in the soil as a result of water uptake by plant roots.
In our approach we work directly with a complex root geometry within a d ≥ 1 domain.
First we consider c = 0. The main well known difficulty when working with Richards equation in any dimension is the nonlinear degenerate parabolic character of the problem due to the behavior of k(S), P c (S). We illustrate this behavior using the well known algebraic model based on experimental data known as van-Genuchten-Mualem model
The set of parameters ϵ = 1 2 , m = 1 -1 ν , α = 10 -4 , and ν = 2.237 characterizes fine soil; see this and others collected in [22]. A change of variables in (7a) reveals that an equivalent model reads
Here D s (S) ≥ 0 but degenerates to 0 when S ↓ 0. Thus the model features a degenerate behavior, while the nonlinear advective term A(S) when D ̸ = const requires careful treatment of this hyperbolic term. Due to these difficulties, the use of lower order numerical scheme such as that in Section 2.1 is appropriate, and so is the use of upwinding. Furthermore, the choice of primary unknowns is delicate, since P c (S) is unbounded near S ↓ 0, but P -1 c (p) is unbounded when p ↓ 0. These features prompted various analyses of numerical schemes and delicate discussion of nonlinear solvers; see, e.g., [5,25,21,22]. In particular, [5,25] prove O(h + τ ) error estimates to the mixed finite element method, but some of this work requires regularization and strong assumptions on smoothness. In turn, [22] evaluates the fully implicit CCFD schemes for nontransformed, nonregularized models with strong heterogeneities and locally large fluxes, and test variants of averaging, implicit or semi-implicit solutions, and different choices of primary unknowns.
When c > 0, the computational model for RS system features an additional difficulty, since most likely |Ω r | << |Ω s |. With the approach we outlined in Section 2.1, this is, however, not an issue. We are able to simulate successfully the flow of water in the overlapping continua system; this is illustrated by our examples below.
[s] [m] Value 0.4 10 -11 5 × 10 -12 10 -3 10 3 9.8066 10 -10 600 1/50 Table 3. Parameter values for Example 5.
Example 5. We consider Ω r = Ω s = (0, L) and L = 1 m, other data in Table 3, and k(S), P c (S) given by (8). To model evaporation and precipitation at the surface of the ground, we impose Neumann boundary conditions at x = 0, and at x = L: q s (0, t) = q s,top , q r (0, t) = q r,top , q s (1, t) = 0, q r (1, t) = 0, ( 9)
The simulation results in Figure 4 show response of the system to the infiltration through boundary (9). Example 6. We continue with the physical data from Example 5 but in Ω s = [0, 1] × [0, 1] and grid 120 × 80 over t ≤ T = 3600 s. We simulate infiltration with the Richards-Darcy model, starting from initial (pressure-saturation) equilibrium (by setting q s = q r = 0, P s = P r , and P s,top = P r,top = 0), adding water only from left top portion of the domain, and no flux boundary conditions elsewhere. We use the coupling coefficient c = 10 -8 .
The simulations results plotted in Figure 5 confirm that the pressure in root is well equilibrated with P s . The water from the left boundary eventually reaches the root. The structure of the solutions features sharp fronts and, as usual for Richards equation, is very rich. The system strongly depends on c. Pr|t=2400 Ps|t=2400 Ss|t=2400
Fig. 5. Numerical solution for Pr, Ps, and Ss, from Example 6. For Ps, Ss, the contour ∂Ωr is superimposed over Ωs to guide the eye. Vertical axis is aligned with gravity.
Summary: We find the RS model based on Richards-Darcy equations quite complex. As long as its Richards equation is well calibrated, the model is robust. In the coupled system, there is high sensitivity to the coefficient c.
Traffic flow models can have continuum (PDE) or discrete (individual based model) form. The former have gained interest since the LWR models [17] and present interesting hyperbolic structure of the underlying PDE, the (macroscopic) transport model u t + ∇ • f (u) = 0 (supplemented by inflow boundary conditions) for the average density u (of cars) with flux function f (u), e.g., f (u) = (1 -u umax )uv. However, these traffic models for cars are inherently one-dimensional, since the vehicles can only travel along trajectories confined to paved roadways. In addition, the spatial scale is much larger than the average length of of car, individualistic behavior is not preserved. In contrast, microsopic models do account for individualistic behavior, but come with a high computational cost and complexity. Individual based models are important for emergency scenarios such as during tsunami or wildfires as well as for urban and architectural design, to determine possible congestion patterns. Large network traffic models can become unmanageable in its complexity, and efforts to manage involve upscaling to a model resembling "Darcy flow" [9].
Computational models of traffic scenarios start with an image such as in Figure 1 (c) which describes the network subset of a domain Ω ⊂ R 2 . One possible model for traffic on this network augments the LWR flux with -K∇u, where K is a diffusion tensor and v captures the preferred direction of motion; both are based on the mean value and standard deviation of transition probabilities between sites at each time step [9]. However, this approach does not include coupling between the different species nor allow travel off the network.
The model we build accounts for the traffic on and off the network Ω P as well as nonlinear interactions between the species. Our simulations might aid in the optimal design and control of robot trajectories as well as of campus pathway design and monitoring.
We consider a traffic problem involving several species numbered m = 1, 2, . . . M . We let u = (u 1 , . . . u M ) and f (u) = (f 1 (u), . . . f M (u)), where u is the vector of average densities of each species and f is the flux function. The M = 3 species are human pedestrians, humans on bicycles, and autonomous delivery robots. We set f m (x, t; u) = v m (x, t; u) fm (u). to separate the trajectories v m (x, t; u) from traffic awareness modelled by fm (u). While robots or humans on bicycles can only travel on paved paths within some Ω P ⊂ Ω, others can alter their trajectory and veer off Ω P when they are aware that the current traffic situation leads to a (possible) congestion.
Flow and trajectories: Each species m has a well defined fixed trajectory given by the velocity field v m = v m (x), x ∈ Ω. Most can only move on the paved paths, i.e. so that supp v m ⊂ Ω P . Typically most species are aware of others and may alter their speed locally in time, but most cannot alter their trajectories. For some species, say m ′ , we allow v
To determine a trajectory v m , assume that species m intends to get from some inlet point x in m ∈ ∂Ω to some outlet point x out m ∈ ∂Ω. We solve a pseudopotential problem ∇ • v m = 0 for ψ m so that v m = -K m ∇ψ m . We set ∂Ω D = Γ in ∪Γ out where each Γ in , Γ out is a small region of ∂Ω around the inlet and outlet, respectively, and require the Dirichlet condition ψ| Γin = 1, and ψ| Γout = 0, and we require Neumann condition v m • n| ∂Ω N = 0 on ∂Ω N = ∂Ω \ ∂Ω D . Solving for v m is a well-posed problem with a scheme for "flow" from Section 2.1.
It remains to set up the species specific "mask" K to determine these paths preferentially and allow, e.g., to set up alternative pathways to stairs, or set up ad-hoc detours. In the simplest scenarios we simply set K| Ω P = κ m and K| Ω\Ω P = 0, and κ m can be calibrated to an average speed of species m. Also, the species "aware" of traffic are allowed to alter their trajectory with K m = K m (x, y, (u m ) m ) built heuristically. For example, if at some point (x, y) ∈ Ω P we have a traffic congestion in some area Ω C , we would set
The transport model: With f m = v m fm and a given v m , fm (u) can be a linear or nonlinear flux function, e.g., the LWR model.
Scheme: To approximate the flow and transport solutions, we follow Section 2.1 and [8,17] and apply Godunov's method, under CFL condition. For M > 1 the situation is more complicated but for the simple case studies we develop the Godunov method suffices. For the simulations in d = 1, we report grid refinement studies; we also confirm grid error for the nonlinear problems to be O( √ h) (not shown).
To build up our intuition, we illustrate the coupled dynamics of the M = 2 species with Ω = (0, 2), assuming they travel in the same direction. The pedestrian density is u 1 , and that of robots is u 2 . We construct the flux functions heuristically to model the interaction between the species, and in particular "traffic awareness". Species 1 (humans) do not feature awareness of other species and f1 (u) = u 1 (1 -u 1 ) (LWR model), but f2 (u) = v 2 (u 1 )u 2 where the robots slow down considerably when they notice humans, with Example 7. We prescribe u init 1 (x) = χ [0.6,0.7] (x) and u init 2 (x) = 0.1χ [0.4,0.8] to illustrate the traffic flow and "traffic awareness": the robots start ahead, alongside and behind human pedestrians on the 1d network; see Figure 6.
We see that u 1 feature a rarefaction typical in LWR traffic flow models [17]. Also, the robots follow linear advection away from pedestrians, but develop two traveling waves for robots ahead and behind the humans according to (10). The snapshot at t = 0.75 shows that u 2 | [0.5,0.65] matches linear advection, but when u 1 > 0.6, the robots slow down to v 2 (u 1 ) = 0.1 causing a large spike.
Example 8. Now we consider the campus network with Ω P shown in Figure 1 (c). We simulate mixed dimensional network traffic on Ω P with the occasional traffic off Ω P . We design trajectories of the species as shown in Figure 7, and we simulate the transport as shown in Figure 8. In one variant, (a) humans stay on pavement. In another, (b) humans venture off the path to avoid congestion. (a) (b) Fig. 8. Results of Example 8 at time t = 0.5 show the average human concentration u1 in case (a). In case (b) due to congestion by u2 and u3 the pedestrians alter v1 off the network Ω P to avoid the bottleneck.
Summary: The computational framework of overlapping continua allows to simulate complex traffic patterns involving species of very different trajectories. More work is needed to identify the data for such simulations and applications to the design of campus network and of, e.g., robot software.
The implementation for each of the applications H, RS, and T is efficient and based on a common robust CCFD framework which allows for approximate enforcement of Dirichlet conditions and easy accounting for complex domains. The model approximations and enhancements were easily introduced; we also identified critical model components and parameters. For H, it is the c b coefficient, and its role in the vasoconstriction. For RS, it is by far also the exchange coefficient c. For T, it is the local velocity, the mask K, and the interaction models between the species. These elements should be validated with data.
However, while parameter identification based on experimental or imaging data is fairly easy for models involving scalar unknowns, it remains quite challenging for coupled systems, where identifying proper relationships and models can be difficult. The prototype models that we presented can serve as useful proof-of-concept first steps to identify the main features of dynamics.