3
$\begingroup$

I'm trying to solve the heat equation in 3D using FEM and back out a slice in the xz plane. My problem is when I define a function that is a slice in the xz plane it uses the full 3D interpolation function instead of creating a new 2D interpolation function. This is a problem because I plan to integrate another formula using the spatial/temporal temperature results on that slice. If it keeps calling the 3D interpolation function then the script I am writing will take forever. There has to be a way to do it in Mathematica but sadly I am stuck.

The setup is: -A rectangular body that is composed of two different materials (one on the bottom and one on the top in the z direction).
-All walls on the sides (xy planes) are insulated. -The bottom is at a constant temperature (z=0). -The top has convection (z=zThkTot). -A forcing function is used to sweep in the x direction.

I can solve the problem with ease and view the results which are correct for my needs. I've cleaned up the code so anyone can copy and paste it. The main part I am struggling with is the:

TempSol2D[x_, z_, t_] = TempSol3D[x, ywidth/2, z, t]

The interpolation function is using the mesh to obtain values. Some of the objects I have been playing with are:

TempSol3D["Methods"]
TempSol3D["Properties"]
TempSol3Dvalues = TempSol3D["ValuesOnGrid"];
TempSol3Dcoord = TempSol3D["Coordinates"]
TempSol3D[[3]][[1]]["Properties"]
TempSol3D[[3]][[1]]["Coordinates"]

Which gives good insight to how it is constructed such as the time array and grid coordinates. Any help would be appreciated.

Here is the all the code.

Needs["NDSolve`FEM`"]
Clear[Derivative]
Clear["Global`*"]

BoolMeshRefine = False; (*Refine the mesh?*)

Tair = 23; (*Fluid temperature of the air*)
Tinit2 = 23; (*Region 1 initial temperature*)
Tinit1 = 80; (*Region 2 initial temperature*)
Tbottom = 80; (*Temperature at the end of z*)
k2 = 0.100; (*Region 1 conductivity*)
rho2 = 1500; (*Region 1 density*)
Cp2 = 1200; (*Region 1 heat capacity*)
k1 = 2.0; (*Region 2 conductivity*)
rho1 = 2200; (*Region 2 density*)
Cp1 = 700;(*Region 2 heat capacity*)
hconv = 6.5; (*Convection coeff*)

xend = 0.005; (*Total length of mesh*)
ywidth = 0.002; (*Width of mesh*)
zThk1 = 0.001;(*Thickness of material 1*)
zThk2 = 0.001;(*Thickness of material 2*)

tend = 5; (*End time for simulation*)

zThkTot = zThk1 + zThk2; (*Total thickness of the materials*)

ForcingVelx = 0.0025; (*Velocity in the x dir m/s*)
ForcingAmp = 1*10^9; (*Power density in W/m^3*)

(*TEMPERATURE solving variables*)
TempMaxStepSize = 0.50; (*in seconds*)

(*MESH variables*)
MeshCellSize = 0.00020; (*Default cell size*)
MeshCellRefSize = 0.00005; (*Refinement size*)
MeshRefDepth = 0.00025; (*Depth from the top of the material 1*)
MeshRefInterfaceWidth = 0.0005; (*Width for refinement where material \
1 meets material 2*)

(*Generate MESH=========================================*)
MeshRefine0 = 
 Function[{vertices, area}, 
  Block[{x, y, z}, {x, y, z} = Mean[vertices];
   If[(z > (zThkTot - MeshRefDepth)) || (zThk1 - 
        MeshRefInterfaceWidth/2 < z < 
       zThk1 + MeshRefInterfaceWidth/2), area > (MeshCellRefSize)^3, 
    area > (MeshCellSize)^3]]];
(*Create the mesh*)
If[BoolMeshRefine == True,
  Mesh0 = 
    ToElementMesh[Cuboid[{0, 0, 0}, {xend, ywidth, zThkTot}], 
     MaxCellMeasure -> MeshCellSize^3, 
     MeshRefinementFunction -> MeshRefine0];,
  Mesh0 = 
    ToElementMesh[Cuboid[{0, 0, 0}, {xend, ywidth, zThkTot}], 
     MaxCellMeasure -> MeshCellSize^3];
  ];

Show[Mesh0["Wireframe"]]

(*Forcing function to apply heat*)
Forcing[x_, y_, z_, t_] = 
  ForcingAmp*Sin[((ForcingVelx*t - x)/xend)*2*Pi];

(*Create functions to give different properties to the different \
materials*)
v1 = rho1*Cp1 (*Combination of density and heat capacity*);
v2 = rho2*Cp2;
C1[z_] := 
  Evaluate[With[{p = zThk1, k1 = k1, k2 = k2}, 
    If[z < p, k1, k2]]] (*Variable conductivity*);
M1[z_] := 
  Evaluate[With[{p = zThk1, v1 = v1, v2 = v2}, 
    If[z < p, v1, v2]]] (*Variable density*heatcapacity*);
A1[z_] := 
  Evaluate[With[{p = zThk1, abs1 = 0, abs2 = 1}, 
    If[z < p, abs1, 
     abs2]]] (*Do we absorb forcing function in the material.  1=yes, \
0=no*);

(*IBC and BCS*)
u03D[z_] := 
  Evaluate[With[{p = zThk1, T1 = Tinit1, T2 = Tinit2}, 
    If[z < p, T1, T2]]] (*Starting temperature profile*);
ibcFor3D = u[x, y, z, 0.] == u03D[z] (*Initial boundary condition*);
bc1For3D = u[x, y, 0, t] == Tbottom(*Temperature at z=0*);
NMBC1 = NeumannValue[hconv*(Tair - u[x, y, z, t]), 
   z == zThkTot] (*Robin Boundary, note:do not include K*);
NMBC2 = NeumannValue[0., x <= 0.](*No heat flux at x=0*);
NMBC3 = NeumannValue[0., x >= xend](*No heat flux at x=xend*);
NMBC4 = NeumannValue[0., y <= 0.](*No heat flux at y=0*);
NMBC5 = NeumannValue[0., y >= ywidth](*No heat flux at y=ywidth*);
opts1 = Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"FiniteElement"}};
heateq3D1 = 
  M1[z]*D[u[x, y, z, t], t] == 
   C1[z]* D[u[x, y, z, t], x, x] + C1[z]* D[u[x, y, z, t], y, y] + 
    C1[z]* D[u[x, y, z, t], z, z] + NMBC1 + NMBC2 + NMBC3 + NMBC4 + 
    NMBC5 + A1[z]*Forcing[x, y, z, t];

eqns1 = {heateq3D1, ibcFor3D, 
   bc1For3D}(*The PDE and conditions to solve for*);

(*Solve 3D heat equation and monitor*)
n = 0;
TempSol3D = 
  Monitor[NDSolveValue[eqns1, 
    u, {x, y, z} \[Element] Mesh0, {t, 0, tend}, opts1, 
    MaxStepSize -> TempMaxStepSize, 
    StepMonitor :> (p = t; ++n)] , {StringForm[
     "Iter:`1`, Time(s):`2`", n, p]}];

(*Create 2D solution using 3D solution*)
TempSol2D[x_, z_, t_] = TempSol3D[x, ywidth/2, z, t]

(*View contour of 3D solution at y=ywidth/2*)
Manipulate[
 ContourPlot[
  TempSol3D[x, ywidth/2, z, t], {x, 0, xend}, {z, 0, zThkTot}, 
  PlotLegends -> Automatic, PlotRange -> Full], {t, 0, tend}]

(*View contour of 2D solution*)
Manipulate[
 ContourPlot[TempSol2D[x, z, t], {x, 0, xend}, {z, 0, zThkTot}, 
  PlotLegends -> Automatic, PlotRange -> Full], {t, 0, tend}]

(------ EDIT on 2018_03_22 ------) Based on user21's suggestion. I tested how long it takes to use the 2D interpolation function derived from a 3D interpolation function and compared it to a 2D interpolation function derived from solving a 2D mesh.

Setup

Needs["NDSolve`FEM`"]
Clear[Derivative]
Clear["Global`*"]

BoolMeshRefine = False; (*Refine the mesh?*)

Tair = 23; (*Fluid temperature of the air*)
Tinit2 = 23; (*Region 1 initial temperature*)
Tinit1 = 80; (*Region 2 initial temperature*)
Tbottom = 80; (*Temperature at the end of z*)
k2 = 0.100; (*Region 1 conductivity*)
rho2 = 1500; (*Region 1 density*)
Cp2 = 1200; (*Region 1 heat capacity*)
k1 = 2.0; (*Region 2 conductivity*)
rho1 = 2200; (*Region 2 density*)
Cp1 = 700;(*Region 2 heat capacity*)
hconv = 6.5; (*Convection coeff*)

xend = 0.005; (*Total length of mesh*)
ywidth = 0.002; (*Width of mesh*)
zThk1 = 0.001;(*Thickness of material 1*)
zThk2 = 0.001;(*Thickness of material 2*)

tend = 10; (*End time for simulation*)

zThkTot = zThk1 + zThk2; (*Total thickness of the materials*)

ForcingVelx = 0.0025; (*Velocity in the x dir m/s*)
ForcingFreq = 0.1; (*Amplitude frequency in Hz*)
ForcingAmp = 1*10^9; (*Power density in W/m^3*)

(*TEMPERATURE solving variables*)
TempMaxStepSize = 0.50; (*in seconds*)

Create Mesh for 3D and 2D control volume

(*MESH variables*)
MeshCellSize = 0.00020; (*Default cell size*)
MeshCellRefSize = 0.00005; (*Refinement size*)
MeshRefDepth = 0.00025; (*Depth from the top of the material 1*)
MeshRefInterfaceWidth = 0.0005; (*Width for refinement where material \
1 meets material 2*)

(*Generate MESH(s)=========================================*)
MeshRefine0 = 
 Function[{vertices, area}, 
  Block[{x, y, z}, {x, y, z} = Mean[vertices];
   If[(z > (zThkTot - MeshRefDepth)) || (zThk1 - 
        MeshRefInterfaceWidth/2 < z < 
       zThk1 + MeshRefInterfaceWidth/2), area > (MeshCellRefSize)^3, 
    area > (MeshCellSize)^3]]];
MeshRefine1 = 
  Function[{vertices, area}, Block[{x, z}, {x, z} = Mean[vertices];
    If[(z > (zThkTot - MeshRefDepth)) || (zThk1 - 
         MeshRefInterfaceWidth/2 < z < 
        zThk1 + MeshRefInterfaceWidth/2), area > (MeshCellRefSize)^2, 
     area > (MeshCellSize)^2]]];
(*Create the mesh(s)*)
If[BoolMeshRefine == True,
  Mesh0 = 
    ToElementMesh[Cuboid[{0, 0, 0}, {xend, ywidth, zThkTot}], 
     MaxCellMeasure -> MeshCellSize^3, 
     MeshRefinementFunction -> MeshRefine0];,
  Mesh0 = 
    ToElementMesh[Cuboid[{0, 0, 0}, {xend, ywidth, zThkTot}], 
     MaxCellMeasure -> MeshCellSize^3];
  ];
If[BoolMeshRefine == True,
  (*Mesh1=ToElementMesh[Cuboid[{0,0},{xend,zThkTot}],MaxCellMeasure->
  MeshCellSize^2,MeshRefinementFunction\[Rule]MeshRefine1];,
  Mesh1=ToElementMesh[Cuboid[{0,0},{xend,zThkTot}],MaxCellMeasure->
  MeshCellSize^2];*)
  Mesh1 = 
    ToElementMesh[Rectangle[{0, 0}, {xend, zThkTot}], 
     MaxCellMeasure -> MeshCellSize^2, 
     MeshRefinementFunction -> MeshRefine1];,
  Mesh1 = 
    ToElementMesh[Rectangle[{0, 0}, {xend, zThkTot}], 
     MaxCellMeasure -> MeshCellSize^2];
  ];

(*Display mesh(s)*)
Show[Mesh0["Wireframe"]] (*Show the 3D mesh*)
Show[Mesh1["Wireframe"]] (*Show the 2D mesh*)

Using a forcing function which varies volumetrically with time (sinusoidal), it is a max at the top of the mesh and decays exponentially through z (vertical dimension).

(*Forcing function to apply heat*)
(*Forcing[x_,y_,z_, \
t_]=ForcingAmp*Sin[((ForcingVelx*t-x)/xend)*2*Pi];*)
Forcing[x_, y_, z_, t_] = (E^(10000*(z - zThkTot)))*ForcingAmp*
   Sin[ForcingFreq*t*2*Pi];
Animate[ContourPlot[
  Forcing[x, ywidth/2, z, t], {x, 0, xend}, {z, 0, zThkTot}, 
  PlotLegends -> Automatic, PlotRange -> All, PlotRange -> Full, 
  PlotPoints -> 100, FrameLabel -> {"x(m)", "z(m)"}, 
  PlotLabel -> "Forcing function (W/m^3)"], {t, 0, tend, 0.1}] 

Setup both the 3D and 2D heat equations

(*Create functions to give different properties to the different \
materials*)
v1 = rho1*Cp1 (*Combination of density and heat capacity*);
v2 = rho2*Cp2;
C1[z_] := 
  Evaluate[With[{p = zThk1, k1 = k1, k2 = k2}, 
    If[z < p, k1, k2]]] (*Variable conductivity*);
M1[z_] := 
  Evaluate[With[{p = zThk1, v1 = v1, v2 = v2}, 
    If[z < p, v1, v2]]] (*Variable density*heatcapacity*);
A1[z_] := 
  Evaluate[With[{p = zThk1, abs1 = 0, abs2 = 1}, 
    If[z < p, abs1, 
     abs2]]] (*Do we absorb forcing function in the material.  1=yes, \
0=no*);

(*IBC and BCS, 3D*)
u0[z_] := 
  Evaluate[With[{p = zThk1, T1 = Tinit1, T2 = Tinit2}, 
    If[z < p, T1, T2]]] (*Starting temperature profile*);
ibcFor3D = u[x, y, z, 0.] == u0[z] (*Initial boundary condition*);
bc1For3D = u[x, y, 0, t] == Tbottom(*Temperature at z=0*);
NMBC1 = NeumannValue[hconv*(Tair - u[x, y, z, t]), 
   z == zThkTot] (*Robin Boundary, note:do not include K*);
NMBC2 = NeumannValue[0., x <= 0.](*No heat flux at x=0*);
NMBC3 = NeumannValue[0., x >= xend](*No heat flux at x=xend*);
NMBC4 = NeumannValue[0., y <= 0.](*No heat flux at y=0*);
NMBC5 = NeumannValue[0., y >= ywidth](*No heat flux at y=ywidth*);
opts = Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"FiniteElement"}};
heateq3D = 
  M1[z]*D[u[x, y, z, t], t] == 
   C1[z]* D[u[x, y, z, t], x, x] + C1[z]* D[u[x, y, z, t], y, y] + 
    C1[z]* D[u[x, y, z, t], z, z] + NMBC1 + NMBC2 + NMBC3 + NMBC4 + 
    NMBC5 + A1[z]*Forcing[x, y, z, t];

eqns3D = {heateq3D, ibcFor3D, 
   bc1For3D}(*The PDE and conditions to solve for*);

(*IBC and BCS, 2D*)
ibcFor2D = u[x, z, 0.] == u0[z] (*Initial boundary condition*);
bc1For2D = u[x, 0, t] == Tbottom(*Temperature at z=0*);
NMBC6 = NeumannValue[hconv*(Tair - u[x, z, t]), 
   z == zThkTot] (*Robin Boundary, note:do not include K*);
NMBC7 = NeumannValue[0., x <= 0.](*No heat flux at x=0*);
NMBC8 = NeumannValue[0., x >= xend](*No heat flux at x=xend*);

heateq2D = 
  M1[z]*D[u[x, z, t], t] == 
   C1[z]* D[u[x, z, t], x, x] + C1[z]* D[u[x, z, t], z, z] + NMBC6 + 
    NMBC7 + NMBC8 + A1[z]*Forcing[x, ywidth/2, z, t];

eqns2D = {heateq2D, ibcFor2D, 
   bc1For2D}(*The PDE and conditions to solve for*);

Solve the heat equations and obtain the temperature profiles

(*Solve 3D heat equation and monitor*)
n = 0;
TempSol3D = 
  Monitor[NDSolveValue[eqns3D, 
    u, {x, y, z} \[Element] Mesh0, {t, 0, tend}, opts, 
    MaxStepSize -> TempMaxStepSize, 
    StepMonitor :> (p = t; ++n)] , {StringForm[
     "Iter:`1`, Time(s):`2`", n, p]}];

(*Solve 2D heat equation and monitor*)
n = 0;
TempSol2D = 
  Monitor[NDSolveValue[eqns2D, 
    u, {x, z} \[Element] Mesh1, {t, 0, tend}, opts, 
    MaxStepSize -> TempMaxStepSize, 
    StepMonitor :> (p = t; ++n)] , {StringForm[
     "Iter:`1`, Time(s):`2`", n, p]}];

(*Create 2D solution using 3D solution*)
TempSol2Dfrom3D[x_, z_, t_] = TempSol3D[x, ywidth/2, z, t]

(*View contour of 3D temperature solution at y=ywidth/2*)
Manipulate[
 ContourPlot[
  TempSol3D[x, ywidth/2, z, t], {x, 0, xend}, {z, 0, zThkTot}, 
  PlotLegends -> Automatic, PlotRange -> Full, PlotPoints -> 100, 
  FrameLabel -> {"x(m)", "z(m)"}, 
  PlotLabel -> "Temperature (C), 3D solution, y=ywidth/2"], {t, 0, 
  tend}]

(*View contour of 2D temperature solution from 3D*)
Manipulate[
 ContourPlot[TempSol2Dfrom3D[x, z, t], {x, 0, xend}, {z, 0, zThkTot}, 
  PlotLegends -> Automatic, PlotRange -> Full, PlotPoints -> 100, 
  FrameLabel -> {"x(m)", "z(m)"}, 
  PlotLabel -> 
   "Temperature (C), 2D solution (from3D), y=ywidth/2"], {t, 0, tend}]

(*View contour of 2D temperature solution*)
Manipulate[
 ContourPlot[TempSol2D[x, z, t], {x, 0, xend}, {z, 0, zThkTot}, 
  PlotLegends -> Automatic, PlotRange -> Full, PlotPoints -> 100, 
  FrameLabel -> {"x(m)", "z(m)"}, 
  PlotLabel -> "Temperature (C), 2D solution"], {t, 0, tend}]

Use both temperature solutions in an Arrhenius accumulation function and time how long it takes to solve for alpha (a value that starts near 0 and rise to 1 depending on the temperature/time effects).

(*Integrate each temperature profile using a simple Arrhenius \
accumulation equation, only operating on the top material*)
Kcoeff = 1*10^6; (*Amplitude coeff*)
Ecoeff = 80000 ;(*Activation energy*)
AccumFnc[\[Alpha]_, T_] = 
 Kcoeff*E^(-Ecoeff/(8.314*
       T)) (1 - \[Alpha]) (*d\[Alpha]/dt function which integrates \
from 0-1 based on temperature profile*)

(*Initial condition for all accumulation integrations*)
IC\[Alpha]3D = \[Alpha][x, z, 0] == 
   1*10^-9; (*\[Alpha] starts out very small*)

(*Integrate the AccumEqn function using the 2D temperature profile \
derived from the 3D temperature profile*)
AccumEqn2Dfrom3D = 
  D[\[Alpha][x, z, t], t] == 
   AccumFnc[\[Alpha][x, z, t], 
    TempSol2Dfrom3D[x, z, t] + 273.15]; (*d\[Alpha]/dt=AccumFnc*)

(*Solve using 2D from 3D*)
n = 0;
TimeBefore1 = AbsoluteTime[];
AccumSolFrom3D = 
  Monitor[NDSolveValue[{AccumEqn2Dfrom3D, IC\[Alpha]3D}, \[Alpha], {x,
      0, xend}, {z, zThk1, zThkTot}, {t, 0, tend},
    MaxStepSize -> 0.50, 
    Method -> {"MethodOfLines", "TemporalVariable" -> t, 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MaxPoints" -> 100, "MinPoints" -> 100}},
    StepMonitor :> (p = t; ++n)], {StringForm["Iter:`1`, Time(s):`2`",
      n, p]}](*solve the PDE*);
TimeAfter1 = AbsoluteTime[];
Print[TimeAfter1 - TimeBefore1]

(*Integrate the AccumEqn function using the 2D temperature profile*)
AccumEqn2D = 
 D[\[Alpha][x, z, t], t] == 
  AccumFnc[\[Alpha][x, z, t], 
   TempSol2D[x, z, t] + 273.15]; (*d\[Alpha]/dt=AccumFnc*)

(*Solve using 2D*)
n = 0;
TimeBefore2 = AbsoluteTime[];
AccumSol2D = 
  Monitor[NDSolveValue[{AccumEqn2D, IC\[Alpha]3D}, \[Alpha], {x, 0, 
     xend}, {z, zThk1, zThkTot}, {t, 0, tend},
    MaxStepSize -> 0.50, 
    Method -> {"MethodOfLines", "TemporalVariable" -> t, 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MaxPoints" -> 100, "MinPoints" -> 100}},
    StepMonitor :> (p = t; ++n)], {StringForm["Iter:`1`, Time(s):`2`",
      n, p]}](*solve the PDE*);
TimeAfter2 = AbsoluteTime[];
Print[TimeAfter2 - TimeBefore2]

Verify the end results for both routines. Alpha should be the same (and it is).

(*Observe \[Alpha]'s for both the 3D derived and direct 2D \[Alpha] \
solution*)
(*View contour of \[Alpha] 2D solution from 3D*)
Manipulate[
 ContourPlot[
  AccumSolFrom3D[x, z, t], {x, 0, xend}, {z, zThk1, zThkTot}, 
  PlotLegends -> Automatic, PlotRange -> Full, PlotPoints -> 100, 
  FrameLabel -> {"x(m)", "z(m)"}, 
  PlotLabel -> "\[Alpha], 2D solution (from3D), y=ywidth/2"], {t, 0, 
  tend}]

(*View contour of \[Alpha] 2D solution from 3D*)
Manipulate[
 ContourPlot[AccumSol2D[x, z, t], {x, 0, xend}, {z, zThk1, zThkTot}, 
  PlotLegends -> Automatic, PlotRange -> Full, PlotPoints -> 100, 
  FrameLabel -> {"x(m)", "z(m)"}, 
  PlotLabel -> "\[Alpha], 2D solution"], {t, 0, tend}]

Using my computer, I obtain a solve time of 104 seconds for alpha using the 2D temperature interpolation function that is "derived" from the 3D interpolation function. If I use the 2D temperature interpolation function which is obtained from solving the heat equation on a 2D mesh, it only takes 21 seconds to solve for alpha. Another thing to note is the massive difference in memory it takes to store the 3D vs 2D interpolation function.

(------ EDIT on 2018_03_23 ------) So I think I've found a solution. It involves recreating the interpolation function using a grid resolution that is the mesh resolution, the "Table" function, and the time array from the 3D interpolation function. It adds a little bit of overhead but will save me a massive amount of time in the end. The routine is:

    (*Grab time array used in interpolation function*)
    Mesh3Dinfo = TempSol3D["Coordinates"] [[1]];
    Mesh3Dtime = TempSol3D["Coordinates"] [[2]];

(*use the time array from 3D mesh*)
    TempSol2DreducedTab = 
     Table[{x, z, Mesh3Dtime[[i]], 
       TempSol3D[x, ywidth/2, z, Mesh3Dtime[[i]]]}, {x, 0, xend, 
       MeshCellSize}, {z, 0, zThkTot, MeshCellSize}, {i, 1, 
       Length[Mesh3Dtime], 1}] 


(*Create new interpolation function*)    
TempSol2Dreduced = 
     Interpolation[Flatten[TempSol2DreducedTab, 2], 
      InterpolationOrder -> {3, 2, 
        2}] (*Flatten the table and create the 2D interpolation function*)

    (*View contour of reduced 2D temperature solution*)
    ContourPlot[
     TempSol2Dreduced[x, z, tend], {x, 0, xend}, {z, 0, zThkTot}, 
     PlotLegends -> Automatic, PlotRange -> Full, PlotPoints -> 100, 
     FrameLabel -> {"x(m)", "z(m)"}, PlotLabel -> "Reduced 2D solution"]

    (*Test to see how long new derived 2D interpolation function takes to \
    solve*)

    (*Integrate the AccumEqn function using the new 2D temperature \
    profile, reduced from 3D*)
    AccumEqn2Dnew = 
     D[\[Alpha][x, z, t], t] == 
      AccumFnc[\[Alpha][x, z, t], 
       TempSol2Dreduced[x, z, t] + 273.15]; (*d\[Alpha]/dt=AccumFnc*)

    (*Solve using new 2D interpolation function derived by creating a \
    table*)
    n = 0;
    TimeBefore3 = AbsoluteTime[];
    AccumSol2Dnew = 
      Monitor[NDSolveValue[{AccumEqn2Dnew, IC\[Alpha]3D}, \[Alpha], {x, 0,
          xend}, {z, zThk1, zThkTot}, {t, 0, tend},
        MaxStepSize -> 0.50, 
        Method -> {"MethodOfLines", "TemporalVariable" -> t, 
          "SpatialDiscretization" -> {"TensorProductGrid", 
            "MaxPoints" -> 100, "MinPoints" -> 100}},
        StepMonitor :> (p = t; ++n)], {StringForm["Iter:`1`, Time(s):`2`",
          n, p]}](*solve the PDE*);
    TimeAfter3 = AbsoluteTime[];
    Print[TimeAfter3 - TimeBefore3]

    (*View contour of \[Alpha] reduced 2D solution*)
    Manipulate[
     ContourPlot[
      AccumSol2Dnew[x, z, t], {x, 0, xend}, {z, zThk1, zThkTot}, 
      PlotLegends -> Automatic, PlotRange -> Full, PlotPoints -> 100, 
      FrameLabel -> {"x(m)", "z(m)"}, 
      PlotLabel -> "\[Alpha], reduced 2D solution"], {t, 0, tend}]

At the end I compared the time it took to use the temperature 2D interpolation functions in another NDSolve routine. To make the comparison more fair I had to match the order of the interpolation functions from the derived 2D function and 2D function from a 2D mesh (Which is 3,2,2 or x,z,t). The reduced 2D interp function is using the Hermite method instead of Unstructured (haven't figured out how to get around that) which might add some overhead. Also all temperature and alpha values are extremely close for all methods. The time it takes to use these functions are:

2D interp derived from 3D interp (is still actually a 3D interp) = 83.76 sec.

2D interp from 2D mesh = 17.92 sec.

2D interp reduced from 3D interp = 27.32 sec.

$\endgroup$
  • $\begingroup$ What’s the time difference between using 2D and slice of 3D? $\endgroup$ – Vsevolod A. Mar 22 '18 at 2:20
  • $\begingroup$ I can't comment on the time it takes with a 2D interpolation function (because I cannot obtain it =/) but when using the 3D interpolation function (using the full script and extended process) it takes about 8 hours for 1 solution. I'm only using the above script as a demo. The problem is I am doing a gradient descent routine which requires using the 2D interpolating function in many NDSolve loops. If I was to guess from creating 2D interpolating functions (just to play with) and then using them in my process, then I would say the time difference is maybe 1/10th. $\endgroup$ – Dr. Quick Mar 22 '18 at 2:40
  • 2
    $\begingroup$ It is probably possible to do something along the lines you think, but it might be a bit involved. However, I am not certain that this will decrease the time for integration. You should conduct a simple experiment, make a 3D interpolation function and integrate over and a similarly coarse 2D version and compare the time it takes to integrate over it and see if such an approach would be worthwhile in the first place. You could add your findings to your post. $\endgroup$ – user21 Mar 22 '18 at 6:43
  • $\begingroup$ user21, I edited my code to do exactly what you said. I'll post it above so that anyone can confirm. The 2D interpolation function (from 3D) took 104 seconds to solve when used in an Arrhenius equation vs 21 seconds for the 2D interpolation function (obtained solving a 2D mesh). Time is not the only problem. My main script will solve for about 30 interpolation functions and put them together in a piecewise function. This piecewise is about 70 GB when solving with refined mesh. If I can strip off the "slice" in a 2D interpolation function, then the memory usage would be drastically reduced. $\endgroup$ – Dr. Quick Mar 22 '18 at 23:41
  • $\begingroup$ I have looked a bit at this. Here is a simpler test: AbsoluteTiming[ res1 = Table[ Block[{x, z, t = 0.1}, {x, z} = c; TempSol2Dfrom3D[x, z, t]], {c, Mesh1["Coordinates"]}];] compared to AbsoluteTiming[ res2 = Table[ Block[{x, z, t = 0.1}, {x, z} = c; TempSol2D[x, z, t]], {c, Mesh1["Coordinates"]}];]. The second one is about a factor of 5 faster. And the norm between the two seems reasonable: Norm[res1 - res2]. Which gives me: 0.000951001. Now, to get a 2D+time interpolating function from a 3D + time function will need some sort of a cut plane mechanism.... $\endgroup$ – user21 Mar 23 '18 at 7:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.