Add lazy-core builders for all 30 COPS instances

ext/COPSBenchmarkExaModels/COPSBenchmarkExaModels.jl

@@ -5,6 +5,12 @@ import COPSBenchmark
 import COPSBenchmark: ExaModelsBackend
 using ExaModels
 
+# `Deferred(f)` (not exported): a function of the resolved args, evaluated
+# once per materialization — used by the *_core builders as a value slot
+# (computed start arrays, bounds) or as an index set (pars records carrying
+# precomputed index data per element).
+const Deferred = ExaModels.Deferred
+
 include("bearing.jl")
 include("camshape.jl")
 include("catmix.jl")

ext/COPSBenchmarkExaModels/bearing.jl

@@ -42,3 +42,42 @@ end
 
 
 
+
+function COPSBenchmark.bearing_core()
+    args = ExaModels.ArgTracer()
+    c = ExaCore(concrete = Val(true))
+    b = 10              # grid is (0,2*pi)x(0,2*b)
+    e = 0.1             # eccentricity
+
+    hx = 2 * pi / (args.nx + 1)  # grid spacing
+    hy = 2 * b / (args.ny + 1)   # grid spacing
+    area = 0.5 * hx * hy         # area of triangle
+
+    wq(i) = (1.0 + e * cos((i - 1) * hx))^3
+    v0 = [max(sin((i - 1) * hx), 0.0) for i in 1:args.nx+2, j in 1:args.ny+2]
+
+    c, v = add_var(c, 1:args.nx+2, 1:args.ny+2; lvar = 0.0, start = v0)
+
+    c, _ = add_obj(
+        c,
+        0.5 * (hx * hy / 6.0) * (wq(i) + 2 * wq(i + 1)) *
+        (((v[i+1, j] - v[i, j]) / hx)^2 + ((v[i, j+1] - v[i, j]) / hy)^2)
+        for i in 1:args.nx+1, j in 1:args.ny+1
+    )
+    c, _ = add_obj(
+        c,
+        0.5 * (hx * hy / 6.0) * (2 * wq(i) + 2 * wq(i - 1)) *
+        (((v[i-1, j] - v[i, j]) / hx)^2 + ((v[i, j-1] - v[i, j]) / hy)^2)
+        for i in 2:args.nx+2, j in 2:args.ny+2
+    )
+    c, _ = add_obj(
+        c,
+        -hx * hy * e * sin((i - 1) * hx) * v[i, j] for i in 1:args.nx+2, j in 1:args.ny+2
+    )
+
+    c, _ = add_con(c, v[i, 1] for i in 1:args.nx+2)
+    c, _ = add_con(c, v[i, args.ny+2] for i in 1:args.nx+2)
+    c, _ = add_con(c, v[1, i] for i in 1:args.ny+2)
+    c, _ = add_con(c, v[args.nx+2, i] for i in 1:args.ny+2)
+    c
+end

ext/COPSBenchmarkExaModels/camshape.jl

@@ -66,3 +66,63 @@
     return ExaModels.ExaModel(core; kwargs...)
 end
 
+
+function COPSBenchmark.camshape_core()
+    args = ExaModels.ArgTracer()
+    c = ExaCore(minimize = false, concrete = Val(true))
+    R_v = 1.0         # design parameter related to the valve shape
+    R_max = 2.0       # maximum allowed radius of the cam
+    R_min = 1.0       # minimum allowed radius of the cam
+    alpha = 1.5       # curvature limit parameter
+
+    d_theta = 2 * pi / (5 * (args.n + 1))   # angle between discretization points
+
+    # radius of the cam at discretization points
+    c, r = add_var(c, 1:args.n; lvar = R_min, uvar = R_max, start = (R_min + R_max) / 2.0)
+
+    c, _ = add_obj(c, (pi * R_v) / args.n * r[i] for i in 1:args.n)
+
+    # Convexity
+    c, _ = add_con(
+        c,
+        (-r[i-1] * r[i] - r[i] * r[i+1] + 2 * r[i-1] * r[i+1] * cos(d_theta) for i = 2:args.n-1);
+        lcon = -Inf, ucon = 0.0,
+    )
+    c, _ = add_con(
+        c,
+        (-R_min * r[1] - r[1] * r[2] + 2 * R_min * r[2] * cos(d_theta) for _ in 1:1);
+        lcon = -Inf, ucon = 0.0,
+    )
+    c, _ = add_con(
+        c,
+        (-R_min^2 - R_min * r[1] + 2 * R_min * r[1] * cos(d_theta) for _ in 1:1);
+        lcon = -Inf, ucon = 0.0,
+    )
+    c, _ = add_con(
+        c,
+        (-r[m-1] * r[m] - r[m] * R_max + 2 * r[m-1] * R_max * cos(d_theta) for m in args.n:args.n);
+        lcon = -Inf, ucon = 0.0,
+    )
+    c, _ = add_con(
+        c,
+        (-2 * R_max * r[m] + 2 * r[m]^2 * cos(d_theta) for m in args.n:args.n);
+        lcon = -Inf, ucon = 0.0,
+    )
+    # Curvature
+    c, _ = add_con(
+        c,
+        ((r[i+1] - r[i]) for i = 1:args.n-1);
+        lcon = -alpha * d_theta, ucon = alpha * d_theta,
+    )
+    c, _ = add_con(
+        c,
+        ((r[1] - R_min) for _ in 1:1);
+        lcon = -alpha * d_theta, ucon = alpha * d_theta,
+    )
+    c, _ = add_con(
+        c,
+        ((R_max - r[m]) for m in args.n:args.n);
+        lcon = -alpha * d_theta, ucon = alpha * d_theta,
+    )
+    c
+end

ext/COPSBenchmarkExaModels/catmix.jl

@@ -79,3 +79,70 @@
 
     return ExaModels.ExaModel(c; kwargs...)
 end
+
+function COPSBenchmark.catmix_core()
+    args = ExaModels.ArgTracer()
+    c = ExaCore(concrete = Val(true))
+    ne = 2
+    nc = 3
+
+    h = 1 / args.nh   # Final time / nh
+
+    rho = [
+        0.11270166537926,
+        0.50000000000000,
+        0.88729833462074,
+    ]
+    bc = [1.0, 0.0]   # Boundary conditions for x
+    alpha = 0.0       # Smoothing parameter
+    neg_one = -1.0
+    ten_T = 10.0
+    one_T = 1.0
+    rho_index = [(i, rho[i]) for i in 1:nc]
+
+    # lvar/uvar/start given as zeros(T,...)/ones(T,...) in the eager
+    # constructor: written as broadcast scalars here (value-identical).
+    c, u = add_var(c, args.nh, nc; lvar = 0.0, uvar = 1.0, start = 0.0)
+    c, v = add_var(c, args.nh, ne; start = [mod(j, ne) for i in 1:args.nh, j in 1:ne])
+    c, w = add_var(c, args.nh, nc, ne; start = 0.0)
+    c, pp = add_var(c, args.nh, nc, ne; start = [mod(k, ne) for i in 1:args.nh, j in 1:nc, k in 1:ne])
+    c, Dpp = add_var(c, args.nh, nc, ne; start = 0.0)
+    c, ppf = add_var(c, ne; start = [mod(i, ne) for i in 1:ne])
+
+    c, _ = add_obj(c, neg_one + ppf[1] + ppf[2] for _ in 1:1)
+    c, _ = add_obj(c, alpha / h * (u[i+1, j] - u[i, j])^2 for i in 1:args.nh-1, j in 1:nc)
+
+    c, _ = add_con(
+        c,
+        pp[i, k, s] - v[i, s] - h * sum(w[i, j, s] * (rho^j / factorial(j)) for j in 1:nc)
+        for i = 1:args.nh, (k, rho) in rho_index, s = 1:ne
+    )
+    c, _ = add_con(
+        c,
+        Dpp[i, k, s] - sum(w[i, j, s] * (rho^(j-1) / factorial(j - 1)) for j in 1:nc)
+        for i = 1:args.nh, (k, rho) in rho_index, s = 1:ne
+    )
+    c, _ = add_con(
+        c,
+        ppf[s] - v[args.nh, s] - h * sum(w[args.nh, j, s] / factorial(j) for j in 1:nc) for s in 1:ne
+    )
+    c, _ = add_con(
+        c,
+        v[i, s] + sum(w[i, j, s] * h / factorial(j) for j in 1:nc) - v[i+1, s]
+        for i in 1:args.nh-1, s in 1:ne
+    )
+    c, _ = add_con(
+        c,
+        Dpp[i, j, 1] - u[i, j] * (ten_T * pp[i, j, 2] - pp[i, j, 1]) for i = 1:args.nh, j = 1:nc
+    )
+    c, _ = add_con(
+        c,
+        Dpp[i, j, 2] - u[i, j] * (pp[i, j, 1] - ten_T * pp[i, j, 2]) + (one_T - u[i, j]) * pp[i, j, 2]
+        for i = 1:args.nh, j = 1:nc
+    )
+    c, _ = add_con(
+        c,
+        v[1, s] - bc for (s, bc) in [(i, bc[i]) for i in 1:ne]
+    )
+    c
+end

ext/COPSBenchmarkExaModels/chain.jl

@@ -78,3 +78,34 @@
     return ExaModels.ExaModel(c; kwargs...)
 end
 
+
+function COPSBenchmark.chain_core()
+    args = ExaModels.ArgTracer()
+    c = ExaCore(concrete = Val(true))
+    nh = max(2, div(args.n - 4, 4))
+
+    L = 4
+    a = 1
+    b = 3
+    tmin = b > a ? 1 / 4 : 3 / 4
+    tf = 1.0
+    h = tf / nh
+
+    c, u = add_var(c, nh + 1; start = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1])
+    c, x1 = add_var(c, nh + 1; start = [4 * abs(b - a) * k / nh * (1 / 2 * k / nh - tmin) + a for k in 1:nh+1])
+    c, x2 = add_var(c, nh + 1; start = [(4 * abs(b - a) * k / nh * (1 / 2 * k / nh - tmin) + a) *
+        (4 * abs(b - a) * (k / nh - tmin)) for k in 1:nh+1])
+    c, x3 = add_var(c, nh + 1; start = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1])
+
+    c, _ = add_obj(c, x2[m] for m in nh+1:nh+1)
+
+    c, _ = add_con(c, x1[j+1] - x1[j] - 1 / 2 * h * (u[j] + u[j+1]) for j in 1:nh)
+    c, _ = add_con(c, x1[1] - a for _ in 1:1)
+    c, _ = add_con(c, x1[m] - b for m in nh+1:nh+1)
+    c, _ = add_con(c, x2[1] for _ in 1:1)
+    c, _ = add_con(c, x3[1] for _ in 1:1)
+    c, _ = add_con(c, x3[m] - L for m in nh+1:nh+1)
+    c, _ = add_con(c, x2[j+1] - x2[j] - 1 / 2 * h * (x1[j] * sqrt(1 + u[j]^2) + x1[j+1] * sqrt(1 + u[j+1]^2)) for j in 1:nh)
+    c, _ = add_con(c, x3[j+1] - x3[j] - 1 / 2 * h * (sqrt(1 + u[j]^2) + sqrt(1 + u[j+1]^2)) for j in 1:nh)
+    c
+end

ext/COPSBenchmarkExaModels/channel.jl

@@ -134,3 +134,130 @@
 
     return ExaModels.ExaModel(core; kwargs...)
 end
+
+function COPSBenchmark.channel_core()
+    args = ExaModels.ArgTracer()
+    nc = 4
+    nd = 4
+    R = 10.0
+    tf = 1.0
+    rho = [0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703]
+
+    # v0 (channel.jl:20-26): mutating fill over 1:nh from the partition t
+    _v0 = a -> begin
+        h = tf / a.nh
+        t = [(i-1)*h for i in 1:a.nh+1]
+        v0 = zeros(a.nh, nd)
+        for i in 1:a.nh
+            v0[i, 1] = t[i]^2*(3 - 2*t[i])
+            v0[i, 2] = 6*t[i]*(1 - t[i])
+            v0[i, 3] = 6*(1 - 2*t[i])
+            v0[i, 4] = -12
+        end
+        v0
+    end
+
+    # fac[k+1] = k!
+    _fac = [factorial(k) for k in 0:nc+nd]
+
+    con1_itr = Deferred(a -> begin
+        h = tf / a.nh
+        fac = _fac
+        [(i, j, s, h*rho[j]^1/fac[2], h*rho[j]^2/fac[3], h*rho[j]^3/fac[4], h*rho[j]^4/fac[5])
+            for i in 1:a.nh, j in 1:nc, s in 1:nd]
+    end)
+    con2_itr = Deferred(a -> begin
+        h = tf / a.nh
+        fac = _fac
+        [(i, j, s,
+          (s <= 1 ? (rho[j]*h)^(1-s)/fac[1-s+1] : 0.0),
+          (s <= 2 ? (rho[j]*h)^(2-s)/fac[2-s+1] : 0.0),
+          (s <= 3 ? (rho[j]*h)^(3-s)/fac[3-s+1] : 0.0),
+          (s <= 4 ? (rho[j]*h)^(4-s)/fac[4-s+1] : 0.0),
+          h^(nd-s+1)*rho[j]^(1+nd-s)/fac[1+nd-s+1],
+          h^(nd-s+1)*rho[j]^(2+nd-s)/fac[2+nd-s+1],
+          h^(nd-s+1)*rho[j]^(3+nd-s)/fac[3+nd-s+1],
+          h^(nd-s+1)*rho[j]^(4+nd-s)/fac[4+nd-s+1])
+            for i in 1:a.nh, j in 1:nc, s in 1:nd]
+    end)
+    cont_itr = Deferred(a -> begin
+        h = tf / a.nh
+        fac = _fac
+        [(i, s,
+          (s <= 1 ? h^(1-s)/fac[1-s+1] : 0.0),
+          (s <= 2 ? h^(2-s)/fac[2-s+1] : 0.0),
+          (s <= 3 ? h^(3-s)/fac[3-s+1] : 0.0),
+          (s <= 4 ? h^(4-s)/fac[4-s+1] : 0.0),
+          h^(nd-s+1)/fac[1+nd-s+1],
+          h^(nd-s+1)/fac[2+nd-s+1],
+          h^(nd-s+1)/fac[3+nd-s+1],
+          h^(nd-s+1)/fac[4+nd-s+1])
+            for i in 1:a.nh-1, s in 1:nd]
+    end)
+    coll_itr = Deferred(a -> begin
+        fac = _fac
+        [(i, j, rho[j]^0/fac[1], rho[j]^1/fac[2], rho[j]^2/fac[3], rho[j]^3/fac[4])
+            for i in 1:a.nh, j in 1:nc]
+    end)
+    # right-BC coefficient records (channel.jl:79-82), one row each
+    bc3_itr = Deferred(a -> begin
+        h = tf / a.nh
+        fac = _fac
+        cv = [h^(k-1)/fac[k] for k in 1:nd]
+        cw = [h^nd/fac[k+nd] for k in 1:nc]
+        [(cv1 = cv[1], cv2 = cv[2], cv3 = cv[3], cv4 = cv[4],
+          cw1 = cw[1], cw2 = cw[2], cw3 = cw[3], cw4 = cw[4])]
+    end)
+    bc4_itr = Deferred(a -> begin
+        h = tf / a.nh
+        fac = _fac
+        cv = [h^(k-2)/fac[k-1] for k in 2:nd]
+        cw = [h^(nd-1)/fac[k+nd-1] for k in 1:nc]
+        [(cv1 = cv[1], cv2 = cv[2], cv3 = cv[3],
+          cw1 = cw[1], cw2 = cw[2], cw3 = cw[3], cw4 = cw[4])]
+    end)
+
+    core = ExaCore(concrete = Val(true))
+
+    core, v = add_var(core, args.nh, nd; start = Deferred(a -> _v0(a)))
+    core, w = add_var(core, args.nh, nc; start = 0.0)
+    core, uc = add_var(core, args.nh, nc, nd;
+        start = Deferred(a -> begin v0 = _v0(a); [v0[i, s] for i in 1:a.nh, j in 1:nc, s in 1:nd] end))
+    core, Duc = add_var(core, args.nh, nc, nd; start = 0.0)
+
+    # Constant objective
+    core, _ = add_obj(core, 1.0 for _i in 1:1)
+
+    core, _ = add_con(core,
+        uc[i, j, s] - v[i, s] - (a1*w[i,1] + a2*w[i,2] + a3*w[i,3] + a4*w[i,4])
+        for (i, j, s, a1, a2, a3, a4) in con1_itr)
+    core, _ = add_con(core,
+        Duc[i, j, s] - (bv1*v[i,1] + bv2*v[i,2] + bv3*v[i,3] + bv4*v[i,4])
+                     - (bw1*w[i,1] + bw2*w[i,2] + bw3*w[i,3] + bw4*w[i,4])
+        for (i, j, s, bv1, bv2, bv3, bv4, bw1, bw2, bw3, bw4) in con2_itr)
+
+    # Boundary conditions (bc = [0.0 1.0; 0.0 0.0])
+    core, _ = add_con(core, v[1, 1] - 0.0 for _ in 1:1)
+    core, _ = add_con(core, v[1, 2] - 0.0 for _ in 1:1)
+    core, _ = add_con(core,
+        p.cv1*v[args.nh,1] + p.cv2*v[args.nh,2] + p.cv3*v[args.nh,3] + p.cv4*v[args.nh,4] +
+        p.cw1*w[args.nh,1] + p.cw2*w[args.nh,2] + p.cw3*w[args.nh,3] + p.cw4*w[args.nh,4] - 1.0
+        for p in bc3_itr)
+    core, _ = add_con(core,
+        p.cv1*v[args.nh,2] + p.cv2*v[args.nh,3] + p.cv3*v[args.nh,4] +
+        p.cw1*w[args.nh,1] + p.cw2*w[args.nh,2] + p.cw3*w[args.nh,3] + p.cw4*w[args.nh,4] - 0.0
+        for p in bc4_itr)
+
+    # Continuity
+    core, _ = add_con(core,
+        cv1*v[i,1] + cv2*v[i,2] + cv3*v[i,3] + cv4*v[i,4] +
+        cw1*w[i,1] + cw2*w[i,2] + cw3*w[i,3] + cw4*w[i,4] - v[i+1, s]
+        for (i, s, cv1, cv2, cv3, cv4, cw1, cw2, cw3, cw4) in cont_itr)
+
+    # Collocation physics
+    core, _ = add_con(core,
+        e1*w[i,1] + e2*w[i,2] + e3*w[i,3] + e4*w[i,4] -
+        R * (Duc[i, j, 2] * Duc[i, j, 3] - Duc[i, j, 1] * Duc[i, j, 4])
+        for (i, j, e1, e2, e3, e4) in coll_itr)
+    core
+end

ext/COPSBenchmarkExaModels/elec.jl

@@ -26,3 +26,34 @@
 
     return ExaModels.ExaModel(core; kwargs...)
 end
+
+function COPSBenchmark.elec_core()
+    args = ExaModels.ArgTracer()
+    c = ExaCore(concrete = Val(true))
+    c, x = add_var(c, 1:args.np; start = Deferred(a -> begin
+        COPSBenchmark.Random.seed!(2713)
+        theta = (2pi) .* rand(a.np)
+        phi = pi .* rand(a.np)
+        [cos(theta[i])*sin(phi[i]) for i = 1:a.np]
+    end))
+    c, y = add_var(c, 1:args.np; start = Deferred(a -> begin
+        COPSBenchmark.Random.seed!(2713)
+        theta = (2pi) .* rand(a.np)
+        phi = pi .* rand(a.np)
+        [sin(theta[i])*sin(phi[i]) for i = 1:a.np]
+    end))
+    c, z = add_var(c, 1:args.np; start = Deferred(a -> begin
+        COPSBenchmark.Random.seed!(2713)
+        theta = (2pi) .* rand(a.np)
+        phi = pi .* rand(a.np)
+        [cos(phi[i]) for i = 1:a.np]
+    end))
+
+    # Coulomb potential
+    itr = Deferred(a -> [(i, j) for i in 1:a.np-1 for j in i+1:a.np])
+    c, _ = add_obj(c, 1.0 / sqrt((x[i] - x[j])^2 + (y[i] - y[j])^2 + (z[i] - z[j])^2) for (i, j) in itr)
+
+    # Unit-ball
+    c, _ = add_con(c, x[i]^2 + y[i]^2 + z[i]^2 - 1 for i = 1:args.np)
+    c
+end