# Recycling solutions of multidimensional NDSolve

Dear wolfram community,

I hope my problem is clear and easy to solve.

I have already solved the following heat equation over a domain:

Clear["Global*"]
Needs["NDSolveFEM"];

pde = 1/r D[r*D[T[t, r, z], {r, 1}], {r, 1}] +
D[D[T[t, r, z], {z, 1}], {z, 1}] - D[T[t, r, z], {t}];

\[CapitalDelta]z = 50*10^(-4);(*[m]*)
R = 5*10^(-2);(*[m]*)
\[Tau] = 3.5*10^(-3)(*[s]*);
Tw = 200;(*[K]*)
\[CapitalOmega] =
ImplicitRegion[0 <= r <= R && -2*\[CapitalDelta]z <= z <= 0, {r, z}];

bc = {DirichletCondition[T[t, r, z] == Tw, z == -2*\[CapitalDelta]z],
DirichletCondition[T[t, r, z] == Tw, r == R]};
(*Note that the unspecified boundaries are set to Neuman zero by \
default*)

ic = {T[0, r, z] == 300 (*[K]*)};

sol = NDSolve[{pde == 0, bc, ic},
T, {t, 0, \[Tau]}, {r, z} \[Element] \[CapitalOmega]]


I optained an Interpolationfunction.

Manipulate[(*Plot mit relativen Farbspectrum*)
Plot3D[T[t, r, z] /. sol, {r, 0, R}, {z, -2*\[CapitalDelta]z, 0},
PlotRange -> {180, Tliq*1.1}, Mesh -> None,
Axes -> True, AxesLabel -> {"r", "z", "T"} ,
AxesStyle -> Thickness[0.001],
ColorFunction -> "TemperatureMap",
(*ColorFunction\[Rule](ColorData["TemperatureMap"][#3]&),*)
ColorFunctionScaling -> True,
ImageSize -> Large,
PlotLegends -> Automatic],
{t, 0, \[Tau], \[Tau]*0.01}]


Now I want to use the temperature distribution T[[Tau],r,z] as part of the initial value (T[0,r,z]) for a new domain (omeganew ) and then solve the heat equation over omeganew. The new domain is essentially the old domain but enlarged by an amount of 4*Δz:

 omeganew = ImplicitRegion[0 <= r <= R && -4*Δz <= z <= 0, {r, z}];


Here is the tricky part: For the new solution I want the Temperature distribution from the previous solution to be "shifted" to the lower part of omeganew and a different function fc[r,z] (lets assume for the purpose of simplicity that fc[r,z]= constant) to be valid within the region that is newly added. So I need a shifted version of the previous solution with respect to the z-coordinate.

In a nutshel, I would like to extract the information (Temperaturevalues at an abitrary time t) obtained by NDsolve (e.g. sol) to construct a Unitstep function (is there a better alternative?) that I would then use as the initial condition for the new domain (omeganew).

• What are the values of tau, R, deltaz, etc? Without them it's not possible to simulate. Please provide them. – user21 Jan 14 at 11:22
• [CapitalDelta]z = 50*10^(-4);([m]) ;[Tau] = 3.5*10^(-2); R = 5*10^(-2); ([m]) – Gustavo Meyagan Jan 14 at 11:39
• ok I have added the additional information. Thanks in advance! – Gustavo Meyagan Jan 14 at 11:48

Here is a way to do it. Let' set up the model:

Δz = 50*10^(-4);(*[m]*)R = 5*10^(-2);(*[m]*)τ =
3.5*10^(-2);(*[s]*)
pde = 1/r D[r*D[T[t, r, z], {r, 1}], {r, 1}] + D[D[T[t, r, z], {z, 1}], {z, 1}] - D[T[t, r, z], {t}];
Ω = ImplicitRegion[0 <= r <= R && -2*Δz <= z <= 0, {r, z}];
bc = {DirichletCondition[T[t, r, z] == Tw, z == -2*Δz], DirichletCondition[T[t, r, z] == Tw, r == R]};


If you now call NDSolveValue you will get a solution (looks like it's zero but that is a different issue)

sol = NDSolveValue[{pde == 0, bc, T[0, r, z] == 0}, T, {t, 0, τ}, {r, z} ∈ Ω];


If you evaluate the inteprolating function out side of the region you will get a warning and an Indetermiante as an answer.

sol[0, -1, 3] Indeterminate


To change that you can use:

sol = NDSolveValue[{pde == 0, bc, T[0, r, z] == 0},
T,
{t, 0, τ},
{r, z} ∈ Ω,
{"ExtrapolationHandler" -> {5 &, "WarningMessage" -> False}}
];


Now you will get the extrapolation value specified (5) and no warning:

sol[0, -1, 3]
5


With this you can then call NDSolveValue on a different domain with a different initial value like so:

sol2 = NDSolveValue[{pde == 0, bc, T[0, r, z] == sol[τ, r, z]},
T,
{t, 0, τ},
{r, z} ∈ ImplicitRegion[ 0 <= r <= R && -10*Δz <= z <= 0, {r, z}]
];

• OHH I also forgot to provide the values for Tw and ic (The "wall temperature" and the initial temperature distribution respectively . They are now included in the post! – Gustavo Meyagan Jan 14 at 12:10
• I want sol[tau, r, z] to be the intial value for the domain1: 0 <= r <= R && -4*Δz <= z <= -Δz (e.g.sol[tau, 0, 0]=T[0,0,-Δz ],sol[tau, 0, -0.5*Δz ]=T[0,0,-0.5*Δz -Δz ] ,and so forth) I want to essentially shift the function downwards. And the initial values of the domain2 0 <= r <= R && -Δz <= z <= 0 to be specified by a completely different function ( lets say a constant) Of course sol2 is to be solved for the union of domain1 and domain2 – Gustavo Meyagan Jan 14 at 12:50
• @GustavoMeyagan, you need to specify everything in the question. Try to do this in a clean, clear way such that people can follow your line of thought. – user21 Jan 14 at 15:14
• @User21I think I've figured out the problem. The interpolation function obtained from NDSolve does not enable you to obtain the value at arbitrary points, for exaple sol[0, -1, 3] will not give the value 5. One actually needs to solve the pde with NDSolveValue instead. From there you can use sol[t,r,z] like any other function including embedding it in a unitstep function.Thank you very much for your help! You guys rock! – Gustavo Meyagan Jan 15 at 10:25
• @GustavoMeyagan, good that you figured it out. NDSolve and NDSolveValue essentially return the same result. NDSolve will give you a Rule where the second argument is the interpolating function. To use the result from NDSolve you first need to replace the dependent variable to get to the interpolating function. With NDSolveValue` there is not such need as it directly returns the Interpolating function. – user21 Jan 15 at 11:19