XP is just a number 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 wandering the Monastery: (3)
As of 2022-01-19 01:50 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
In 2022, my preferred method to securely store passwords is:

Results (54 votes). Check out past polls.

Notices?