Keep It Simple, Stupid PerlMonks

### Re: PDEs oo

by perlfan (Vicar)
 on Nov 03, 2021 at 06:34 UTC Need Help??

Do you want to share the code?

Here's the 2d version, in C/OpenMP and wrapped in Perl;

If you want to do this in pure Perl, this can be done rather easily. You're computing a "stencil" (pattern) over space, from the previous timestep. You iterate timesteps until the points in space cease to change by some difference (epsilon). You can also do this in PDL. Other relevant terms: Jacobi, Gauss-Seidel, SOR.

Replies are listed 'Best First'.
Re^2: PDEs oo
by Lucero (Novice) on Nov 03, 2021 at 16:34 UTC
```from numpy import *
import matplotlib.pylab as p
from mpl_toolkits.mplot3d import Axes3D

Nx = 101; Nt = 3000; Dx = 0.03 ; Dt = 0.9
KAPPA = 210. ; SPH = 900. ; RHO = 2700. #conductivity,specfheat ,densi
+ty
T = zeros( (Nx , 2) , float ); Tpl =zeros( (Nx , 31) , float )
print ("Working, wait for figure after count to 10")

for  ix in range ( 1 , Nx - 1) : T[ ix , 0] = 100.0 ; # initial temper
+ature
T [ 0 , 0 ] = 0.0 ; T [ 0 , 1 ] = 0. #first and last points a t0
T[Nx - 1 ,0] = 0. ; T[Nx - 1 ,1] = 0.0
cons = KAPPA / ( SPH*RHO) *Dt / ( Dx*Dx) ; # constant
m=1 # counter for rows , one every 300 time steps

for t in range ( 1 , Nt ) : # time iteration
for ix in range ( 1 , Nx - 1) : #Finitedifferences
T[ ix , 1] = T[ ix , 0] + cons *(T[ ix + 1 , 0] + T[ ix - 1, 0]
+ - 2.*T[ ix , 0 ] )
if t %300 == 0 or t == 1: #for t = 1 and every 300 time steps
for ix in range ( 1 , Nx - 1 , 2) : Tpl [ ix , m] = T[ ix , 1]
print (m)
m=m+1 # increasemevery 300 time steps
for ix in range ( 1 , Nx - 1): T[ ix , 0] = T[ ix , 1]# 100 p o s i
+ t o n s a t t =m
x = list ( range ( 1 , Nx - 1 , 2) ) # p l o t every o t h e r x p o i
+ n t
y = list ( range ( 1 , 30) ) # every 10 p o i n t s i n y ( time )
X, Y = p.meshgrid ( x , y ) # g r i d f o r p o s i t i o n and time

def functz ( Tpl ) : # Function r e t u r n s t e m p e r a t u r e
z = Tpl [X, Y]
return z

Z = functz ( Tpl )
fig = p.figure( ) # c r e a t e f i g u r e
ax = Axes3D ( fig ) # p l o t s a x i s
ax.plot_wireframe(X, Y, Z , color = 'r') # red wireframe
ax.set_xlabel('Position') # l a b e l axes
ax.set_ylabel('time')
ax.set_zlabel('Temperature')
p.show( ) # shows f i g u r e , c l o s e Python s h e l l
print ("finished")
Switching to Perl, I can change the matrix elements but I dont know how to change the dimensions of the matrix.

Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: note [id://11138362]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others meditating upon the Monastery: (6)
As of 2022-08-18 11:04 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?

No recent polls found

Notices?