5
$\begingroup$

Here is my Mathematica code which implements the Newmark method to solve an equation of motion. The variable ag contains the acceleration values from an earthquake record. Is it possible to optimize this piece of code? It takes a long time to calculate a spectral response, which runs this code about 500 times, taking around 20 minutes.

The "Data.txt" can be found here or here.

Parallelize[Axel = Import["https://pastebin.com/raw/t7itJ7Q3"]];
Data = Array[{Axel[[#]][[1]], Axel[[#]][[2]]} &, Length[Axel] - 1, 2];
ag = 10 Data[[All, 2]];


M[y_] := N[
  Parallelize[
   u = fs = 
     kt = kthat = rhat = phat = v = a = ConstantArray[1., Length[ag]];
   k = 1.;
   m = 1.;
   p = -m ag;
   c = 0.1;
   fy = y;
   gamma = 1/2.; beta = 1/6.;
   v[[1]] = 0.;
   fs[[1]] = 0.;
   u[[1]] = 0.;
   a[[1]] = (p[[1]] - c v[[1]] - fs[[1]])/m;
   dt = 0.01;
   a1 = 1./(beta dt^2) m + gamma/(beta dt) c;
   a2 = 1./(beta dt) m + (gamma/beta - 1) c;
   a3 = (1/(2 beta) - 1.) m + dt (gamma/(2 beta) - 1) c;

   For[i = 1, i < Length[ag], i++,
    u[[i + 1]] = u[[i]];
    fs[[i + 1]] = fs[[i]];
    phat[[i + 1]] = p[[i + 1]] + a1 u[[i]] + a2 v[[i]] + a3 a[[i]];
    While[
     Abs[rhat[[i + 1]] = 
        phat[[i + 1]] - fs[[i + 1]] - a1 u[[i + 1]]] > 0.000001,
     u[[i + 1]] = u[[i + 1]] + rhat[[i + 1]]/(k + a1);
     If[
      (fz = fs[[i]] + k (u[[i + 1]] - u[[i]])) >= 0,
      fs[[i + 1]] = Min[fz, fy];,
      fs[[i + 1]] = Max[fz, -fy];
      ];
     ];
    v[[i + 1]] = 
     gamma/(beta dt) (u[[i + 1]] - u[[i]]) + (1. - gamma/
         beta) v[[i]] + dt (1 - gamma/(2 beta)) a[[i]];
    a[[i + 1]] = 
     1/(beta dt^2) (u[[i + 1]] - u[[i]]) - 
      1/(beta dt) v[[i]] - (1/(2 beta) - 1.) a[[i]];

    ];
   uy = fy/k;
   um = Max[Abs[u]];
   miu = Max[um/uy, 1]
   ]]

Here's a single run which takes quarter of a second, but implementing the function into a solver algorithm runs this around thousand times.

In[778]:= M[0.04343505938959662`] // AbsoluteTiming

Out[778]= {0.259419, 5.97369}

single solver run:

In[734]:= y = 6;
x = 1/y;
While[Abs[az = y - M[x]] > 0.05, 
   x = (1 - 0.1 Sign[az]) x]; // AbsoluteTiming
M[x]

Out[736]= {7.78515, Null}

Out[737]= 5.97369
$\endgroup$
  • 3
    $\begingroup$ There is a Newmark time integrator here and the documentation is here. Perhaps that's useful to you. $\endgroup$ – user21 May 3 '16 at 21:26
  • 1
    $\begingroup$ It might be good to add some example runs of this and some example data so that one can try it out; else one does not know if the optimization works correctly. $\endgroup$ – user21 May 3 '16 at 21:28
3
$\begingroup$

I have tried compiling my code in C, as well as storing the Length function into a variable. It has increased the speed by a factor of 100.

cf = Compile[{{x, _Real}}, Block[{u, fs, kt, kthat, rhat, phat, v, a,
    dt, k, c, m, gamma, beta, fy, a1, a2, a3, p, um, uy, i, fz, l},
   l = Length[ag]; 
   u = fs = kt = kthat = rhat = phat = v = a = ConstantArray[0., l];
   k = 1.;
   m = 1.;
   p = -m ag;
   c = 0.1;
   fy = x;
   gamma = 1/2.; beta = 1/6.;
   v[[1]] = 0.;
   fs[[1]] = 0.;
   u[[1]] = 0.;
   a[[1]] = (p[[1]] - c v[[1]] - fs[[1]])/m;
   dt = 0.01;
   a1 = 1./(beta dt^2) m + gamma/(beta dt) c;
   a2 = 1./(beta dt) m + (gamma/beta - 1) c;
   a3 = (1/(2 beta) - 1.) m + dt (gamma/(2 beta) - 1) c;

   For[i = 1, i < l, i++,
    u[[i + 1]] = u[[i]];
    fs[[i + 1]] = fs[[i]];
    phat[[i + 1]] = p[[i + 1]] + a1 u[[i]] + a2 v[[i]] + a3 a[[i]];
    While[
     Abs[rhat[[i + 1]] = 
        phat[[i + 1]] - fs[[i + 1]] - a1 u[[i + 1]]] > 0.000001,
     u[[i + 1]] = u[[i + 1]] + rhat[[i + 1]]/(k + a1);
     If[
      (fz = fs[[i]] + k (u[[i + 1]] - u[[i]])) >= 0,
      fs[[i + 1]] = Min[fz, fy];,
      fs[[i + 1]] = Max[fz, -fy];
      ];
     ];
    v[[i + 1]] = 
     gamma/(beta dt) (u[[i + 1]] - u[[i]]) + (1. - gamma/
         beta) v[[i]] + dt (1 - gamma/(2 beta)) a[[i]];
    a[[i + 1]] = 
     1/(beta dt^2) (u[[i + 1]] - u[[i]]) - 
      1/(beta dt) v[[i]] - (1/(2 beta) - 1.) a[[i]];

    ];
   uy = fy/k;
   um = Max[Abs[u]];
   miu = Max[um/uy, 1]
   ], CompilationTarget -> "C", Parallelization -> True]

Running times for the C-compiled version:

In[870]:= cf[0.0434] // AbsoluteTiming

Out[870]= {0.00213977, 5.9806}

In[887]:= y = 6;
x = 1/y;
While[Abs[az = y - cf[x]] > 0.05,x = (1 - 0.1 Sign[az]) x]; // AbsoluteTiming
cf[x]

Out[889]= {0.0474829, Null}

Out[890]= 5.97369

I'd appreciate it if you could point out any other possible optimizations ;)

| improve this answer | |
$\endgroup$

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.