Initialize some functions at every remake in PDEProblem

I looked into this earlier this year, and while it is a known missing feature, it seems to still be in progress.

While your problem may be different than mine, since I only ever needed on the order of 10 boundary condition sets at a time, my solution was to initialize all runs with different datasets at once in a multithreaded loop with the interpolation defined within the loop and to save the discretized problems in a vector for easy access. I could then modify non-BC parameters across all problems without re-discretizing. Some generic code similar to what I used is below and it works well enough, though it’s not very elegant.

#define some input_data matrix or data structure such that it can be indexed for each BC set.
probvec = []
solvec = []
Threads.@threads for i = 1:n_sims
     #just need to somehow index with loop such that t_input, x_input change with i
     t_input = input_data[i,1,:] 
     x_input = input_data[i,2,:]
     x_BC = LinearInterpolation(x_input,t_input)

     #rest of code
     #relevant BC
     x(t,0) ~ x_BC(t)
     

     #more code, discretization
      prob = ModelingToolkit.discretize(pdesys,discretization);
     push!(probvec, (i, prob)) #saved as tuple for indexing since multithreading means code can finish in random orders


     sol = solve(prob) #done to run compilation and speed up future solves with different non-BC parameters
     push!(solvec, (i, sol)) #saved to check whether solve ran properly, but not really necessary
end

#strip probs out of tuple and put in correct order
#definitely could be done more elegantly, but I really just pretend to be a programmer
ordered_prob_tuples = sort(probvec, rev=false)
ordered_probs = []
for i = 1:size(ordered_prob_tuples,1) 
    push!(ordered_probs, ordered_prob_tuples[i][2])
end

#remake for non-BC parameters is fast, uses i to find correct problem 
#from problem vector and remake with some set of parameters from param_vector
newprob = remake(ordered_probs[i], p = param_vector)
sol = solve(newprob)

I still hope to figure out a better way to pass in BCs since I will eventually want to modify those parametrically, but I haven’t had much luck yet beyond some extremely low-quality experiments with series approximations. If you end up figuring something better out, please post about it, I’d be very interested to hear!