NDSolve Transient PDE Mesh

Question

I apologize in advance if this is a somewhat basic question but I can't seem to find the solution I need and my mathematica skills are limited. I am looking to solve a 2D heat equation (for diffusion in a semiconductor) but I want to use a non-uniform mesh to improve my simulation time. I want a dense mesh at the top and a sparse mesh in lower regions of the simulations. Essentially, my question is, how do I specify the mesh for ND solve for a transient PDE?

Thanks, in advance for any help you can provide.

Here is my code with my own mesh defined:

(*Initial Parameters*)
Needs["NDSolve`FEM`"];
N0 = 1*^19;
D0 = 1.49*^-6; (*Diff Coeff for Sn in InAs*)
Q = 1.17; (*diffusion activation energy for temp dep diff*)
kBoltz = 8.617*^-5;(*Boltzmann Constant (eV)*)
T = 800 + 273; (*Diffusion Temp*)
tmax = 120*60;
xmax = 28*^-4;
ymin = -40*^-4;
meshStep = xmax*ymin/(4*^3);
Deff = D0*Exp[-Q/(kBoltz*T)];

(*Define the source Profile Conditions*)
tab = Table[Mod[.0371*x, 1], {x, 0, 28, 1}];
dcR = Round[dc, .1];
sourceTab = {};
Do[
 dc = tab[[i]];
 dcR = Round[dc, .1];
 dutyTab = Table[If[i <= dcR*10, 1, 0], {i, 1, 10, 1}];
 AppendTo[sourceTab, dutyTab],
 {i, 1, Length[tab] - 1}
 ]
sourceTab = Flatten[sourceTab];
AppendTo[sourceTab, 0];
xlist = Table[x, {x, 0, 28*10^-4, 28*^-4/(Length[sourceTab] - 1)}];
xyList = Thread[{xlist, sourceTab}];

iFun = Interpolation[xyList, InterpolationOrder -> 1, 
   Method -> "Piecewise"];

(*Create the Mesh Region*)
bmesh1 = ToBoundaryMesh[
   "Coordinates" -> {{0, 0}, {xmax, 0}, {xmax, ymin}, {0, ymin}, {0, 
      ymin/5}, {xmax, ymin/5}}, 
   "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 
        1}, {5, 6}}]}];
mesh = ToElementMesh[bmesh1, 
   "RegionMarker" -> {{{xmax/2, ymin/10}, 
      1, .0001}, {{xmax/2, ymin/2}, 2, .01}}];


(*Setup and Solve Diff EQ*)
un = NDSolveValue[{
   D[u[x, y, t], t] - (D[u[x, y, t], x, x] + D[u[x, y, t], y, y])*
      Deff == NeumannValue[0, 
     y == 0 && (0 < x <= 2*xmax/5 || 3*xmax/5 <= x < xmax)],
   DirichletCondition[u[x, y, t] == N0, y == 0 && iFun[x] >= 1],
   u[x, y, 0] == 0,
   u[0, y, t] == u[xmax, y, t],
   PeriodicBoundaryCondition[u[x, y, t], x == 0, 
    TranslationTransform[{0, xmax}]]
   }, u,
  {x, 0, xmax}, {y, ymin, 0}, {t, 0, tmax},
  Method -> {"MethodOfLines", "TemporalVariable" -> t, 
    "SpatialDiscretization" -> {"FiniteElement"}}, 
  StartingStepSize -> xmax/10000
  ]

(*Plot Result*)
    Show[
     ContourPlot[un[x, y, tmax], {x, 0, xmax}, {y, ymin, 0},
      PlotPoints -> 50,
      Contours -> Flatten[Table[{1., 5., 10.}*10^n, {n, 15, 20, 1}]],
      ColorFunction -> Function[{f},
        If[f < 10^15, ColorData["Rainbow"][.0],
         ColorData["Rainbow"][(Log[10, f] - 15)/(20 - 15)]]
        ],
      ColorFunctionScaling -> False,
      PlotLegends -> BarLegend[{Function[{f},
          If[f < 10^15,
           ColorData["Rainbow"][.5],
           ColorData["Rainbow"][(Log[10, f] - 15)/(20 - 15)]]],
         {1*^15, 1*^20}}, 
        Flatten[Table[{1., 5., 10.}*10^n, {n, 15, 19, 1}]]],

      PlotRange -> All,
      ImageSize -> Large
      ], un["ElementMesh"]["Wireframe"]]

Answer 1

First initialise the problem data:

(*Initial Parameters*)Needs["NDSolve`FEM`"];
N0 = 1*^19;
D0 = 1.49*^-6;(*Diff Coeff for Sn in InAs*)Q = 1.17;(*diffusion \
activation energy for temp dep diff*)kBoltz = 8.617*^-5;(*Boltzmann \
Constant (eV)*)T = 800 + 273;(*Diffusion Temp*)tmax = 120*60;
xmax = 28*^-4;
ymin = 40*^-4;
meshStep = xmax*ymin/(4*^3);
Deff = D0*Exp[-Q/(kBoltz*T)];

(*Define the source Profile Conditions*)
tab = Table[Mod[.0371*x, 1], {x, 0, 28, 1}];
dcR = Round[dc, .1];
sourceTab = {};
Do[dc = tab[[i]];
 dcR = Round[dc, .1];
 dutyTab = Table[If[i <= dcR*10, 1, 0], {i, 1, 10, 1}];
 AppendTo[sourceTab, dutyTab], {i, 1, Length[tab] - 1}]
sourceTab = Flatten[sourceTab];
AppendTo[sourceTab, 0];
xlist = Table[x, {x, 0, 28*10^-4, 28*^-4/(Length[sourceTab] - 1)}];
xyList = Thread[{xlist, sourceTab}];

iFun = Interpolation[xyList, InterpolationOrder -> 1, 
   Method -> "Piecewise"];

Then create the finite element mesh:

mesh = ToElementMesh[
ImplicitRegion[0 <= x && 0 <= y, {x, y}], {{0, xmax}, {0, ymin}}, 
MeshRefinementFunction -> 
Function[{vertices, area}, 
Block[{x, y}, {x, y} = Mean[vertices]; 
If[y < ymin 0.1 , area > 0.000000001, area > 1]]]];
mesh["Wireframe"]

Solve the problem:

un = NDSolveValue[{D[u[x, y, t], 
      t] - (D[u[x, y, t], x, x] + D[u[x, y, t], y, y])*Deff == 
    NeumannValue[0, 
     y == 0 && (0 < x <= 2*xmax/5 || 3*xmax/5 <= x < xmax)], 
   DirichletCondition[u[x, y, t] == N0, y == 0 && iFun[x] >= 1], 
   u[x, y, 0] == 0, 
   PeriodicBoundaryCondition[u[x, y, t], x == 0, 
    TranslationTransform[{0, xmax}]]}, 
  u, {x, y} \[Element] mesh, {t, 0, tmax}, 
  Method -> {"MethodOfLines", "TemporalVariable" -> t, 
    "SpatialDiscretization" -> {"FiniteElement"}}, 
  StartingStepSize -> xmax/10000]

and plot the solution:

(*Plot Result*)
Show[ContourPlot[un[x, y, tmax], {x, 0, xmax}, {y, ymin, 0}, 
  PlotPoints -> 50, 
  Contours -> Flatten[Table[{1., 5., 10.}*10^n, {n, 15, 20, 1}]], 
  ColorFunction -> 
   Function[{f}, 
    If[f < 10^15, ColorData["Rainbow"][.0], 
     ColorData["Rainbow"][(Log[10, f] - 15)/(20 - 15)]]], 
  ColorFunctionScaling -> False, 
  PlotLegends -> 
   BarLegend[{Function[{f}, 
      If[f < 10^15, ColorData["Rainbow"][.5], 
       ColorData["Rainbow"][(Log[10, f] - 15)/(20 - 15)]]], {1*^15, 
      1*^20}}, Flatten[Table[{1., 5., 10.}*10^n, {n, 15, 19, 1}]]], 
  PlotRange -> All, ImageSize -> Large], 
 un["ElementMesh"]["Wireframe"]]