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!