9
$\begingroup$

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?

$\endgroup$
3
  • 3
    $\begingroup$ Where's the epsilon in the new PDE? Have you read the "Formal Partial Differential Equations" section in the document page FEMDocumentation/tutorial/FiniteElementOverview? $\endgroup$
    – xzczd
    Mar 27, 2020 at 5:18
  • $\begingroup$ @xzczd just edited. I basically don't understand why the way I write the equation (splitting it into the x and y derivatives instead of using grad and div) doesn't work. $\endgroup$
    – Physicist
    Mar 27, 2020 at 5:49
  • 1
    $\begingroup$ I think this is a good question. My guess is Inactive is only for building the "formal PDE", but it's not explicitly mentioned in the document. Let's wait for user21's clarification :) . $\endgroup$
    – xzczd
    Mar 27, 2020 at 5:51

1 Answer 1

10
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ Any comment on the Inactive[D] part? $\endgroup$
    – xzczd
    Mar 27, 2020 at 11:21
  • $\begingroup$ @xzczd, basically because it is not implemented. You should be able to do everything with Inactive[Div] and Inactive[Grad], no? $\endgroup$
    – user21
    Mar 27, 2020 at 11:58
  • $\begingroup$ Yeah, but Inactive[D] will make coding more flexible, right? :) I think the more important part is, this isn't explicitly documented. Do you mean currently only Inactive[Div] and Inactive[Grad] are supported by FiniteElement method of NDSolve? $\endgroup$
    – xzczd
    Mar 27, 2020 at 12:06
  • 1
    $\begingroup$ @xzczd, yes only Inactive Div/Grad (for differential operators) are supported. And, yes, it might make sense to mention that in the documentation. $\endgroup$
    – user21
    Mar 27, 2020 at 12:42
  • 1
    $\begingroup$ @xzczd, I tried to clarify this in the documentation and this should be in 12.1.1. Thanks. $\endgroup$
    – user21
    Mar 30, 2020 at 5:59

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.