Possibility to restart solving a system starting from a solution (ics override)

[Qfl3x]

Basically the possibility to do solve(discretized_problem, previous_solution, p = new_parameters, otherkwargs...) by overriding the previous ics (or in some other way).

[xtalax]

I'd recommend to do this with an interpolation, There is an idea in the works to to do this automatically in some cases but in needs changes in Interpolations.jl

[Qfl3x]

Any hint how this could be implemented?

[Qfl3x]

remake in SciMLBase supports applying new u0's: https://github.com/SciML/SciMLBase.jl/blob/03cdaedfd3af1ece7c897f6d28ad2d4236ead4d9/src/remake.jl#L52

[xtalax]

Ah ok I have an idea, yes I will implement a helper that will discretize a new u0 on to the grid, or you can use it directly if you have the same discretization

[Qfl3x]

Perhaps something like this (seems to work):

function resolve(ps, sol)
    ps = # Make a vector of pairs here
    _prob = remake(prob, p=ps, u0=vcat([val[end,:] for val in values(sol.u)]...))
    
    solve(_prob, Tsit5(), saveat=0.1)
end

[Qfl3x]

Full MWE (for 2 dependent variables):

using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets

# Parameters, variables, and derivatives
@parameters t x
@variables u(..) v(..)
@parameters α β
Dt = Differential(t)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
eq  = [Dt(u(t, x)) ~ (α + β) * Dxx(u(t, x)),
        Dt(v(t, x)) ~ (α + β) * Dxx(v(t, x))]
bcs = [u(0, x) ~ cos(x),
        u(t, 0) ~ exp(-t),
        u(t, 1) ~ exp(-t) * cos(1),
        v(0,x) ~ sin(x),
        v(t,0) ~ exp(-t),
        v(t,1) ~ 0.,]

# Space and time domains
domains = [t ∈ Interval(0.0, 1.0),
           x ∈ Interval(0.0, 1.0)]

# Parameters
ps = [α => 1.2, β => 2.1]
# PDE system
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x), v(t,x)], ps)

# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x => dx], t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

function pde_solution(ps)
    ps = [α => ps[1], β => ps[2]]
    _prob = remake(prob, p=ps)
    
    solve(_prob, Tsit5(), saveat=0.1)
end

function resolve(ps, sol)
    ps = [α => ps[1], β => ps[2]]
    _prob = remake(prob, p=ps, u0=vcat([val[end,:] for val in values(sol.u)]...))
    
    solve(_prob, Tsit5(), saveat=0.1)
end

sol = pde_solution([10.,10.]);
sol2 = resolve([-1.,-1.], sol);

all(sol2.prob.u0[1:11] .== sol.u[u(t,x)][end,:])