1
$\begingroup$

I was having a look at the following Wolfram demonstration that is pertinent to my current work.

Demonstration: 2D Lid-driven Stokes flow

(And here is a quick report that outlines the problem and solution method: https://www.researchgate.net/publication/341090559_2-D_Stokes_Solution_for_Lid-Driven_Cavity_Flow)

It plots the streamlines for a certain problem in fluid dynamics, e.g

plot

What I would like to do is get the raw data that is used to generate this plot. However, if you look at the source code used to make this visualization, all of the relevant data is wrapped inside a Manipulate[Module ... ] ] environment, so they are all local variables. When I tried to extract the code out of the module, it didn't work (Ran it for 10 minutes and it didn't terminate whereas the original takes a few seconds.)

I would like to be able to get the data stored inside the sol array using Nx=101, Ny=101, and gamma=1, and export it to some useful format, like .csv. How can I do that?

Note: I know my question is similar to this one, but I don't just want the data used to make the plot. I want the full array named sol.

Source code is below.


Manipulate[
 Module[{a, b, Nx, Ny, \[Psi], \[Omega], StreamfuncEqn, VorticityEqn, 
   SFBC1, SFBC2, SFBC3, SFBC4, VortBC1, VortBC2, VortBC3, VortBC4, 
   Eqns, systemEqns, var, SFData1, \[Gamma], h, sol, A, b1, myarrow, 
   plt1, Nymax, Nxmax, contours},
  Nxmax = 101; Nymax = 101;
  a = 1; b = 1; Nx = Nxmax; Ny = Nymax; \[Gamma] = gamma;
  h = a/(Nx - 1) // N;
  
  contours = 
   Which[gamma == 1, {-0.0015, -0.01, -0.03, -0.05, -0.09, -5 10^-7 , 
     10^-7 , 10^-6, 5 10^-6}, 
    gamma == 2, {-0.0015, -0.02, -0.05, -0.07, -5 10^-7, 10^-6, 
     10^-11},
    gamma == 
     3, {-0.0015, -0.0075, -0.02, -0.03698, -0.045, -0.049, -5 10^-7, 
     10^-6, 10^-11},
    gamma == 
     4, {-0.0015, -0.0075, -0.02, -0.03, -0.03698, -0.045, -0.049, -5 \
10^-7, 10^-6, 10^-11}, 
    gamma == 0.75, {-0.0010, -0.0075, -0.05, -0.09, -5 10^-7, -10^-6, 
     10^-6, 10^-11},
    gamma == 0.5, {-0.001, -0.005, -0.02, -0.05, -0.08, -5 10^-7, 
     10^-6, 10^-4, 2 10^-4, 2 10^-5, -2 10^-9}, gamma == 0.25, 
    Sort[{-0.05, -0.005, -0.0015, -0.0002, -5 10^-7, -10^-8, -10^-9, \
- 10^-7, 10^-6, 10^-4, 10^-5, 2 10^-4, 10^-9, 3.5 10^-9, 7 10^-9, 
      1.5 10^-9, 10^-10}], gamma == 0.20, 
    Sort[{-0.05, -0.005, -0.0015, -0.0002, -5 10^-7, -10^-8, -10^-9, \
- 10^-7, 10^-6, 10^-4, 10^-5, 2 10^-4, 10^-9, 3.5 10^-9, 7 10^-9, 
      1.5 10^-9, 10^-10}]]; 
  StreamfuncEqn[i_, 
    j_] := (\[Psi][i + 1, j] + \[Psi][i - 1, j] - 
      2 \[Psi][i, 
        j] + \[Gamma]^2 (\[Psi][i, j - 1] + \[Psi][i, j + 1] - 
         2 \[Psi][i, j]) + h^2 \[Omega][i, j]) == 0;
  VorticityEqn[i_, 
    j_] := (\[Omega][i + 1, j] + \[Omega][i - 1, j] - 
      2 \[Omega][i, 
        j] + \[Gamma]^2 (\[Omega][i, j - 1] + \[Omega][i, j + 1] - 
         2 \[Omega][i, j])) == 0;
  SFBC1 = Table[\[Psi][i, 1] == 0, {i, 1, Nx}];
  SFBC2 = Table[\[Psi][Nx, j] == 0, {j, 2, Ny - 1}];
  SFBC3 = Table[\[Psi][i, Ny] == 0, {i, 1, Nx}];
  SFBC4 = Table[\[Psi][1, j] == 0, {j, 2, Ny - 1}];
  VortBC1 = 
   Table[\[Omega][1, j] == -2 \[Psi][2, j]/h^2, {j, 2, Ny - 1}];
  VortBC2 = 
   Table[\[Omega][i, 1] == -2 \[Gamma]^2 \[Psi][i, 2]/h^2, {i, 1, Nx}];
  VortBC3 = 
   Table[\[Omega][Nx, j] == -2 \[Psi][Nx - 1, j]/h^2, {j, 2, Ny - 1}];
  VortBC4 = 
   Table[\[Omega][i, 
      Ny] == -2 \[Gamma]^2 (\[Psi][i, Ny - 1] + h/\[Gamma])/h^2, {i, 
     1, Nx}];
  Eqns = Table[{StreamfuncEqn[i, j], VorticityEqn[i, j]}, {i, 2, 
      Nx - 1}, {j, 2, Ny - 1}] // Flatten;
  systemEqns = 
   Join[Eqns, SFBC1, SFBC2, SFBC3, SFBC4, VortBC1, VortBC2, VortBC3, 
    VortBC4];
  var = Union[
    Cases[systemEqns, \[Psi][_, _] | \[Omega][_, _], \[Infinity]]];
  {b1, A} = CoefficientArrays[systemEqns, var];
  sol = LinearSolve[A, -b1];
  SFData1 = Transpose[Partition[Take[sol, Length[var]/2], Ny]];
  plt1 = ListContourPlot[SFData1, Contours -> contours, 
    ColorFunction -> "Pastel", ContourShading -> Automatic, 
    ContourStyle -> Black, DataRange -> {{0, 1}, {0, 1/\[Gamma]}}, 
    FrameLabel -> {Style[x, 16], Style[y, 16]}, 
    AspectRatio -> Automatic, MaxPlotPoints -> Infinity, 
    InterpolationOrder -> 2, PlotRange -> {Full, Full, All}, 
    Frame -> True, ImageSize -> 400 {1, 1}, 
    Epilog -> {{Transparent, EdgeForm[{Thickness[0.007], Black}], 
       Rectangle[{0, 0}, {1, 1/gamma}]}}]],
 {{gamma, 1, "cavity aspect ratio (W/H) ="}, {0.20, 0.25, 0.5, 0.75, 
   1, 2, 3, 4}, ControlType -> Setter},
 TrackedSymbols -> {gamma}, ControlPlacement -> Top,
 ContentSize -> {450, 450}]
Improve this question
$\endgroup$
2

1 Answer 1

Reset to default
1
$\begingroup$

You can remove the Manipulate and the code to set gamma at the bottom, then just set gamma to some desired value within the Module. However, examining the Which, each value of gamma is required to be '==' to a specific number, so setting it to 101 won't work. There is no default clause, so unless you modify the code, you need to choose from one of the values in the code: gamma = 1,2,3,4,0.75,0.5,0.25,0.20 are the only supported values in the demo.I chose to set it to '1' here instead of 101. Then just add a semicolon after the plot, and type sol to dump the values.

Module[{a, b, Nx, Ny, \[Psi], \[Omega], StreamfuncEqn, VorticityEqn, 
  SFBC1, SFBC2, SFBC3, SFBC4, VortBC1, VortBC2, VortBC3, VortBC4, 
  Eqns, systemEqns, var, SFData1, \[Gamma], h, sol, A, b1, myarrow, 
  plt1, Nymax, Nxmax, contours},
 Nxmax = 101; Nymax = 101;
 a = 1; b = 1; Nx = Nxmax; Ny = Nymax; \[Gamma] = gamma = 1;
 h = a/(Nx - 1) // N;
 contours = Which[
   gamma == 1, {-0.0015, -0.01, -0.03, -0.05, -0.09, -5 10^-7, 10^-7, 
    10^-6, 5 10^-6},
   gamma == 2, {-0.0015, -0.02, -0.05, -0.07, -5 10^-7, 10^-6, 
    10^-11}, 
   gamma == 
    3, {-0.0015, -0.0075, -0.02, -0.03698, -0.045, -0.049, -5 10^-7, 
    10^-6, 10^-11},
   gamma == 
    4, {-0.0015, -0.0075, -0.02, -0.03, -0.03698, -0.045, -0.049, -5 \
10^-7, 10^-6, 10^-11},
   gamma == 0.75, {-0.0010, -0.0075, -0.05, -0.09, -5 10^-7, -10^-6, 
    10^-6, 10^-11}, 
   gamma == 0.5, {-0.001, -0.005, -0.02, -0.05, -0.08, -5 10^-7, 
    10^-6, 10^-4, 2 10^-4, 2 10^-5, -2 10^-9},
   gamma == 0.25, 
   Sort[{-0.05, -0.005, -0.0015, -0.0002, -5 10^-7, -10^-8, -10^-9, \
-10^-7, 10^-6, 10^-4, 10^-5, 2 10^-4, 10^-9, 3.5 10^-9, 7 10^-9, 
     1.5 10^-9, 10^-10}],
   gamma == 0.20, 
   Sort[{-0.05, -0.005, -0.0015, -0.0002, -5 10^-7, -10^-8, -10^-9, \
-10^-7, 10^-6, 10^-4, 10^-5, 2 10^-4, 10^-9, 3.5 10^-9, 7 10^-9, 
     1.5 10^-9, 10^-10}]];
 StreamfuncEqn[i_, 
   j_] := (\[Psi][i + 1, j] + \[Psi][i - 1, j] - 
     2 \[Psi][i, 
       j] + \[Gamma]^2 (\[Psi][i, j - 1] + \[Psi][i, j + 1] - 
        2 \[Psi][i, j]) + h^2 \[Omega][i, j]) == 0;
 VorticityEqn[i_, 
   j_] := (\[Omega][i + 1, j] + \[Omega][i - 1, j] - 
     2 \[Omega][i, 
       j] + \[Gamma]^2 (\[Omega][i, j - 1] + \[Omega][i, j + 1] - 
        2 \[Omega][i, j])) == 0;
 SFBC1 = Table[\[Psi][i, 1] == 0, {i, 1, Nx}];
 SFBC2 = Table[\[Psi][Nx, j] == 0, {j, 2, Ny - 1}];
 SFBC3 = Table[\[Psi][i, Ny] == 0, {i, 1, Nx}];
 SFBC4 = Table[\[Psi][1, j] == 0, {j, 2, Ny - 1}];
 VortBC1 = 
  Table[\[Omega][1, j] == -2 \[Psi][2, j]/h^2, {j, 2, Ny - 1}];
 VortBC2 = 
  Table[\[Omega][i, 1] == -2 \[Gamma]^2 \[Psi][i, 2]/h^2, {i, 1, 
    Nx}];
 VortBC3 = 
  Table[\[Omega][Nx, j] == -2 \[Psi][Nx - 1, j]/h^2, {j, 2, Ny - 1}];
 VortBC4 = 
  Table[\[Omega][i, 
     Ny] == -2 \[Gamma]^2 (\[Psi][i, Ny - 1] + h/\[Gamma])/h^2, {i, 1,
     Nx}];
 Eqns = Table[{StreamfuncEqn[i, j], VorticityEqn[i, j]}, {i, 2, 
     Nx - 1}, {j, 2, Ny - 1}] // Flatten;
 systemEqns = 
  Join[Eqns, SFBC1, SFBC2, SFBC3, SFBC4, VortBC1, VortBC2, VortBC3, 
   VortBC4];
 var = Union[
   Cases[systemEqns, \[Psi][_, _] | \[Omega][_, _], \[Infinity]]];
 {b1, A} = CoefficientArrays[systemEqns, var];
 sol = LinearSolve[A, -b1];
 SFData1 = Transpose[Partition[Take[sol, Length[var]/2], Ny]];
 plt1 = ListContourPlot[SFData1, Contours -> contours, 
   ColorFunction -> "Pastel", ContourShading -> Automatic, 
   ContourStyle -> Black, DataRange -> {{0, 1}, {0, 1/\[Gamma]}}, 
   FrameLabel -> {Style[x, 16], Style[y, 16]}, 
   AspectRatio -> Automatic, MaxPlotPoints -> Infinity, 
   InterpolationOrder -> 2, PlotRange -> {Full, Full, All}, 
   Frame -> True, ImageSize -> 400 {1, 1}, 
   Epilog -> {{Transparent, EdgeForm[{Thickness[0.007], Black}], 
      Rectangle[{0, 0}, {1, 1/gamma}]}}]
 ; sol
 ]

And once you verify that it works, just replace sol at the bottom with Export[sol,].

Improve this answer
$\endgroup$
3
  • $\begingroup$ Whoops, my mistake. I want gamma=1, not 101. $\endgroup$
    – K.defaoite
    Commented Jan 28, 2023 at 19:04
  • $\begingroup$ Good, that's what I set it to. $\endgroup$
    – user87932
    Commented Jan 28, 2023 at 19:15
  • $\begingroup$ Turns out what I wanted is actually the SFData1 object. Anyway, your method worked a treat! Thanks. $\endgroup$
    – K.defaoite
    Commented Jan 28, 2023 at 20:35

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.