Why can't I use Inactive[D] in NDSolve?

Question

I have a 2D rectangular domain with -L/2 < x < L/2 and 0 < y < D0, which is divided into two parts, with a sinusoidal curve y == 1 + δ Cos[2 Pi x/L] separating the two.

I want to solve Laplace equation D[κ[x, y] D[T[x, y], x], x] + D[κ[x, y] D[T[x, y], y], y] == 0, where κ[x, y] takes different values in the two regions (equals 1 below the curve, and 0.2422 above the curve). After copying some codes in StackExchange, I have the following codes which I think give the right answer:

D0 = 3;
k = 0.2422;
L = 3;
δ = 0.2;
sol = NDSolveValue[{Inactive[
   Div][-If[y <= 1 - δ Cos[2 π x/L], 1, 
     k].Inactive[Grad][T[x, y], {x, y}], {x, y}] == 0, 
PeriodicBoundaryCondition[T[x, y], x == -L/2, 
 TranslationTransform[{L, 0}]],
DirichletCondition[T[x, y] == 1 + δ Cos[2 π x/L], 
 y == 0 && -L/2 < x < L/2],
DirichletCondition[T[x, y] == 0, y == D0 && -L/2 < x < L/2]},
   T, {x, -L/2, L/2}, {y, 0, D0}, 
   Method -> {"FiniteElement", "InterpolationOrder" -> {T -> 2}, 
 "MeshOptions" -> {"MaxCellMeasure" -> 0.0005}, 
 "IntegrationOrder" -> 5}];
Show[ContourPlot[sol[x, y], {x, -L/2, L/2}, {y, 0, D0}, 
  PlotLegends -> Automatic, Contours -> 20],
 Plot[1 - δ Cos[2 π x/L], {x, -L/2, L/2}]]

Now I want to solve a slightly different equation D[κ[x, y] D[T[x, y], x], x] + epsilon^2 D[κ[x, y] D[T[x, y], y], y] == 0. This is basically Laplace equation where the two spatial coordinates are scaled differently, and I want to study the effect of ignoring the second term with epsilon^2 in my analytical work.

I try to define the equation differently as Inactive[D][-If[y <= 1-δ Cos[2 π x/L], 1, k] Inactive[D][T[x, y], {x}], {x}] + (epsilon^2) Inactive[D][-If[y <= 1-δ Cos[2 π x/L], 1, k] Inactive[D][T[x, y], {y}], {y}] == 0

But it gives an error NDSolveValue::derivs: No derivatives of dependent variables were found in the equations. NDSolveValue is designed to solve differential or differential algebraic equations. Use NSolve or FindRoot to numerically solve algebraic equations.

What's the problem? In fact, I don't think I fully understand the use of Inactive. What's the best way to define the PDE (or the domain) in this case?

Answer 1

The coefficient you give for the diffusion operator is a scalar but you can give a matrix of coefficients. I think this is what you are looking for:

D0 = 3;
k = 0.2422;
L = 3;
δ = 0.2;
epsilon = 10^-1;

Note how the coeff now returns a matrix of coefficients:

coeff = If[
   y <= 1 - δ Cos[2 π x/L], -{{1, 0}, {0, 
      epsilon^2*1}}, -{{k, 0}, {0, epsilon^2*k}}];

Well just plug that in:

sol = NDSolveValue[{Inactive[Div][
      coeff.Inactive[Grad][T[x, y], {x, y}], {x, y}] == 0, 
    PeriodicBoundaryCondition[T[x, y], x == -L/2, 
     TranslationTransform[{L, 0}]], 
    DirichletCondition[T[x, y] == 1 + δ Cos[2 π x/L], 
     y == 0 && -L/2 < x < L/2], 
    DirichletCondition[T[x, y] == 0, y == D0 && -L/2 < x < L/2]}, 
   T, {x, -L/2, L/2}, {y, 0, D0}, 
   Method -> {"FiniteElement", 
     "MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}];

As a side note, for single equations you never need to specify "InterpolationOrder"; that is only relevant for coupled PDE where you want to interpolate some of the dependent variables not at the highest level.

The new in version 12.1 HeatTransfer tutorial has a section on isotropic, orthotropic and anisotropic diffusion that explain the concepts and implementation more.