Creating PDEs from ODEs

Modeling the Brusselator PDE:

\[\begin{split} \begin{align} \frac{\partial u}{\partial t} &= 1 + u^2v - 4.4u + \alpha (\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}) + f(x, y, t) \\ \frac{\partial v}{\partial t} &= 3.4u - u^2 v + \alpha (\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}) \end{align} \end{split}\]

import ComponentArrays as CA
using SimpleUnPack
using OrdinaryDiffEq
using LinearSolve
using Plots

function build_brussilator(; N::Int=32, x_min=0, y_min=0, t_min=0, x_max=1, y_max=1, t_max=11.5, α = 10.0)
    @assert (x_min < x_max) && (y_min < y_max) && (t_min < t_max) && (N > 1)
    u0 = CA.ComponentArray(u = zeros(N, N), v = zeros(N, N))
    xx = range(x_min, x_max, length=N)
    yy = range(y_min, y_max, length=N)
    dx = step(xx)
    dy = step(yy)
    ps = CA.ComponentArray(; α, xx, yy, N, dx, dy)

    for i in 1:N, j in 1:N
        x = xx[j]
        y = yy[i]
        u0.u[i, j] = 22 * (y * (1 - y))
        u0.v[i, j] = 27 * (x * (1 - x))
    end

    odef = (ds, s, p, t) -> begin
        @unpack α, xx, yy, dx, dy = p
        @unpack u, v = s
        brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5
        for (i, y) in enumerate(yy), (j, x) in enumerate(xx)
            uterm = 1.0 + v[i, j] * u[i, j]^2 - 4.4 * u[i, j] + brusselator_f(x, y, t)
            vterm = 3.4 * u[i, j] - v[i, j] * u[i, j]^2
            # No flux boundary conditions
            i_prev = clamp(i - 1, 1, N)
            i_next = clamp(i + 1, 1, N)
            j_prev = clamp(j - 1, 1, N)
            j_next = clamp(j + 1, 1, N)
            unabla = (u[i, j_prev] - 2u[i, j] + u[i, j_next]) / dx^2 + (u[i_prev, j] - 2u[i, j] + u[i_next, j]) / dy^2
            vnabla = (v[i, j_prev] - 2v[i, j] + v[i, j_next]) / dx^2 + (v[i_prev, j] - 2v[i, j] + v[i_next, j]) / dy^2
            ds.u[i, j] = α * unabla + uterm
            ds.v[i, j] = α * vnabla + vterm
        end
    end
    return (; u0, odef, ps)
end

build_brussilator (generic function with 1 method)

@unpack u0, odef, ps = build_brussilator()
tspan = (0.0, 11.5)
prob = ODEProblem(odef, u0, tspan, ps)
@time sol = solve(prob, KenCarp47(linsolve = KrylovJL_GMRES()), saveat=0.1)

 36.019777 seconds (1.25 G allocations: 33.748 GiB, 10.73% gc time, 10.77% compilation time)

retcode: Success
Interpolation: 1st order linear
t: 116-element Vector{Float64}:
  0.0
  0.1
  0.2
  0.3
  0.4
  0.5
  0.6
  0.7
  0.8
  0.9
  ⋮
 10.7
 10.8
 10.9
 11.0
 11.1
 11.2
 11.3
 11.4
 11.5
u: 116-element Vector{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(u = ViewAxis(1:1024, ShapedAxis((32, 32))), v = ViewAxis(1025:2048, ShapedAxis((32, 32))))}}}}:
 ComponentVector{Float64}(u = [0.0 0.0 … 0.0 0.0; 0.6867845993756503 0.6867845993756503 … 0.6867845993756503 0.6867845993756503; … ; 0.6867845993756498 0.6867845993756498 … 0.6867845993756498 0.6867845993756498; 0.0 0.0 … 0.0 0.0], v = [0.0 0.8428720083246618 … 0.8428720083246611 0.0; 0.0 0.8428720083246618 … 0.8428720083246611 0.0; … ; 0.0 0.8428720083246618 … 0.8428720083246611 0.0; 0.0 0.8428720083246618 … 0.8428720083246611 0.0])
 ComponentVector{Float64}(u = [6.773290649233472 6.773290725936031 … 6.773290725936014 6.773290649233484; 6.77328532685449 6.773285414416952 … 6.773285414416986 6.773285326854452; … ; 6.773285326854472 6.773285414416982 … 6.773285414416978 6.773285326854461; 6.773290649233485 6.7732907259360235 … 6.773290725936022 6.773290649233483], v = [0.6709604463847637 0.6709539992600472 … 0.6709539992600495 0.6709604463847645; 0.6709605815409581 0.6709541235532056 … 0.6709541235532068 0.670960581540964; … ; 0.6709605815409634 0.6709541235531926 … 0.6709541235531861 0.670960581540975; 0.6709604463847623 0.6709539992600636 … 0.6709539992600494 0.670960446384754])
 ComponentVector{Float64}(u = [6.356141376219994 6.35614137371218 … 6.3561413737121795 6.3561413762199885; 6.356141371068119 6.356141368581267 … 6.356141368581257 6.35614137106812; … ; 6.3561413710681185 6.356141368581258 … 6.3561413685812695 6.356141371068112; 6.3561413762199885 6.356141373712178 … 6.35614137371218 6.356141376219993], v = [0.5264757088675984 0.526475702969669 … 0.5264757029696688 0.5264757088675975; 0.5264757071909414 0.5264757012721475 … 0.5264757012721479 0.5264757071909418; … ; 0.5264757071909406 0.5264757012721485 … 0.5264757012721483 0.5264757071909414; 0.5264757088675993 0.5264757029696683 … 0.5264757029696681 0.5264757088675988])
 ComponentVector{Float64}(u = [5.8039450312330025 5.803945032670937 … 5.803945032670955 5.8039450312330105; 5.803945008829085 5.803945010331297 … 5.803945010331298 5.803945008829071; … ; 5.803945008829085 5.803945010331292 … 5.803945010331283 5.803945008829103; 5.803945031233014 5.803945032670929 … 5.803945032670947 5.80394503123299], v = [0.5710329001041721 0.5710328713798584 … 0.5710328713798626 0.5710329001041705; 0.5710329002710766 0.5710328714823919 … 0.5710328714823888 0.5710329002710778; … ; 0.5710329002710796 0.5710328714823935 … 0.5710328714823905 0.571032900271078; 0.5710329001041679 0.5710328713798644 … 0.5710328713798642 0.571032900104169])
 ComponentVector{Float64}(u = [5.297578205206687 5.297578209392803 … 5.297578209392773 5.297578205206653; 5.297578175531061 5.297578179783483 … 5.297578179783501 5.297578175531086; … ; 5.297578175531074 5.297578179783508 … 5.297578179783506 5.29757817553104; 5.29757820520665 5.297578209392824 … 5.297578209392806 5.2975782052067], v = [0.6226807503070436 0.6226807123779619 … 0.6226807123779533 0.622680750307051; 0.6226807524560014 0.6226807144605842 … 0.6226807144605878 0.6226807524559981; … ; 0.6226807524559957 0.6226807144605788 … 0.6226807144605859 0.6226807524559975; 0.622680750307053 0.6226807123779476 … 0.62268071237795 0.6226807503070518])
 ComponentVector{Float64}(u = [4.835550075851494 4.835550077670883 … 4.835550077670894 4.835550075851434; 4.835550062182636 4.835550064035003 … 4.835550064035016 4.835550062182632; … ; 4.835550062182667 4.83555006403501 … 4.835550064034994 4.8355500621826435; 4.835550075851446 4.835550077670894 … 4.835550077670941 4.835550075851473], v = [0.6783892452223519 0.6783892276259822 … 0.6783892276259791 0.67838924522236; 0.6783892460201488 0.6783892283908307 … 0.6783892283908279 0.6783892460201475; … ; 0.6783892460201474 0.678389228390825 … 0.6783892283908277 0.6783892460201472; 0.6783892452223589 0.678389227625972 … 0.6783892276259778 0.6783892452223582])
 ComponentVector{Float64}(u = [4.413411517650755 4.413411521170437 … 4.413412575364365 4.413412580913129; 4.413411503749828 4.413411507213682 … 4.413412571571171 4.413412577172818; … ; 4.413410275050868 4.413410275384847 … 4.413412128038978 4.413412137957223; 4.4134102875010655 4.41341028831333 … 4.413412127953709 4.413412137828075], v = [0.7384241955416518 0.7384241947306054 … 0.7384241563610594 0.7384241568103019; 0.7384241958614747 0.7384241950496698 … 0.7384241565335218 0.738424156982771; … ; 0.7384242283224687 0.7384242274903339 … 0.7384241743468276 0.7384241747148992; 0.7384242282496445 0.7384242274163991 … 0.7384241743429102 0.7384241747104546])
 ComponentVector{Float64}(u = [4.0252218880435775 4.025222368295847 … 4.025369779908983 4.025370567726263; 4.025220010016473 4.025220482358175 … 4.02536931517563 4.025370110494713; … ; 4.025048196758811 4.025048231437623 … 4.025307294522808 4.025308693447834; 4.025049871923686 4.0250499735892875 … 4.025307216918114 4.025308609598616], v = [0.8033001849117587 0.8033001538675842 … 0.8032947885203122 0.8032947689722305; 0.803300221364924 0.8033001903211177 … 0.8032948044737356 0.803294784820357; … ; 0.8033047604989224 0.8033047266120706 … 0.80329729537036 0.8032972643655582; 0.8033047585844018 0.8033047244358746 … 0.8032973029851302 0.8032972720128221])
 ComponentVector{Float64}(u = [3.667604663956293 3.667606079043101 … 3.6680410197742037 3.6680433461501067; 3.6675991332963944 3.6676005250291137 … 3.668039659050114 3.668042007574265; … ; 3.6670921950549653 3.667092295454984 … 3.6678566659501826 3.6678607954252174; 3.667097127165844 3.6670974252256965 … 3.6678564264958817 3.667860537532007], v = [0.8733020382987315 0.8733019596132584 … 0.8732861290554518 0.8732860584675328; 0.8733021443355864 0.8733020656652141 … 0.8732861746214222 0.8732861037088211; … ; 0.8733155371355844 0.8733154500767333 … 0.8732935240586448 0.8732934196534524; 0.873315533005648 0.873315445160705 … 0.8732935480310696 0.8732934437357457])
 ComponentVector{Float64}(u = [3.3392384252934137 3.3392407389810654 … 3.339952212316578 3.3399560188901236; 3.3392293843418255 3.3392316598191023 … 3.339949992512025 3.339953835323569; … ; 3.338400137926853 3.338400301050081 … 3.339650653540949 3.3396574096209806; 3.338408199775564 3.3384086862375626 … 3.3396502557837895 3.3396569816933726], v = [0.9483680454734608 0.9483679242184885 … 0.9483420286907633 0.9483419057652567; 0.9483682180440728 0.9483680968218225 … 0.9483421023514855 0.9483419788868531; … ; 0.948390125902637 0.9483899909585555 … 0.9483541245150184 0.9483539462634124; 0.948390120030725 0.9483899837927395 … 0.9483541646047681 0.9483539865408929])
 ⋮
 ComponentVector{Float64}(u = [0.4168802797269945 0.4168741647838941 … 0.4150735538034345 0.4150639740405664; 0.4169015964043824 0.4168955656681689 … 0.41507907083134615 0.4150694069099402; … ; 0.41891794665429755 0.41891642859239336 … 0.415848254909682 0.4158317324390199; 0.4189003028908161 0.4188980638980493 … 0.4158495479781351 0.415833093326587], v = [4.146044990330784 4.146044992898402 … 4.146045615883768 4.146045618971389; 4.146044982735447 4.146044985307038 … 4.146045614175036 4.146045617291338; … ; 4.146044092179732 4.146044096180869 … 4.146045381686316 4.1460453880296; 4.146044094218854 4.146044098298231 … 4.146045381305585 4.146045387641745])
 ComponentVector{Float64}(u = [0.42160485892286464 0.4215987329878546 … 0.41979560791035864 0.41978601533623144; 0.4216261908855238 0.4216201491611287 … 0.41980113189461304 0.41979145511566024; … ; 0.4236445658251053 0.4236430356107171 … 0.4205715064970574 0.42055496655345104; 0.4236269225509242 0.42362467129595244 … 0.42057280391971247 0.42055633181195284], v = [4.214532470506935 4.214532484230108 … 4.214535655800392 4.214535671869559; 4.2145324475016235 4.214532461225123 … 4.21453564704525 4.214535663189571; … ; 4.214529510182169 4.214529526515946 … 4.214534206536349 4.2145342305473665; 4.214529511648382 4.214529528170273 … 4.214534201722006 4.214534225709187])
 ComponentVector{Float64}(u = [0.4264496863106328 0.4264435428940228 … 0.424636496230486 0.42462688379550795; 0.4264710400333443 0.426464980827362 … 0.42464203090367286 0.424632334202573; … ; 0.428492458606303 0.42849090878072355 … 0.42541428159958483 0.4253977153299919; 0.42847481844080143 0.42847254744139623 … 0.42541558624726117 0.42539908783114905], v = [4.282093935904858 4.282093967368634 … 4.282101114976782 4.282101151175451; 4.2820938909447115 4.282093922408954 … 4.282101095391377 4.282101131726782; … ; 4.282087874924879 4.2820879111580705 … 4.282097751262803 4.282097801908266; 4.2820878731440555 4.282087909698825 … 4.282097739097958 4.282097789702016])
 ComponentVector{Float64}(u = [0.4317006854308529 0.43169452369264544 … 0.4298833656619009 0.4298737323956933; 0.4317220620331911 0.43171598450499715 … 0.42988891154332814 0.42987919394638513; … ; 0.43374668104447744 0.4337451105880318 … 0.4306631262475786 0.4306465323035243; 0.433729044198553 0.43372675242960046 … 0.43066443844082997 0.4306479123683452], v = [4.348550802198796 4.348550852248083 … 4.348562166761299 4.34856222406542; 4.348550734158045 4.348550784208577 … 4.348562135823456 4.348562193329223; … ; 4.348541481677921 4.348541538837029 … 4.348556799680815 4.348556878317282; 4.348541476436427 4.34854153405603 … 4.348556779843163 4.348556858420203])
 ComponentVector{Float64}(u = [0.4378002584832121 0.4377940813695529 … 0.43597943640557013 0.435969785400775; 0.43782165554957025 0.4378155626470363 … 0.4359849918703541 0.4359752564744689; … ; 0.4398490565950912 0.43984746882656356 … 0.436760858923081 0.43674424096258885; 0.4398314215909431 0.43982911237665134 … 0.4367621772706929 0.43674562720033683], v = [4.413591250545334 4.413591316179601 … 4.413606161870376 4.413606237131992; 4.413591161875608 4.413591227510145 … 4.413606121233814 4.413606196757692; … ; 4.413579098331624 4.413579173037963 … 4.4135991105046575 4.413599213412576; 4.413579091138869 4.413579166439169 … 4.413599084412596 4.413599187242562])
 ComponentVector{Float64}(u = [0.44481239406421086 0.4448062029367525 … 0.44298835263967795 0.4429786852761522; 0.44483381071913436 0.4448277038049812 … 0.44299391697077795 0.44298416515660377; … ; 0.4468638189252345 0.44686221537514803 … 0.4437712930410321 0.4437546526017796; 0.44684618490592026 0.4468438597762696 … 0.44377261686533337 0.4437560443356255], v = [4.4770558806525065 4.477055960481953 … 4.477074049929438 4.477074141740255; 4.477055772247006 4.477055852074668 … 4.477074000325992 4.477074092458875; … ; 4.477041076054254 4.477041166746501 … 4.477065461328566 4.47706558693906; 4.477041067783887 4.477041159204095 … 4.477065429674004 4.477065555187166])
 ComponentVector{Float64}(u = [0.4528010807703504 0.4527948754062516 … 0.4509737585424087 0.4509640744872246; 0.45282251767853643 0.45281639653033395 … 0.45097933193062717 0.4509695633625692; … ; 0.454855201701998 0.4548535821126145 … 0.45175824201637216 0.4517415785156908; 0.4548375684314332 0.45483522712217656 … 0.4517595713532354 0.4517429757827135], v = [4.538785292228352 4.5387853864739025 … 4.538806780563911 4.538806889255296; 4.5387851634218475 4.538785257663148 … 4.538806721803528 4.538806830879466; … ; 4.538767766014077 4.538767872948275 … 4.538796629746775 4.538796778641648; 4.538767756903987 4.538767864706658 … 4.538796592495138 4.538796741272092])
 ComponentVector{Float64}(u = [0.4618303071981311 0.4618240857900572 … 0.45999929829194625 0.4599895954993913; 0.4618517665644288 0.4618456293745959 … 0.46000488183592864 0.4599950944621444; … ; 0.4638874385924718 0.4638858009177917 … 0.46078551926404115 0.46076882999891616; 0.46386980645543036 0.4638674469080928 … 0.4607868548628823 0.46077023355010216], v = [4.598620084980905 4.598620195474215 … 4.598645303399222 4.598645431042206; 4.598619933549734 4.598620044036595 … 4.598645234369958 4.598645362466149; … ; 4.5985995193793565 4.598599644628915 … 4.598633393353519 4.598633568265254; 4.598599509031673 4.598599635302711 … 4.598633349743743 4.59863352451541])
 ComponentVector{Float64}(u = [0.4719640619440537 0.4719578211001766 … 0.47012861606647505 0.47011889077805186; 0.47198554751346433 0.4719793908892691 … 0.4701342117727099 0.4701244018251088; … ; 0.4740247632637462 0.47402310366950895 … 0.4709169381989794 0.47090021834605006; 0.47400713326586036 0.4740047516277393 … 0.47091828152275794 0.4709016296462936], v = [4.656400858618204 4.656400988801655 … 4.656430568060796 4.65643071846607; 4.656400680780273 4.656400810956011 … 4.656430486728823 4.656430637665617; … ; 4.656376687318359 4.6563768347740515 … 4.656416529743033 4.656416735554792; 4.656376674699452 4.656376823348102 … 4.656416478287569 4.656416683935198])

---

This notebook was generated using Literate.jl.