Failed to retrieve raw data from importing InterpolatingFunction derived from NDSolve

Question

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

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?

Very likely, the problem is that you did not define Tend – you see, it is displayed in blue color. — Domen, Feb 16 '24 at 13:52
Tend was alreday defined and set to various values for trial. In the code, 0.0001 was set. — dopey, Feb 16 '24 at 13:55
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? — Domen, Feb 16 '24 at 14:11
I revised the screenshot as shown above. Evaluate[v[0.0001, x, 0] /. First[isol]] also returned the same as shown above. Thanks for your response~ — dopey, Feb 16 '24 at 14:20
We don't know what you export exactly, but maybe first you need to define interpolation function like U=u /. First[isol]] and then evaluate U[.001,0,0]. — Alex Trounev, Feb 16 '24 at 14:34
@Alex, I also tried to define interpolation function like you suggested and returned the same result. What I did is like "sol=NDSolve[....]; Export["Test.m",sol];". — dopey, Feb 16 '24 at 14:44
@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). — Domen, Feb 16 '24 at 15:18
@Domen, Thanks for your help and reediting my question. I tried your example, it works well. However, it's not the case as running my code. — dopey, Feb 16 '24 at 17:14
@dopey Your code generated file Test1.m of about 76.9Mb. What actually do you try to simulate? — Alex Trounev, Feb 17 '24 at 15:38
@ Alex, The code is to simulate the initial charging process of a battery with 2 electrodes located to the cell peripheral and central region. u and v are concentrations for cations and anions respectively and w is electrical potential. The Tend here is only set to 0.0001 for demonstration, 5 is used instead in real simulation. The output size of data file is about 300-400mb. Several parametric condidtions were tuned through ParallelDo. — dopey, Feb 18 '24 at 07:01
@dopey It looks like InterpolatingFunction retrieved from the text file Test1.m is some data available to use. If you like I can show how to use it. — Alex Trounev, Feb 18 '24 at 11:58
@Alex, yes, Please show how to use it. I have to store data in files because I got other simulations with even much longer comuptations. Thanks a lot~ — dopey, Feb 18 '24 at 13:01
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. — Goofy, Feb 18 '24 at 17:27
As a side question: Why do you Inactivate your PDE only to Activate it later again? — user21, Feb 18 '24 at 20:12
This looks a bit like a bug, you should send it to WRI: support@wolfram.com — user21, Feb 18 '24 at 20:13
@Goofy, it's unlikely that the bugs are related, the one is arbitrary precision code this one is machine precision. — user21, Feb 18 '24 at 20:15

Alex Trounev · Accepted Answer · 2024-02-18T17:01:30.410

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

And in v.14.0.0 it looks like

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.610^-7, 
 2.410^-7, 3.210^-7, 4.810^-7, 6.410^-7, 8.10^-7, 1.1210^-6, 
 1.4410^-6, 2.0810^-6, 2.65610^-6, 3.23210^-6, 3.80810^-6, 
 4.38410^-6, 5.53610^-6, 6.68810^-6, 8.0925510^-6, 
 9.497110^-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]}

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

ListPlot[points]

Triangle elements

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

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]}

score 2 · Answer 2 · answered Feb 18 '24 at 17:55

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}]

Failed to retrieve raw data from importing InterpolatingFunction derived from NDSolve

2 Answers2