Skip to content

[COPS] add PDE instances #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
frapac wants to merge 3 commits into
base: main
Choose a base branch
from
Draft

[COPS] add PDE instances #3

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
7156779
[COPS] add PDE instances
frapac May 21, 2024
40f5869
add missing files
frapac May 21, 2024
91b687b
Fixed pde_models and added models (#7)
jerryfung110 Nov 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion src/ExaModelsExamples.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,13 +15,14 @@ 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")

const NAMES = filter(names(ExaModelsExamples; all = true)) do x
str = string(x)
Expand Down
230 changes: 230 additions & 0 deletions src/pde_models.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,230 @@

include("shapes/circle.jl")
Copy link
Copy Markdown
Member

@sshin23 sshin23 May 21, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

these files are not committed

include("shapes/circle_rec.jl")
include("shapes/rectangle.jl")

struct PDEProblem
a::Float64
b::Vector{Float64}
c::Vector{Float64}
d::Vector{Float64}
p::Vector{Float64}
end

struct PDEDiscretizationDomain
NODES::Int
ELEM::Int
DIMEN::Int
BREAK::Int
AREA::Vector{Float64}
TRIANG::Matrix{Int}
COORDS::Matrix{Float64}
BNDRY::Vector{Int}
EDGE::Array{Float64, 3}
US::Vector{Float64}
UE::Vector{Float64}
end

function PDEDiscretizationDomain(nh, domain::Dict)
NODES = domain[:NODES]
ELEMS = domain[:ELEMS]
COORDS = domain[:COORDS]
ELEMS = domain[:ELEMS]
BNDRY = domain[:BNDRY]
DIMEN = 2
BREAK = nh

# Description of triangular elements
TRIANG = domain[:TRIANG]
# Edge lengths
EDGE = [
COORDS[TRIANG[e, mod(d1, DIMEN+1)+1], d2] - COORDS[TRIANG[e, d1], d2]
for e in 1:ELEMS, d1 in 1:DIMEN+1, d2 in 1:DIMEN
]
# Area of element
AREA = [(EDGE[e, 1, 1]*EDGE[e, 2, 2] - EDGE[e, 1, 2]*EDGE[e, 2, 1]) / 2.0 for e in 1:ELEMS]
US = domain[:US] # starting point
UE = domain[:UE] # ending point

return PDEDiscretizationDomain(
NODES, ELEMS, DIMEN, BREAK,
AREA, TRIANG, COORDS, BNDRY, EDGE, US, UE,
)
end

function _update_values!(integral, energy, u, problem::PDEProblem, dom::PDEDiscretizationDomain)
# Unpack values
a, b, c, d, p = problem.a, problem.b, problem.c, problem.d, problem.p
TRIANG, EDGE = dom.TRIANG, dom.EDGE

for b1 in 1:dom.BREAK+2, e1 in 1:dom.ELEM
integral[b1, e1] =
dom.AREA[e1]*(
1 / (dom.DIMEN+1) *
(sum((b[TRIANG[e1,c1]]*u[b1,TRIANG[e1,c1]]^2/2-
c[TRIANG[e1,c1]]*u[b1,TRIANG[e1,c1]]^(p[TRIANG[e1,c1]]+1)/(p[TRIANG[e1,c1]]+1)+
d[TRIANG[e1,c1]]*u[b1,TRIANG[e1,c1]]) for c1 in 1:dom.DIMEN+1)) +
a / (8*dom.AREA[e1]^2)*(
u[b1,TRIANG[e1,1]]^2*(EDGE[e1,2,1]^2 + EDGE[e1,2,2]^2) +
u[b1,TRIANG[e1,2]]^2*(EDGE[e1,3,1]^2 + EDGE[e1,3,2]^2) +
u[b1,TRIANG[e1,3]]^2*(EDGE[e1,1,1]^2 + EDGE[e1,1,2]^2) +
2*u[b1,TRIANG[e1,1]]*u[b1,TRIANG[e1,2]]*(EDGE[e1,2,1]*EDGE[e1,3,1] + EDGE[e1,2,2]*EDGE[e1,3,2]) +
2*u[b1,TRIANG[e1,1]]*u[b1,TRIANG[e1,3]]*(EDGE[e1,2,1]*EDGE[e1,1,1] + EDGE[e1,2,2]*EDGE[e1,1,2]) +
2*u[b1,TRIANG[e1,2]]*u[b1,TRIANG[e1,3]]*(EDGE[e1,1,1]*EDGE[e1,3,1] + EDGE[e1,1,2]*EDGE[e1,3,2])
)
)
end
for b1 in 1:dom.BREAK+2
energy[b1] = sum(integral[b1, e1] for e1 in 1:dom.ELEM)
end
end

function _initial_position!(problem::PDEProblem, d::PDEDiscretizationDomain, niter)
ubar = zeros(d.NODES)
u0 = zeros(d.BREAK+2, d.NODES)
energy0 = zeros(d.BREAK+2)
integral0 = zeros(d.BREAK+2, d.ELEM)

# Calculate an ending point on the other side of the barrier
for i in 1:niter
for b1 in 1:d.BREAK+2, n in 1:d.NODES
u0[b1, n] = (1 - (b1-1)/(d.BREAK+1))*d.US[n] + ((b1-1)/(d.BREAK+1))*d.UE[n]
end
_update_values!(integral0, energy0, u0, problem, d)
if energy0[end] < 0
break
end
d.UE .*= 2.0
end

# Backtrack to the barrier to get a better representation
ubar .= d.US
for i in 1:niter
for b1 in 1:d.BREAK+2, n in 1:d.NODES
u0[b1, n] = (1 - (b1-1)/(d.BREAK+1))*ubar[n] + ((b1-1)/(d.BREAK+1))*d.UE[n]
end
_update_values!(integral0, energy0, u0, problem, d)
ebar = maximum(energy0[b1] for b1 in 1:d.BREAK+2 if energy0[b1] < 0.0)
bbar = minimum(b1 for b1 in 1:d.BREAK+2 if energy0[b1] == ebar)
ubar .= u0[bbar-1, :]
d.UE .= u0[bbar, :]
end

for b1 in 1:d.BREAK+2, n in 1:d.NODES
u0[b1, n] = (1 - (b1-1)/(d.BREAK+1))*d.US[n] + ((b1-1)/(d.BREAK+1))*d.UE[n]
end
_update_values!(integral0, energy0, u0, problem, d)
z0 = maximum(energy0[b1] for b1 in 1:d.BREAK+2)
return (
u=u0, energy=energy0, integral=integral0, z=z0,
)
end

function _transition_state_model(problem, dom::PDEDiscretizationDomain; T = Float64, backend = nothing, kwargs...)
a, b, c, d, p = problem.a, problem.b, problem.c, problem.d, problem.p
x0 = _initial_position!(problem, dom, 10)
_proto_model()

ALPHA = 2.0
H = ALPHA / (dom.BREAK+1) * sqrt(sum((dom.US[n] - dom.UE[n])^2 for n in 1:dom.NODES))

# Build optimization problem
core = ExaModels.ExaCore(T; backend= backend)

u = ExaModels.variable(core, 1:dom.BREAK+2, 1:dom.NODES; start=x0.u)
integral = ExaModels.variable(core, 1:dom.BREAK+2, 1:dom.ELEM)
z = ExaModels.variable(core, 1; start=x0.z)

ExaModels.objective(core, z[1])

ExaModels.constraint(
core,
sum(integral[b1+1, e1] for e1 in 1:n) - z[1]
for b1 in 2:dom.ELEMS+1;
lcon=-Inf,
ucon=0.0,
)
Copy link
Copy Markdown
Member

@sshin23 sshin23 May 21, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These constraints can benefit significantly from constraint! if n is large.

Consider reformulating this as

c = ExaModels.constraint(
    core, 
    - z[1] for b1 in 2:dom.ELEMS+1;
    lcon=-Inf,
    ucon=0.0,
) 
ExaModels.constraint!(
    core,
    b1 -1 => integral[b1+1, e1] for b1 in 2:dom.ELEMS+1
)

ExaModels.constraint(
core,
sum((u[b1+1, n] - u[b1, n])^2 for n in 1:dom.NODES)
for b1 in 1:dom.BREAK+1;
lcon=-Inf,
ucon=H^2,
)
ExaModels.constraint(
core,
dom.AREA[e1]*(
1 / (dom.DIMEN+1) *
(sum((b[dom.TRIANG[e1,c1]]*u[b1,dom.TRIANG[e1,c1]]^2/2-
c[dom.TRIANG[e1,c1]]*u[b1,dom.TRIANG[e1,c1]]^(p[dom.TRIANG[e1,c1]]+1)/(p[dom.TRIANG[e1,c1]]+1)+
d[dom.TRIANG[e1,c1]]*u[b1,dom.TRIANG[e1,c1]]) for c1 in 1:dom.DIMEN+1)) +
a / (8*dom.AREA[e1]^2)*(
u[b1,dom.TRIANG[e1,1]]^2*(dom.EDGE[e1,2,1]^2 + dom.EDGE[e1,2,2]^2) +
u[b1,dom.TRIANG[e1,2]]^2*(dom.EDGE[e1,3,1]^2 + dom.EDGE[e1,3,2]^2) +
u[b1,dom.TRIANG[e1,3]]^2*(dom.EDGE[e1,1,1]^2 + dom.EDGE[e1,1,2]^2) +
2*u[b1,dom.TRIANG[e1,1]]*u[b1,dom.TRIANG[e1,2]]*(dom.EDGE[e1,2,1]*dom.EDGE[e1,3,1] + dom.EDGE[e1,2,2]*dom.EDGE[e1,3,2]) +
2*u[b1,dom.TRIANG[e1,1]]*u[b1,dom.TRIANG[e1,3]]*(dom.EDGE[e1,2,1]*dom.EDGE[e1,1,1] + dom.EDGE[e1,2,2]*dom.EDGE[e1,1,2]) +
2*u[b1,dom.TRIANG[e1,2]]*u[b1,dom.TRIANG[e1,3]]*(dom.EDGE[e1,1,1]*dom.EDGE[e1,3,1] + dom.EDGE[e1,1,2]*dom.EDGE[e1,3,2])
)
) - integral[b1, e1]
for b1 in 1:dom.BREAK+2, e1 in 1:dom.ELEM
)

# Boundary
boundary_nodes = findall(isequal(1), dom.BNDRY)
ExaModels.constraint(
core,
u[b1+1, n] for b1 in 1:dom.BREAK, n in boundary_nodes
)
ExaModels.constraint(
core,
u[1, n] for n in 1:dom.NODES;
lcon=dom.US,
ucon=dom.US,
)
ExaModels.constraint(
core,
u[dom.BREAK+2, n] for n in 1:dom.NODES
lcon=dom.UE,
ucon=dom.UE,
)
return core
end

function dirichlet_model(nh)
dom = PDEDiscretizationDomain(nh, CIRCLE_DOMAIN)
println(dom.ELEM)
pb = PDEProblem(
0.01,
fill(1.0, dom.NODES),
fill(1.0, dom.NODES),
fill(0.0, dom.NODES),
fill(3.0, dom.NODES),
)
return _transition_state_model(pb, dom)
end

function henon_model(nh)
dom = PDEDiscretizationDomain(nh, CIRCLE_REC_DOMAIN)
pb = PDEProblem(
1.0,
fill(0.0, dom.NODES),
sqrt.(dom.COORDS[:, 1].^2 .+ dom.COORDS[:, 2].^2),
fill(0.0, dom.NODES),
fill(3.0, dom.NODES),
)
return _transition_state_model(pb, dom)
end

function lane_emden_model(nh)
dom = PDEDiscretizationDomain(nh, RECTANGLE_DOMAIN)
pb = PDEProblem(
1.0,
fill(0.0, dom.NODES),
fill(1.0, dom.NODES),
fill(0.0, dom.NODES),
fill(3.0, dom.NODES),
)
return _transition_state_model(pb, dom)
end