Skip to main content

plotting - Solving Cahn-Hilliard equation: LinearSolve: Linear equation encountered that has no solution


I have built the Cahn-Hilliard Eqs. in MMA (Mixed Formulation, second order), However, it doesnot work in MMA using Finite Element.


LinearSolve: Linear equation encountered that has no solution.


And "... are not the same shape".


Theory & numerical formulation based on this FEniCS Benchmark Test enter link description here


My code:



(*Initial Parameters*)Needs["NDSolve`FEM`"];
Mobi = 1.0; lame = 0.01; noise = 0.02; conu0 = 0.63;
xmax = 1.0;
ymax = 1.0;
tmax = 1.0;

Ω = Rectangle[{0, 0}, {a, b}] /. {a -> 1, b -> 1};
RegionPlot[Ω, AspectRatio -> Automatic]
mesh = ToElementMesh[Ω, "MaxCellMeasure" -> 1/1000, "MeshElementType" -> QuadElement];
mesh["Wireframe"]

n = Length[mesh["Coordinates"]]
u0 = ElementMeshInterpolation[{mesh}, conu0 + noise*(0.5 - RandomReal[{0, 1}, n])];
Plot3D[u0[x, y], {x, y} ∈ mesh]

op1 = D[u[t, x, y], t] - Laplacian[v[t, x, y], {x, y}] Mobi

op2 = v[t, x, y] - 200 u[t, x, y] (1 - 3 u[t, x, y] + 2 u[t, x, y]^2) + 
      lame Laplacian[u[t, x, y], {x, y}]

{unn, vnn} = 

  NDSolve[{op1 == 0, op2 == 0, u[0, x, y] == u0[x, y], 
    v[0, x, y] == 0}, {u, v}, {t, 0, tmax}, {x, y} ∈ mesh];

Answer



I can offer an easy-to-implement explicit method of Euler using FEM and NDSolve. Here we used a test example like on Python from https://fenicsproject.org/olddocs/dolfin/1.3.0/python/demo/documented/cahn-hilliard/python/documentation.html#. The output picture is about the same. These are the initial data, equations, and parameters.


<< NDSolve`FEM`
Lx = 1; Ly = 1; nn = 50; t0 = 5*10^-6;
reg = Rectangle[{0, 0}, {1, 1}];

f[x_] := 100 x^2 (1 - x)^2
lambd = 1/100; noise = 0.02; conu0 = 0.63;

M = 1;
thet = 1/2;
eq1 = D[c[t, x, y], t] - Div[M Grad[u[t, x, y], {x, y}], {x, y}] == 0;
eq2 = u[t, x, y] - D[f[c[t, x, y]], c[t, x, y]] + 
    lambd Laplacian[c[t, x, y], {x, y}] == 0;
mesh = ToElementMesh[reg, "MaxCellMeasure" -> 1/1000, 
      "MeshElementType" -> QuadElement];
mesh["Wireframe"]
n = Length[mesh["Coordinates"]];
u0 = ElementMeshInterpolation[{mesh}, 

      conu0 + noise*(0.5 - RandomReal[{0, 1}, n])];
uf[0][x_, y_] := 0
cf[0][x_, y_] := u0[x, y]
Plot3D[u0[x, y], {x, y} \[Element] mesh]

This is the implementation of the explicit Euler.


eq = {-Laplacian[u[x, y], {x, y}] + (c[x, y] - cf[i - 1][x, y])/t0 == 
   NeumannValue[0, True], -200 (1 - cf[i - 1][x, y])^2 c[x, y] + 
    200 (1 - c[x, y]) cf[i - 1][x, y]^2 + u[x, y] + 
    1/100 Laplacian[c[x, y], {x, y}] == 

   NeumannValue[0, True]}; Do[{cf[i], uf[i]} = 
   NDSolveValue[eq, {c, u}, {x, y} \[Element] mesh] // Quiet;, {i, 1, 
  nn}]

This is an animation and 3D image.


frame = Table[
   DensityPlot[cf[i][x, y], {x, y} \[Element] mesh, 
    ColorFunction -> "Rainbow", Frame -> False, 
    PlotLabel -> Row[{"t = ", i t0 1.}]], {i, 0, nn, 2}];


ListAnimate[frame]
Plot3D[cf[50][x, y], {x, y} \[Element] mesh, PlotRange -> All, 
 Mesh -> None, ColorFunction -> "Rainbow"]

Figure 1


I managed to debug code @Henrik Schumacher, so that with equal parameters and the same input data, similar results are obtained with code above and with code @Henrik Schumacher. Thus, code @Henrik Schumacher passed the test for Python.


Henrik Schumacher debugged code:


Needs["NDSolve`FEM`"];
Mobi = 1.0; lame = 0.01; noise = 0.02; conu0 = 0.63;
xmax = 1.0;

ymax = 1.0;
tmax = 1.0;
a = 1.;
b = 1.;

\[CapitalOmega] = Rectangle[{0, 0}, {a, b}];
mesh = ToElementMesh[\[CapitalOmega], "MaxCellMeasure" -> 1/5000, 
  "MeshElementType" -> QuadElement, "MeshOrder" -> 1]

ClearAll[x, y, u];

vd = NDSolve`VariableData[{"DependentVariables", 
     "Space"} -> {{u}, {x, y}}];
sd = NDSolve`SolutionData[{"Space"} -> {mesh}];
cdata = InitializePDECoefficients[vd, sd, 
   "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}, 
   "MassCoefficients" -> {{1}}];
bcdata = InitializeBoundaryConditions[vd, 
   sd, {{DirichletCondition[u[x, y] == 0., True]}}];
mdata = InitializePDEMethodData[vd, sd];


(*Discretization*)
dpde = DiscretizePDE[cdata, mdata, sd];
dbc = DiscretizeBoundaryConditions[bcdata, mdata, sd];
{load, A, damping, M} = dpde["All"];
(*DeployBoundaryConditions[{load,A},dbc];*)
(*DeployBoundaryConditions[{load,M},dbc];*)
\[Theta] = 1;
\[Tau] = 0.000005;
\[Mu] = Mobi;
\[Lambda] = lame;

L = ArrayFlatten[{{M, \[Tau] \[Mu] \[Theta] A}, {-\[Lambda] A, M}}];
n = Length[mesh["Coordinates"]];
m = 50;
f = x \[Function] 100. x^2 (1. - x^2);
Df = x \[Function] Evaluate[f'[x]];
rhs[u_, v_] := 
  Join[M.u - (\[Mu] \[Tau] (1. - \[Theta])) A.v, 
   M.(200 (1 - u)^2 u - 200 (1 - u) u^2)];
S = LinearSolve[L, Method -> "Pardiso"];


u0 = conu0 + noise*(0.5 - RandomReal[{0, 1}, n]);
ulist = ConstantArray[0., {m, n}];
ulist[[1]] = u = u0;

v0 = 0. rhs[u0, 0. u0][[n + 1 ;; 2 n]];
v = v0;
Do[sol = S[rhs[u, v]];
  ulist[[k]] = u = sol[[1 ;; n]];
  v = sol[[n + 1 ;; 2 n]];, {k, 2, m}];
frames = Table[

   Image[Map[ColorData["Rainbow"], 
     Partition[ulist[[k]], Sqrt[n]], {2}], Magnification -> 3], {k, 1,
     m, 1}];
Manipulate[frames[[k]], {k, 1, Length[frames], 1}, 
 TrackedSymbols :> {k}]

My code (for comparison):


u0i = ElementMeshInterpolation[{mesh}, 
      u0];
uf[0][x_, y_] := 0

cf[0][x_, y_] := u0i[x, y]
DensityPlot[u0i[x, y], {x, y} \[Element] mesh, 
 ColorFunction -> "Rainbow", PlotLegends -> Automatic]
nn = 50; t0 = 
 5*10^-6; eq = {-Laplacian[
      u1[x, y], {x, y}] + (c[x, y] - cf[i - 1][x, y])/t0 == 
   NeumannValue[0, True], -200 (1 - cf[i - 1][x, y])^2 c[x, y] + 
    200 (1 - c[x, y]) cf[i - 1][x, y]^2 + u1[x, y] + 
    1/100 Laplacian[c[x, y], {x, y}] == 
   NeumannValue[0, True]}; Do[{cf[i], uf[i]} = 

   NDSolveValue[eq, {c, u1}, {x, y} \[Element] mesh] // Quiet;, {i, 1,
   nn}]

frame = Table[
   DensityPlot[cf[i][x, y], {x, y} \[Element] mesh, 
    ColorFunction -> "Rainbow", Frame -> False, 
    PlotLabel -> Row[{"t = ", i t0 1.}]], {i, 0, nn, 1}];

ListAnimate[frame]


Comparison of two results


ul = ElementMeshInterpolation[{mesh}, 
     ulist[[nn]]]; {Plot3D[ul[x, y], {x, y} \[Element] mesh, 
  ColorFunction -> "Rainbow", Mesh -> None, 
  PlotLabel -> Row[{"\[Theta] = ", \[Theta]}]], 
 Plot3D[cf[nn][x, y], {x, y} \[Element] mesh, 
  ColorFunction -> "Rainbow", Mesh -> None]}

Figure 2 For $\theta=\frac {1}{2}$ matching is better Figure 3


Another method using NDSolveValue and "MethodOfLines". The code is very slow and with a warning NDSolveValue::ibcinc: Warning: boundary and initial conditions are inconsistent. The result does not match Python and FEM.



<< NDSolve`FEM`
Lx = 1; Ly = 1; nn = 50; t0 = 5*10^-6; tmax = t0 nn;
reg = Rectangle[{0, 0}, {1, 1}];

f[x_] := 100 x^2 (1 - x)^2
lambd = 1/100; noise = 0.02; conu0 = 0.63;
M = 1;
thet = 1/2;
eq1 = D[c[t, x, y], t] - Div[M Grad[u[t, x, y], {x, y}], {x, y}] == 0;
eq2 = u[t, x, y] - D[f[c[t, x, y]], c[t, x, y]] + 

    lambd Laplacian[c[t, x, y], {x, y}] == 0;

mesh = ToElementMesh[reg, "MaxCellMeasure" -> 1/1000, 
      "MeshElementType" -> QuadElement];
mesh["Wireframe"]
n = Length[mesh["Coordinates"]];
u0 = ElementMeshInterpolation[{mesh}, 
      conu0 + noise*(0.5 - RandomReal[{0, 1}, n])];
ic = {c[0, x, y] == u0[x, y], u[0, x, y] == 0};
bc = {Derivative[0, 1, 0][c][t, 0, y] == 0, 

   Derivative[0, 1, 0][c][t, 1, y] == 0, 
   Derivative[0, 1, 0][u][t, 0, y] == 0, 
   Derivative[0, 1, 0][u][t, 1, y] == 0, 
   Derivative[0, 0, 1][c][t, x, 0] == 0, 
   Derivative[0, 0, 1][c][t, x, 1] == 0, 
   Derivative[0, 0, 1][u][t, x, 0] == 0, 
   Derivative[0, 0, 1][u][t, x, 1] == 0};

Monitor[{csol, usol} = 
  NDSolveValue[{eq1, eq2, ic, bc}, {c, u}, {x, 0, 1}, {y, 0, 1}, {t, 

    0, tmax}, 
   Method -> {"IndexReduction" -> Automatic, 
     "EquationSimplification" -> "Residual", 
     "PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"TensorProductGrid", 
         "MinPoints" -> 41, "MaxPoints" -> 81, 
         "DifferenceOrder" -> "Pseudospectral"}}}, 
   EvaluationMonitor :> (monitor = 
      Row[{"t=", CForm[t], " csol=", CForm[c[t, .5, .5]]}])], monitor]


Compare the result with FEM (my code)


uf[0][x_, y_] := 0
cf[0][x_, y_] := u0[x, y]

eq = {-Laplacian[u[x, y], {x, y}] + (c[x, y] - cf[i - 1][x, y])/t0 == 
   NeumannValue[0, True], -200 (1 - cf[i - 1][x, y])^2 c[x, y] + 
    200 (1 - c[x, y]) cf[i - 1][x, y]^2 + u[x, y] + 
    1/100 Laplacian[c[x, y], {x, y}] == 
   NeumannValue[0, True]}; Do[{cf[i], uf[i]} = 
   NDSolveValue[eq, {c, u}, {x, y} \[Element] mesh] // Quiet;, {i, 1, 

  nn}]
{Plot3D[csol[tmax, x, y], {x, 0, 1}, {y, 0, 1}, Mesh -> None, 
  ColorFunction -> "Rainbow"], 
 Plot3D[cf[50][x, y], {x, y} \[Element] mesh, PlotRange -> All, 
  Mesh -> None, ColorFunction -> "Rainbow"]}

On the left fig. 4 the "MethodOfLines", on the right FEM. It can be seen that in the `"MethodOfLines" high-frequency harmonics are added. Figure 4


Comments

Post a Comment

Popular posts from this blog

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

plotting - How to draw lines between specified dots on ListPlot?

I would like to create a plot where I have unconnected dots and some connected. So far, I have figured out how to draw the dots. My code is the following: ListPlot[{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4,13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full] I have thought using ListLinePlot command, but I don't know how to specify to the command to draw only selected lines between the dots. Do have any suggestions/hints on how to do that? Thank you. Answer One possibility would be to use Epilog with Line : ListPlot[ {{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4, 13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full, Epilog -> { Line[ ...
Post a Comment
Read more