2
$\begingroup$

I tried numerically solving a set of time-dependent PDEs with variables {u, v, w} by NDSolve over 2 regions with 2D grids, generated by FEM. The solution was exported to hard drive in the format of InterpolatingFunction. But as I import from the saved file, the code fails to retrieve the stored data.

Needs["NDSolve`FEM`"]
dr = 0.5 10^-0; rx = 2./2  10^-0; xs = -(dr + rx); ys = 0.0;
ec = RegionUnion[Disk[{-rx, 0}, dr], Rectangle[{-rx, -dr}, {rx, dr}], 
   Disk[{rx, 0}, dr]];
bmesh = ToBoundaryMesh[ec, "MaxBoundaryCellMeasure" -> .05];
bmesh["Wireframe"];
f = Function[{vertices, area}, 
   area > 0.002 (1. - 0.5 Norm[Mean[vertices]])];
(mesh = ToElementMesh[bmesh, MeshRefinementFunction -> f])["Wireframe"]
ndr = 0.3 10^-0; nrx = 
 1.2/2  10^-0; nxs = -(ndr + nrx); nys = 0.0; nlx = 0.4/1  10^-0;
nc = RegionUnion[Disk[{-nrx, 0}, ndr], 
   Rectangle[{-nrx, -ndr}, {(-nrx + nlx), ndr}], 
   Disk[{(-nrx + nlx), 0}, ndr], Disk[{(nrx - nlx), 0}, ndr], 
   Rectangle[{(nrx - nlx), -ndr}, {nrx, ndr}], Disk[{nrx, 0}, ndr]];
rPC = 0.00  10^-0;
xPC1 = rx + (dr - 1*rPC)*Cos[1 Pi/4]; yPC1 = (dr - 1*rPC)*
  Sin[1 Pi/4];
xPC2 = -rx - (dr - 1*rPC)*Cos[-Pi/3]; yPC2 = (dr - 1*rPC)*Sin[-Pi/3];
PtChg1 = Disk[{xPC1, yPC1}, rPC];
PtChg2 = Disk[{xPC2, yPC2}, rPC];
PtChg = RegionUnion[PtChg1, PtChg2];
Show[Region[Style[nc, Cyan]], mesh["Wireframe"], 
 Graphics[{Magenta, Thick, PointSize[0.01], Point[{xPC1, yPC1}], 
   Point[{xPC2, yPC2}]}]]
Ck0 = 139.20  10^9; Ccl0 = 12.43  10^9; 
u0 = Ck0/Ck0;
v0 = Ccl0/Ck0;
dk = 1.95  10^3; dcl = 2.02  10^3; 
d0 = dcl/dk;
Tend = 0.0001 ;
V2 = -0.05*37.436;  
q = 516.75  10^9; rho = (q/(\[Pi]  ndr^2 + 4  ndr  nrx) )/( Ck0);
Eqs = Inactivate[{D[u[t, x, y], 
      t] - (Laplacian[u[t, x, y], {x, y}] + 
       Div[(u[t, x, y])  Grad[w[t, x, y], {x, y}], {x, y}]), 
    D[v[t, x, y], t] - 
     d0 (Laplacian[v[t, x, y], {x, y}] - 
        Div[(v[t, x, y])  Grad[w[t, x, y], {x, y}], {x, y}]), 
    0.001  D[w[t, x, y], t] - Laplacian[w[t, x, y], {x, y}] - 
     4  \[Pi]  (u[t, x, y] - v[t, x, y] - 
        rho   Boole[{x, y} \[Element] nc])}, Laplacian | D];
ics = {u[0, x, y] == u0, v[0, x, y] == v0, w[0, x, y] == 0.0};
bcs = DirichletCondition[w[t, x, y] == V2, True];
sol = NDSolve[{Activate[Eqs] == {NeumannValue[0, True], 
      NeumannValue[0, True], 0}, ics, bcs}, {u, v, 
    w}, {x, y} \[Element] mesh, {t, 0, Tend}, 
   InterpolationOrder -> All];
Export[FileNameJoin[{"C:\\Efield\\Expt_20240214\\a1\\K.vs.Cl-1_11.20\\", "Test1.m"}], sol];

isol = Import[FileNameJoin[{"C:\\Efield\\Expt_20240214\\a1\\K.vs.Cl-1_11.20\\", 
"Test.m"}]];
Plot[Evaluate[{u[0.0001, x, 0], v[0.0001, x, 0]} /.First[isol]], {x, -1.5, 1.5}]
Evaluate[u[0.0001, 0, 0] /. First[isol]]

Evaluation the value of specific points also returned no results, but shows InterpolatingFunction like

![enter image description here

I checked the data file and it looks like

( Created with the Wolfram Language : www.wolfram.com ) {{u -> InterpolatingFunction[{{0., 0.0001}, {-1.4999999999999998, 1.5}, {-0.5, 0.5}}, {5, 4225, 5, {42, 8018, 0}, {4, 3, 3}, 0, 0, 0, 0, Indeterminate & , {}, {}, False}, {{0., 1.^-8, 2.^-8, 4.^-8, 8.^-8, 1.6^-7, 2.4000000000000003*^-7, 3.2*^-7, 4.800000000000001*^-7, 6.4*^-7, 8.^-7, 1.1199999999999999^-6, 1.4399999999999998*^-6, 2.0799999999999996*^-6, 2.6559999999999996*^-6, 3.2319999999999997*^-6, 3.808*^-6, 4.384*^-6, 5.536*^-6, 6.688*^-6, 8.09254755698846*^-6, .... *

What is wrong with my code?

Improve this question
$\endgroup$
16
  • $\begingroup$ Very likely, the problem is that you did not define Tend – you see, it is displayed in blue color. $\endgroup$
    – Domen
    Commented Feb 16 at 13:52
  • $\begingroup$ Tend was alreday defined and set to various values for trial. In the code, 0.0001 was set. $\endgroup$
    – dopey
    Commented Feb 16 at 13:55
  • $\begingroup$ I'm looking at the screenshot you provided. In there, Tend is not defined. Also, does Evaluate[v[0.0001, x, 0] /. First[isol]] return anything? $\endgroup$
    – Domen
    Commented Feb 16 at 14:11
  • 2
    $\begingroup$ @dopey, I strongly suggest you include the whole code for solving the equation and exporting to the question. Otherwise, it will be difficult for others to help and find the problem. For example, I've tried to make a simple example, similar to your, and my code works as expected. Also, try your code with a fresh kernel (by using Quit). $\endgroup$
    – Domen
    Commented Feb 16 at 15:18
  • 1
    $\begingroup$ Use the MX format, "Test1.mx".....There was another problem with InterpolatingFunction and the .m format here. It seems a different problem, but the .m format and InterpolatingFunction do not seem to work well together in both cases. Maybe the bugs are related. $\endgroup$
    – Goofy
    Commented Feb 18 at 17:27

2 Answers 2

Reset to default
2
$\begingroup$

As we know file with extension .m is a plain text file, therefore all data in it are available to use. For example we generate file Test1.m with code

Needs["NDSolve`FEM`"]
dr = 0.5 10^-0; rx = 2./2  10^-0; xs = -(dr + rx); ys = 0.0;
ec = RegionUnion[Disk[{-rx, 0}, dr], Rectangle[{-rx, -dr}, {rx, dr}], 
   Disk[{rx, 0}, dr]];
bmesh = ToBoundaryMesh[ec, "MaxBoundaryCellMeasure" -> .05];
bmesh["Wireframe"];
f = Function[{vertices, area}, 
   area > 0.002 (1. - 0.5 Norm[Mean[vertices]])];
(mesh = ToElementMesh[bmesh, MeshRefinementFunction -> f])["Wireframe"]
ndr = 0.3 10^-0; nrx = 
 1.2/2  10^-0; nxs = -(ndr + nrx); nys = 0.0; nlx = 0.4/1  10^-0;
nc = RegionUnion[Disk[{-nrx, 0}, ndr], 
   Rectangle[{-nrx, -ndr}, {(-nrx + nlx), ndr}], 
   Disk[{(-nrx + nlx), 0}, ndr], Disk[{(nrx - nlx), 0}, ndr], 
   Rectangle[{(nrx - nlx), -ndr}, {nrx, ndr}], Disk[{nrx, 0}, ndr]];
rPC = 0.00  10^-0;
xPC1 = rx + (dr - 1*rPC)*Cos[1 Pi/4]; yPC1 = (dr - 1*rPC)*
  Sin[1 Pi/4];
xPC2 = -rx - (dr - 1*rPC)*Cos[-Pi/3]; yPC2 = (dr - 1*rPC)*Sin[-Pi/3];
PtChg1 = Disk[{xPC1, yPC1}, rPC];
PtChg2 = Disk[{xPC2, yPC2}, rPC];
PtChg = RegionUnion[PtChg1, PtChg2];
Show[Region[Style[nc, Cyan]], mesh["Wireframe"], 
 Graphics[{Magenta, Thick, PointSize[0.01], Point[{xPC1, yPC1}], 
   Point[{xPC2, yPC2}]}]]
Ck0 = 139.20  10^9; Ccl0 = 12.43  10^9; 
u0 = Ck0/Ck0;
v0 = Ccl0/Ck0;
dk = 1.95  10^3; dcl = 2.02  10^3; 
d0 = dcl/dk;
Tend = 0.0001 ;
V2 = -0.05*37.436;  
q = 516.75  10^9; rho = (q/(\[Pi]  ndr^2 + 4  ndr  nrx) )/( Ck0);
Eqs = Inactivate[{D[u[t, x, y], 
      t] - (Laplacian[u[t, x, y], {x, y}] + 
       Div[(u[t, x, y])  Grad[w[t, x, y], {x, y}], {x, y}]), 
    D[v[t, x, y], t] - 
     d0 (Laplacian[v[t, x, y], {x, y}] - 
        Div[(v[t, x, y])  Grad[w[t, x, y], {x, y}], {x, y}]), 
    0.001  D[w[t, x, y], t] - Laplacian[w[t, x, y], {x, y}] - 
     4  \[Pi]  (u[t, x, y] - v[t, x, y] - 
        rho   Boole[{x, y} \[Element] nc])}, Laplacian | D];
ics = {u[0, x, y] == u0, v[0, x, y] == v0, w[0, x, y] == 0.0};
bcs = DirichletCondition[w[t, x, y] == V2, True];
sol = NDSolve[{Activate[Eqs] == {NeumannValue[0, True], 
      NeumannValue[0, True], 0}, ics, bcs}, {u, v, 
    w}, {x, y} \[Element] mesh, {t, 0, Tend}, 
   InterpolationOrder -> All];
Export[FileNameJoin[{"C:\\...\\", "Test1.m"}], sol];

Then we load this file with

isol = Import["C:\\...\\Test1.m"]

As output we have for example in v.12.2 Figure 1

And in v.14.0.0 it looks like Figure 2

Obviously in different versions InterpolatingFunction looks different. Nevertheless it is just data stored in it, and we can try to retrieve data as follows

{U, V, W} = {u, v, w} /. isol[[1]];

First we check time table

time = U[[3, 1]]

(*Out[]= {0., 1.*10^-8, 2.*10^-8, 4.*10^-8, 8.*10^-8, 1.6*10^-7, 
 2.4*10^-7, 3.2*10^-7, 4.8*10^-7, 6.4*10^-7, 8.*10^-7, 1.12*10^-6, 
 1.44*10^-6, 2.08*10^-6, 2.656*10^-6, 3.232*10^-6, 3.808*10^-6, 
 4.384*10^-6, 5.536*10^-6, 6.688*10^-6, 8.09255*10^-6, 
 9.4971*10^-6, 0.0000109016, 0.0000123062, 0.0000137107, \
0.0000151153, 0.0000179244, 0.0000207335, 0.0000235426, 0.0000263517, \
0.0000291608, 0.000034779, 0.0000403971, 0.0000460153, 0.0000516335, \
0.0000572517, 0.0000628699, 0.0000684881, 0.0000741063, 0.0000841063, \
0.0000941063, 0.0001}*)

Now we can prepare InterpolatingFunction for example at t=0.0001, or last element time[[-1]] as follows

u0001 = ElementMeshInterpolation[{mesh}, U[[4]][[-1, 1]]]

v0001 = ElementMeshInterpolation[{mesh}, V[[4]][[-1, 1]]]

w0001 = ElementMeshInterpolation[{mesh}, W[[4]][[-1, 1]]]

Finally plot data

{Plot3D[u0001[x, y], Element[{x, y}, mesh], ColorFunction -> Hue, 
  PlotPoints -> 50, PlotRange -> All], 
 Plot3D[v0001[x, y], Element[{x, y}, mesh], ColorFunction -> Hue, 
  PlotPoints -> 50, PlotRange -> All], 
 Plot3D[w0001[x, y], Element[{x, y}, mesh], ColorFunction -> Hue, 
  PlotPoints -> 50, PlotRange -> All]}

Figure 3

Please note, in different versions mesh generated in different ways, so that mesh from v.12.x.x is not same as in v.14.0.0. If we have no mesh available as above, then we can retrieve mesh elements from isol as follows

points=U[[3, 2]][[1]]

In v.14 we have output Figure 4

ListPlot[points]

Figure 5

Triangle elements

 triel=U[[3, 2]][[2]]

Figure6

Line elements

linel=U[[3, 2]][[3]]

Point Elements

pointel=U[[3, 2]][[4]]

In v.14 run first

Needs["NDSolve`FEM`"];

Then we can generate mesh as

mesh1= U[[3, 2]]

(*Out[]= ElementMesh[{{-1.5, 1.5}, {-0.5, 0.5}}, {TriangleElement[
   "<" 3917 ">"]}]*)

Then we can use mesh1 as above, for example for t=time[[-2]] we have

u0002 = ElementMeshInterpolation[{mesh1}, U[[4]][[-2, 1]]];

v0002 = ElementMeshInterpolation[{mesh1}, V[[4]][[-2, 1]]];

w0002 = ElementMeshInterpolation[{mesh1}, W[[4]][[-2, 1]]];

Visualization

{Plot3D[u0002[x, y], Element[{x, y}, mesh1], ColorFunction -> Hue, 
  PlotPoints -> 50, PlotRange -> All], 
 Plot3D[v0002[x, y], Element[{x, y}, mesh1], ColorFunction -> Hue, 
  PlotPoints -> 50, PlotRange -> All], 
 Plot3D[w0002[x, y], Element[{x, y}, mesh1], ColorFunction -> Hue, 
  PlotPoints -> 50, PlotRange -> All]}

Figure 7

Improve this answer
$\endgroup$
2
$\begingroup$

Export/import with ".m" files does not (it seems) preserve packed arrays. But InterpolatingFunction relies (it seems) on some of its data arrays being packed. I suggest that one consider using ".mx" files instead; they do not have this problem.

The following repairs the bug. (It seems to pack more arrays than are packed in the original solution, but it produces a working solution.)

jsol = isol /. 
   a_?(ArrayQ[#, _, NumericQ] &) :> Developer`ToPackedArray[a];

ByteCount /@ {isol, jsol}

(* {80590000, 27442504} *)

u[0.0001, 0, 0] /. First[jsol]

(* 1.00143 *)

Plot[Evaluate[{u[0.0001, x, 0], v[0.0001, x, 0]} /. 
   First[jsol]], {x, -1.5, 1.5}]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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