[COPS] add PDE instances

Project.toml

@@ -9,6 +9,7 @@ ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+LuksanVleckBenchmark = "dc0393f0-98a3-4386-a65b-465e1307889b"
 
 [compat]
 CUDA = "5.2.0"

data/tetra.jl

@@ -0,0 +1,25 @@
+xe_tetra = [
+    0 0 0
+    1 0 0
+    0 1 0
+    0 0 1
+    5//10 2//10 1//10
+];
+
+Tets_tetra = [
+    1 2 3 5
+    5 4 3 1
+    5 4 1 2
+    5 2 3 4
+];
+
+Tets_tetra = vec(reshape(Tets_tetra, 16, 1));
+xe_tetra = vec(reshape(xe_tetra, 15, 1));
+Constants_tetra = [1, 2, 3, 4];
+
+include("tetra_duct12.jl")
+include("tetra_duct15.jl")
+include("tetra_duct20.jl")
+include("tetra_hook.jl")
+include("tetra_foam5.jl")
+include("tetra_gear.jl")

data/triangle.jl

@@ -0,0 +1,20 @@
+xe = [
+    0 0
+    1 0
+    0 1
+    64//100 25//100
+];
+
+Tr = [
+    1 2 4
+    4 2 3
+    4 3 1
+];
+
+Tr = vec(reshape(Tr, 9, 1));
+xe = vec(reshape(xe, 8, 1));
+Constants = [1, 2, 3];
+
+include("triangle_deer.jl")
+include("triangle_pacman.jl")
+include("triangle_turtle.jl")

src/ExaModelsExamples.jl

@@ -5,6 +5,7 @@ import Downloads
 import ExaModels: ExaModels, NLPModels
 import PowerModels: PowerModels, silence
 import Random
+import LuksanVleckBenchmark
 
 include("opf.jl")
 include("scopf.jl")
@@ -15,13 +16,24 @@ include("quadrotor.jl")
 include("goddard.jl")
 include("robot.jl")
 include("rocket.jl")
-include("bearing.jl") 
+include("bearing.jl")
 include("camshape.jl")
 include("elec.jl")
 include("steering.jl")
 include("pinene.jl")
 include("marine.jl")
 include("gasoil.jl")
+include("pde_models.jl")
+include("catmix.jl")
+include("chain.jl")
+include("channel.jl")
+include("glider.jl")
+include("minsurf.jl")
+include("polygon.jl")
+include("tetra.jl")
+include("torsion.jl")
+include("triangle.jl")
+include("methanol.jl")
 
 const NAMES = filter(names(ExaModelsExamples; all = true)) do x
     str = string(x)

src/bearing.jl

@@ -5,7 +5,7 @@
 # COPS 3.1 - March 2004
 
 
-function bearing_model(nx, ny; T = Float64, backend = nothing, kwargs...)
+function bearing_model(nx, ny; T = Float64, backend = CUDABackend(), kwargs...)
     b = 10              # grid is (0,2*pi)x(0,2*b)
     e = 0.1             # eccentricity
 
@@ -38,7 +38,4 @@ function bearing_model(nx, ny; T = Float64, backend = nothing, kwargs...)
     ExaModels.constraint(core, v[nx+2, i] - 0.0 for i in 1:nx+2)
 
     return ExaModels.ExaModel(core; kwargs...)
-end
-
-
-
+end

src/catmix.jl

@@ -0,0 +1,74 @@
+# Catalyst Mixing Problem
+# Collocation formulation
+# COPS 3.0 - November 2002
+# COPS 3.1 - March 2004
+
+function catmix_model(nh; T = Float64, backend = nothing, kwargs...)
+    ne = 2
+    nc = 3
+
+    tf = 1
+    h = tf / nh   # Final time
+
+    rho = [
+        0.11270166537926,
+        0.50000000000000,
+        0.88729833462074,
+    ]
+    bc = [1.0, 0.0]  # Boundary conditions for x
+    alpha = 0.0      # Smoothing parameter
+    rho_index = [(i, rho[i]) for i in 1:nc]
+
+    c = ExaModels.ExaCore(T; backend = backend)
+    u = ExaModels.variable(c, nh, nc; lvar = zeros(nh, nc), uvar = ones(nh, nc), start = zeros(nh, nc))
+    v = ExaModels.variable(c, nh, ne; start = [mod(j, ne) for i in 1:nh, j in 1:ne])
+    w = ExaModels.variable(c, nh, nc, ne; start = zeros(nh, nc, ne))
+    pp = ExaModels.variable(c, nh, nc, ne; start = [mod(k, ne) for i in 1:nh, j in 1:nc, k in 1:ne])
+    Dpp = ExaModels.variable(c, nh, nc, ne; start = zeros(nh, nc, ne))
+    ppf = ExaModels.variable(c, ne; start = [mod(i,ne) for i in 1:ne])
+
+    ExaModels.objective(c, -1.0 + ppf[1] + ppf[2])
+    ExaModels.objective(c, alpha/h*(u[i+1, j] - u[i, j])^2 for i in 1:nh-1, j in 1:nc)
+
+    ExaModels.constraint(
+        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:nh,  (k, rho) in rho_index, s=1:ne
+    )
+
+    ExaModels.constraint(
+        c,
+        Dpp[i, k, s] - sum(w[i, j, s]*(rho^(j-1)/factorial(j-1)) for j in 1:nc) for i=1:nh, (k, rho) in rho_index, s=1:ne
+    )
+
+    ExaModels.constraint(
+        c,
+        ppf[s] - v[nh, s] - h * sum(w[nh, j, s] / factorial(j) for j in 1:nc) for s in 1:ne
+    )
+
+    ExaModels.constraint(
+        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:nh-1, s in 1:ne
+    )
+
+
+
+    ExaModels.constraint(
+        c,
+        Dpp[i,j,1] - u[i,j] * (10.0*pp[i,j,2] - pp[i,j,1]) for i=1:nh, j=1:nc
+    )
+
+    ExaModels.constraint(
+        c, 
+        Dpp[i,j,2] - u[i,j] * (pp[i,j,1] - 10.0*pp[i,j,2]) + (1 - u[i,j])*pp[i,j,2] for i=1:nh, j=1:nc
+    )
+
+
+    ExaModels.constraint(
+        c,
+        v[1, s] - bc for (s, bc) in [(i, bc[i]) for i in 1:ne]
+    )
+
+    return ExaModels.ExaModel(c; kwargs...)
+end
+
+

src/chain.jl

@@ -0,0 +1,72 @@
+# Hanging Chain
+
+# Find the chain (of uniform density) of length L suspended between two points with minimal
+# potential energy.
+
+#   This is problem 4 in the COPS (Version 3) collection of 
+#   E. Dolan and J. More'
+#   see "Benchmarking Optimization Software with COPS"
+#   Argonne National Labs Technical Report ANL/MCS-246 (2004)
+
+function chain_model(n; T = Float64, backend = nothing, kwargs...)
+    nh = max(2, div(n - 4, 4))
+
+    L = 4
+    a = 1
+    b = 3
+    tmin = b > a ? 1 / 4 : 3 / 4
+    tf = 1.0
+    h = tf / nh
+
+    c = ExaModels.ExaCore(T; backend = backend)
+    u = ExaModels.variable(c, nh + 1; start = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1])
+    x1 = ExaModels.variable(c, nh + 1; start = [4 * abs(b - a) * k / nh * (1 / 2 * k / nh - tmin) + a for k in 1:nh+1])
+    x2 = ExaModels.variable(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])
+    x3 = ExaModels.variable(c, nh + 1;  start = [4 * abs(b - a) * (k / nh - tmin) for k in 1:nh+1])
+
+    ExaModels.objective(c, x2[nh + 1])
+
+    ExaModels.constraint(
+        c,
+        x1[j + 1] - x1[j] - 1 / 2 * h * (u[j] + u[j + 1]) for j in 1:nh
+    )
+
+    ExaModels.constraint(
+        c,
+        x1[1] - a
+    )
+
+    ExaModels.constraint(
+        c,
+        x1[nh + 1] - b
+    )
+
+    ExaModels.constraint(
+        c,
+        x2[1]
+    )
+
+    ExaModels.constraint(
+        c,
+        x3[1]
+    )
+
+    ExaModels.constraint(
+        c,
+        x3[nh+1] - L
+    )
+
+    ExaModels.constraint(
+        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
+    )
+
+    ExaModels.constraint(
+        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
+    )
+
+    return ExaModels.ExaModel(c; kwargs...)
+end
+

src/channel.jl

@@ -15,7 +15,14 @@ function channel_model(nh; T = Float64, backend = nothing, kwargs...)
     bc = [0.0 1.0; 0.0 0.0]
     rho = [0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703]
     t = [(i-1)*h for i in 1:nh+1]
-
+    tuples = Vector{Tuple{Int, Int, Int}}()
+    for s in 1:nd
+        for j in s:nd
+            push!(tuples, (s, j, factorial(j-s)))
+        end
+    end    
+    array = [(i,j, factorial(j+nd-i)) for i in 1:nd, j in 1:nc]
+    rho_j = [(j, rho[j]) for j in 1:nc]
     # Initial value
     v0 = zeros(nh, nd)
     for i in 1:nh
@@ -27,55 +34,71 @@ function channel_model(nh; T = Float64, backend = nothing, kwargs...)
 
     core = ExaModels.ExaCore(T; backend= backend)
 
-    v = ExaModels.variable(model, 1:nh, 1:nd)
-    w = ExaModels.variable(model, 1:nh, 1:nc; start=0.0)
+    v = ExaModels.variable(core, nh, nd;)
+    w = ExaModels.variable(core, nh, nc; start=0.0)
 
-    uc = ExaModels.variable(model, 1:nh, 1:nc, 1:nd; start=[v0[i, s] for i=1:nh, j=1:nc, s=1:nd])
-    Duv = ExaModels.variable(model, 1:nh, 1:nc, 1:nd; start=0.0)
+    uc = ExaModels.variable(core, nh, nc, nd; start=[v0[i, s] for i=1:nh, j=1:nc, s=1:nd])
+    Duc = ExaModels.variable(core, nh, nc, nd; start=0.0)
+    y = ExaModels.variable(core, 1; lvar =1, uvar =1, start = 1)
 
     # Constant objective
-    ExaModels.objective(model, Min, 1.0)
+    ExaModels.objective(core, y[1])
 
     # Collocation model
     ExaModels.constraint(
-        model,
-        [i=1:nh, j=1:nc, s=1:nd],
-        uc[i, j, s] == v[i,s] + h*sum(w[i,k]*(rho[j]^k/factorial(k)) for k in 1:nc),
+        core,
+        uc[i, j, s] - v[i,s] - h*sum(w[i,k]*(rho^k/factorial(k)) for k in 1:nc) for i in 1:nh, (j,rho) in rho_j, s in 1:nd
     )
-    ExaModels.constraint(
-        model,
-        [i=1:nh, j=1:nc, s=1:nd],
-        Duc[i, j, s] == sum(v[i,k]*((rho[j]*h)^(k-s)/factorial(k-s)) for k in s:nd) +
-                            h^(nd-s+1) * sum(w[i, k]*(rho[j]^(k+nd-s)/factorial(k+nd-s)) for k in 1:nc)
+    c1 = ExaModels.constraint(
+        core,
+        Duc[i, j, s] for i in 1:nh, j in 1:nc, s in 1:nd
+    )
+
+    ExaModels.constraint!(
+        core,
+        c1,
+        (s-1)*nc*nh + (j-1)*nh + i=> - v[i,k]*((rho*h)^(k-s))/fac_k for i in 1:nh, (j,rho) in rho_j, (s, k, fac_k) in tuples
+    )
+
+    ExaModels.constraint!(
+        core,
+        c1,
+        (s-1)*nc*nh + (j-1)*nh + i => - h^(nd-s+1) * w[i, k]*(rho^(k+nd-s)/fac_nd) for i in 1:nh, (j,rho) in rho_j, (s, k, fac_nd) in array
     )
+
     # Boundary
-    ExaModels.constraint(model, bc_1, v[1, 1] == bc[1, 1])
-    ExaModels.constraint(model, bc_2, v[1, 2] == bc[2, 1])
+    ExaModels.constraint(core, v[1, 1] - bc[1, 1])
+    ExaModels.constraint(core, v[1, 2] -  bc[2, 1])
     ExaModels.constraint(
-        model,
-        bc_3,
+        core,
         sum(v[nh, k]*(h^(k-1)/factorial(k-1)) for k in 1:nd) +
-        h^nd * sum(w[nh, k]/factorial(k+nd-1) for k in 1:nc) == bc[1, 2]
+        h^nd * sum(w[nh, k]/factorial(k+nd-1) for k in 1:nc) - bc[1, 2]
     )
     ExaModels.constraint(
-        model,
-        bc_4,
+        core,
         sum(v[nh, k]*(h^(k-2)/factorial(k-2)) for k in 2:nd) +
-        h^(nd-1) * sum(w[nh, k]/factorial(k+nd-2) for k in 1:nc) == bc[2, 2]
+        h^(nd-1) * sum(w[nh, k]/factorial(k+nd-2) for k in 1:nc) - bc[2, 2]
     )
-    ExaModels.constraint(
-        model,
-        continuity[i=1:nh-1, s=1:nd],
-        sum(v[i, k]*(h^(k-s)/factorial(k-s)) for k in s:nd)
-        + h^(nd-s+1)* sum(w[i, k]/factorial(k+nd-s) for k in 1:nc) == v[i+1, s]
+    c2 = ExaModels.constraint(
+        core,
+        - v[i+1, s] for i in 1:nh-1, s in 1:nd
+    )
+    ExaModels.constraint!(
+        core,
+        c2,
+        (s-1)*(nh-1) + i => v[i, k]*(h^(k-s)/fac_k) for i in 1:nh-1, (s, k, fac_k) in tuples
+    )
+    ExaModels.constraint!(
+        core,
+        c2,
+        (s-1)*(nh-1) + i =>  h^(nd-s+1)* w[i, k]/fac_nd for i in 1:nh-1, (s, k, fac_nd) in array
+
     )
     ExaModels.constraint(
-        model,
-        collocation[i=1:nh, j=1:nc],
-        sum(w[i, k] * (rho[j]^(k-1)/factorial(k-1)) for k in 1:nc) ==
-        R * (Duc[i, j, 2] * Duc[i, j, 3] - Duc[i, j, 1] * Duc[i, j, 4])
+        core,
+        sum(w[i, k] * (rho^(k-1)/factorial(k-1)) for k in 1:nc) -
+        R * (Duc[i, j, 2] * Duc[i, j, 3] - Duc[i, j, 1] * Duc[i, j, 4]) for i in 1:nh, (j,rho) in rho_j
     )
 
     return ExaModels.ExaModel(core; kwargs...)
 end
-