Parallel Spectral Numerical Methods/Incompressible Magnetohydrodynamics

Incompressible Magnetohydrodynamics Background The equations of incompressible magnetohydrodynamics are similar to the Navier-Stokes equations, but can show a variety of other phenomena. This section assumes that the reader is familiar with the material covered for the Navier-Stokes case here. Their mathematical investigation of the incompressible magnetohydrodynamics equations is also in as poor a state as that of the Navier Stokes equations. A physical introduction to these equations can be found in Carbone and PouquetCarbone and Pouquet (2009), ChandrasekharChandrasekhar (1961), DavidsonDavidson (2001), Gerbeau, Le Bris and LelièvreGerbeau, Le Bris and Lelièvre (2006), SchnakSchnak (2009) and VermaVerma (2004). A more general overview can be found here Magnetohydrodynamics. Under the assumption that fluid velocities are significantly smaller than the speed of light, and the fluid is nearly incompressible, dimensionless equations are given as : 1 : 2 : . 3 : . 4 In these equations, is the fluid velocity with components in the , and directions, is the magnetic field with components in the , and directions, is pressure field. Equation 1 represents conservation of momentum, eq. 2 corresponds to Maxwell's equations for nonrelativistic flows for which the time evolution of the electric field is neglected (see VermaVerma (2004)), eq. 3 is the continuity equation which represents conservation of mass for an incompressible fluid and eq. 4 represents the fact that there are no magnetic monopole sources. The Two Dimensional Case We will first consider the two-dimensional case. A difficulty in simulating the incompressible magnetohydrodynamic equations is the numerical satisfaction of the constraints in eqs. 3 and 4 , these are sometimes referred to as a divergence free conditions or a solenoidal constraints. To automatically satisfy these constraints in two dimensions, where and it is possible to re-write the equations using a potential formulation. In this case, we let and where is the streamfunction and the magnetic potential. Note that , and , so eqs. 3 and 4 are automatically satisfied. Making this change of variables, we obtain a reduced system of two coupled partial differential equation by taking the curl of the momentum equation, eq. 1 and the magnetic transport eq. 2 . We define the fluid vorticity , so that and the magnetic vorticity so that and eqs. 1 and 2 become : 5 : 6 : . 7 and : . 8 where and represent the and components of the force . The Three Dimensional Case In the three dimensional case, a potential formulation of the equations is not as convenient. Consequently, it is convenient to introduce the notion of a magnetic pressure to allow for a simple algorithm for enforcing the divergence free magnetic field constraint. The equations then become, : , 9 and : , 10 both of which can be discretized by the implicit midpoint rule. Here denotes the fluids Reynolds number and the magnetic Reynolds number. The fluid pressure is given by and the magnetic ``pressure" is given by Serial Programs We first write Matlab programs to demonstrate how to solve these equations on a single processor. The first program uses implicit midpoint rule timestepping to solve the two-dimensional magnetohydrodynamic equations and is in listing A . In this program, an additional equation is solved for a passive scalar with a diffusion constant : 11 A A Matlab program which finds a numerical solution to the 2D magnetohydrodynamic equations Code download % Numerical solution of the 2D incompressible Magnetohydrodynamics equations % on a Square Domain [0,1]x[0,1] using a Fourier pseudo-spectral method % and Implicit midpoint rule timestepping. A passive scalar is also advected % by the flow % %Periodic free-slip boundary conditions and Initial conditions: %u(x,y,0)=sin(2*pi*x)cos(2*pi*y) %v(x,y,0)=-cos(2*pi*x)sin(2*pi*y) %Analytical Solution: %u(x,y,t)=sin(2*pi*x)cos(2*pi*y)exp(-8*pi^2*t/Re) %v(x,y,t)=-cos(2*pi*x)sin(2*pi*y)exp(-8*pi^2*t/Re) % also consider an advection field % \theta_t+u\theta_x+v\theta_y=Dtheta(\theta_{xx}+\theta_{yy}) clear all; format compact; format short; clc; clf; Re=100;%Reynolds number Rem=100; % Magnetic Reynolds number Dtheta=0.01;%scalar diffusion constant %grid Nx=64; h=1/Nx; x=h*(1:Nx); Ny=64; h=1/Ny; y=h*(1:Ny)'; [xx,yy]=meshgrid(x,y); %initial conditions for velocity field u=sin(2*pi*xx).*cos(2*pi*yy); v=-cos(2*pi*xx).*sin(2*pi*yy); u_y=-2*pi*sin(2*pi*xx).*sin(2*pi*yy); v_x=2*pi*sin(2*pi*xx).*sin(2*pi*yy); omega=v_x-u_y; % initial magnetic potential field alpha=sin(4*pi*xx).*cos(6*pi*yy); % initial passive scalar field theta=exp(-2.0*((4*pi*(xx-0.25)).^2+(4*pi*(yy-0.25)).^2)); dt=0.0025; t(1)=0; tmax=1.0; nplots=ceil(tmax/dt); %wave numbers for derivatives k_x=2*pi*(1i*[(0:Nx/2-1) 0 1-Nx/2:-1]'); k_y=2*pi*(1i*[(0:Ny/2-1) 0 1-Ny/2:-1]); k2x=k_x.^2; k2y=k_y.^2; %wave number grid for multiplying matricies [kxx,kyy]=meshgrid(k2x,k2y); [kx,ky]=meshgrid(k_x,k_y); % use a high tolerance so time stepping errors % are not dominated by errors in solution to nonlinear % system tol=10^(-10); %compute \hat{Phi}^{n+1,k+1} alphahat=fft2(alpha); phihat=-alphahat./(kxx+kyy); %NOTE: kxx+kyy has to be zero at the following points to avoid a % discontinuity. However, we suppose that the streamfunction has % mean value zero, so we set them equal to zero phihat(1,1)=0; phihat(Nx/2+1,Ny/2+1)=0; phihat(Nx/2+1,1)=0; phihat(1,Ny/2+1)=0; %computes {\psi}_x by differentiation via FFT dphix = real(ifft2(phihat.*kx)); %computes {\psi}_y by differentiation via FFT dphiy = real(ifft2(phihat.*ky)); % components of magnetic field bx=dphiy; by=-dphix; %compute \hat{\omega}^{n+1,k} omegahat=fft2(omega); alphahat=fft2(alpha); thetatotal(1)=sum(sum(theta.^2))*h*h; for i=1:nplots chg=1; % save old values uold=u; vold=v; omegaold=omega; omegacheck=omega; omegahatold=omegahat; thetaold=theta; thetacheck=theta; alphahatold=alphahat; alphaold=alpha; alphacheck=alpha; bxold=bx; byold=by; while chg>tol % Fluid field %nonlinear {n+1,k} nonlinhat=0.25*fft2((u+uold).*ifft2((omegahat+omegahatold).*kx)+... (v+vold).*ifft2((omegahat+omegahatold).*ky)-... (bx+bxold).*ifft2((alphahat+alphahatold).*kx)-... (by+byold).*ifft2((alphahat+alphahatold).*ky)); %Implicit midpoint rule timestepping omegahat=((1/dt + 0.5*(1/Re)*(kxx+kyy)).*omegahatold-nonlinhat)... ./(1/dt -0.5*(1/Re)*(kxx+kyy)); %compute \hat{\psi}^{n+1,k+1} psihat=-omegahat./(kxx+kyy); %NOTE: kxx+kyy has to be zero at the following points to avoid a % discontinuity. However, we suppose that the streamfunction has % mean value zero, so we set them equal to zero psihat(1,1)=0; psihat(Nx/2+1,Ny/2+1)=0; psihat(Nx/2+1,1)=0; psihat(1,Ny/2+1)=0; %computes {\psi}_x by differentiation via FFT dpsix = real(ifft2(psihat.*kx)); %computes {\psi}_y by differentiation via FFT dpsiy = real(ifft2(psihat.*ky)); u=dpsiy; %u^{n+1,k+1} v=-dpsix; %v^{n+1,k+1} % magnetic field nonlinhat=0.25*fft2((u+uold).*ifft2((alphahat+alphahatold).*kx)+... (v+vold).*ifft2((alphahat+alphahatold).*ky)-... (bx+bxold).*ifft2((omegahat+omegahatold).*kx)-... (by+byold).*ifft2((omegahat+omegahatold).*ky)); %Implicit midpoint rule timestepping alphahat=((1/dt + 0.5*(1/Re)*(kxx+kyy)).*alphahatold-nonlinhat)... ./(1/dt -0.5*(1/Rem)*(kxx+kyy)); %compute \hat{\psi}^{n+1,k+1} phihat=-alphahat./(kxx+kyy); %NOTE: kxx+kyy has to be zero at the following points to avoid a % discontinuity. However, we suppose that the streamfunction has % mean value zero, so we set them equal to zero phihat(1,1)=0; phihat(Nx/2+1,Ny/2+1)=0; phihat(Nx/2+1,1)=0; phihat(1,Ny/2+1)=0; %computes {\psi}_x by differentiation via FFT dphix = real(ifft2(phihat.*kx)); %computes {\psi}_y by differentiation via FFT dphiy = real(ifft2(phihat.*ky)); bx=dphiy; %u^{n+1,k+1} by=-dphix; %v^{n+1,k+1} thetax=0.5*ifft2(kx.*fft2(theta+thetaold)); thetay=0.5*ifft2(ky.*fft2(theta+thetaold)); theta=ifft2((fft2(thetaold-dt*0.5*((uold+u).*thetax+(vold+v).*thetay))+... dt*0.5*Dtheta*(kxx+kyy).*fft2(thetaold))./(1-dt*0.5*Dtheta*(kxx+kyy))); %\omega^{n+1,k+1} omega=ifft2(omegahat); % check for convergence chg=max(max(abs(omega-omegacheck))) + max(max(abs(theta-thetacheck)))+... max(max(abs(alpha-alphacheck))) % store omega and theta to check for convergence of next iteration omegacheck=omega; thetacheck=theta; alphacheck=alpha; end t(i+1)=t(i)+dt; thetatotal(i+1)=sum(sum(theta.^2))*h*h; figure(1); subplot(2,2,1); pcolor(xx,yy,omega); shading interp; xlabel x; ylabel y; title(['Fluid Vorticity, Time ',num2str(t(i+1))]); colorbar; drawnow; subplot(2,2,2); pcolor(xx,yy,alpha); shading interp; xlabel x; ylabel y; title(['Alpha, Time ',num2str(t(i+1))]); colorbar; drawnow; subplot(2,2,3); pcolor(xx,yy,real(ifft2(psihat))); shading interp; xlabel x; ylabel y; title(['Psi, Time ',num2str(t(i+1))]); colorbar; drawnow; subplot(2,2,4); pcolor(xx,yy,real(theta)); shading interp; xlabel x; ylabel y; title(['Theta, Time ',num2str(t(i+1))]); colorbar; drawnow; end figure(2); clf; subplot(2,1,1); plot(t,thetatotal,'r-'); xlabel time; ylabel('Total theta'); subplot(2,1,2); semilogy(t,abs(1-thetatotal./thetatotal(1)),'r-'); xlabel time; ylabel('Total theta Error'); A1 A Python program which finds a numerical solution to the 2D magnetohydrodynamic equations Code download #!/usr/bin/env python """ Numerical solution of the 2D incompressible Magnetohydrodynamics equations on a Square Domain [0,1]x[0,1] using a Fourier pseudo-spectral method and Implicit Midpoint rule timestepping. A scalar field is also advected by the flow \theta_t+u\theta_x+v\theta_y=Dtheta(\theta_{xx}+\theta_{yy}) Periodic free-slip boundary conditions and Initial conditions: u(x,y,0)=sin(2*pi*x)cos(2*pi*y) v(x,y,0)=-cos(2*pi*x)sin(2*pi*y) """ import math import numpy import matplotlib.pyplot as plt from mayavi import mlab import time # Grid N=64; h=1.0/N x = [h*i for i in xrange(1,N+1)] y = [h*i for i in xrange(1,N+1)] numpy.savetxt('x.txt',x) xx=numpy.zeros((N,N), dtype=float) yy=numpy.zeros((N,N), dtype=float) for i in xrange(N): for j in xrange(N): xx[i,j] = x[i] yy[i,j] = y[j] dt=0.0025; t=0.0; tmax=1.0 #nplots=int(tmax/dt) Re=100 Rem=100 Dtheta=0.001 u=numpy.zeros((N,N), dtype=float) v=numpy.zeros((N,N), dtype=float) u_y=numpy.zeros((N,N), dtype=float) v_x=numpy.zeros((N,N), dtype=float) omega=numpy.zeros((N,N), dtype=float) theta=numpy.zeros((N,N), dtype=float) thetaold=numpy.zeros((N,N), dtype=float) thetax=numpy.zeros((N,N), dtype=float) thetay=numpy.zeros((N,N), dtype=float) thetacheck=numpy.zeros((N,N), dtype=float) alpha=numpy.zeros((N,N), dtype=float) alphaold=numpy.zeros((N,N), dtype=float) bx=numpy.zeros((N,N), dtype=float) bxold=numpy.zeros((N,N), dtype=float) by=numpy.zeros((N,N), dtype=float) byold=numpy.zeros((N,N), dtype=float) # Initial conditions for i in range(len(x)): for j in range(len(y)): u[i][j]=numpy.sin(2*math.pi*x[i])*numpy.cos(2*math.pi*y[j]) v[i][j]=-numpy.cos(2*math.pi*x[i])*numpy.sin(2*math.pi*y[j]) u_y[i][j]=-2*math.pi*numpy.sin(2*math.pi*x[i])*numpy.sin(2*math.pi*y[j]) v_x[i][j]=2*math.pi*numpy.sin(2*math.pi*x[i])*numpy.sin(2*math.pi*y[j]) omega[i][j]=v_x[i][j]-u_y[i][j] theta[i][j]=numpy.exp(-2.0*((4*math.pi*(x[i]-0.3))**2+(4*math.pi*(y[j]-0.3))**2)) alpha[i][j]=numpy.sin(4*math.pi*x[i])*numpy.cos(6*math.pi*y[j]) src = mlab.imshow(xx,yy,theta,colormap='jet') mlab.scalarbar(object=src) mlab.xlabel('x',object=src) mlab.ylabel('y',object=src) # Wavenumber k_x = 2*math.pi*numpy.array([complex(0,1)*n for n in range(0,N/2) \ + [0] + range(-N/2+1,0)]) k_y=k_x kx=numpy.zeros((N,N), dtype=complex) ky=numpy.zeros((N,N), dtype=complex) kxx=numpy.zeros((N,N), dtype=complex) kyy=numpy.zeros((N,N), dtype=complex) for i in xrange(N): for j in xrange(N): kx[i,j] = k_x[i] ky[i,j] = k_y[j] kxx[i,j] = k_x[i]**2 kyy[i,j] = k_y[j]**2 tol=10**(-10) psihat=numpy.zeros((N,N), dtype=complex) phihat=numpy.zeros((N,N), dtype=complex) alphahat=numpy.zeros((N,N), dtype=complex) alphahatold=numpy.zeros((N,N), dtype=complex) omegahat=numpy.zeros((N,N), dtype=complex) omegahatold=numpy.zeros((N,N), dtype=complex) thetahat=numpy.zeros((N,N), dtype=complex) thetahatold=numpy.zeros((N,N), dtype=complex) nlhat=numpy.zeros((N,N), dtype=complex) dpsix=numpy.zeros((N,N), dtype=float) dpsiy=numpy.zeros((N,N), dtype=float) dphix=numpy.zeros((N,N), dtype=float) dphiy=numpy.zeros((N,N), dtype=float) omegacheck=numpy.zeros((N,N), dtype=float) omegaold=numpy.zeros((N,N), dtype=float) temp=numpy.zeros((N,N), dtype=float) omegahat=numpy.fft.fft2(omega) thetahat=numpy.fft.fft2(theta) alphahat=numpy.fft.fft2(alpha) while (t<=tmax): chg=1.0 # Save old values uold=u vold=v omegaold=omega omegacheck=omega omegahatold = omegahat thetaold=theta thetahatold=thetahat thetacheck=theta bxold=bx byold=by alphahatold=alphahat alphaold=alpha alphacheck=alpha while(chg>tol): # fluid field # nolinear {n+1,k} nlhat=0.25*numpy.fft.fft2((u+uold)*numpy.fft.ifft2((omegahat+omegahatold)*kx)+\ (v+vold)*numpy.fft.ifft2((omegahat+omegahatold)*ky)-\ (bx+bxold)*numpy.fft.ifft2((alphahat+alphahatold)*kx)-\ (by+byold)*numpy.fft.ifft2((alphahat+alphahatold)*ky)) # Implicit midpoint rule timestepping omegahat=((1/dt + 0.5*(1/Re)*(kxx+kyy))*omegahatold \ -nlhat) \ /(1/dt -0.5*(1/Re)*(kxx+kyy)) psihat=-omegahat/(kxx+kyy) psihat[0][0]=0 psihat[N/2][N/2]=0 psihat[N/2][0]=0 psihat[0][N/2]=0 dpsix = numpy.real(numpy.fft.ifft2(psihat*kx)) dpsiy = numpy.real(numpy.fft.ifft2(psihat*ky)) u=dpsiy v=-1.0*dpsix omega=numpy.real(numpy.fft.ifft2(omegahat)) # magnetic field nlhat=0.25*numpy.fft.fft2((u+uold)*numpy.fft.ifft2((alphahat+alphahatold)*kx)+\ (v+vold)*numpy.fft.ifft2((alphahat+alphahatold)*ky)-\ (bx+bxold)*numpy.fft.ifft2((omegahat+omegahatold)*kx)-\ (by+byold)*numpy.fft.ifft2((omegahat+omegahatold)*ky)) # Implicit midpoint rule timestepping alphhat=((1/dt + 0.5*(1/Rem)*(kxx+kyy))*alphahatold \ -nlhat) \ /(1/dt -0.5*(1/Rem)*(kxx+kyy)) phihat=-alphahat/(kxx+kyy) phihat[0][0]=0 phihat[N/2][N/2]=0 phihat[N/2][0]=0 phihat[0][N/2]=0 dphix = numpy.real(numpy.fft.ifft2(phihat*kx)) dphiy = numpy.real(numpy.fft.ifft2(phihat*ky)) bx=dphiy by=-1.0*dphix alpha=numpy.real(numpy.fft.ifft2(alphahat)) # passive scalar thetax=0.5*numpy.real(numpy.fft.ifft2(kx*(thetahat+thetahatold))) thetay=0.5*numpy.real(numpy.fft.ifft2(ky*(thetahat+thetahatold))) thetahat= numpy.fft.fft2(thetaold-dt*0.5*((uold+u)*thetax+(vold+v)*thetay) ) theta=numpy.real(numpy.fft.ifft2( thetahat)) temp=abs(omega-omegacheck) chg=numpy.max(temp) temp=abs(theta-thetacheck) chg=chg+numpy.max(temp) print(chg) omegacheck=omega thetacheck=theta t+=dt src.mlab_source.scalars = theta The second program uses the implicit midpoint rule to do timestepping for the three-dimensional incompressible Magnetohydrodynamics equations and it is in listing B . B A Matlab program which finds a numerical solution to the 3D magnetohydrodynamic equations Code download % A program to solve the 3D incompressible magnetohydrodynamic equations % with periodic boundary % conditions. The program is based on the Orszag-Patterson algorithm as % documented on pg. 98 of C. Canuto, M.Y. Hussaini, A. Quarteroni and % T.A. Zhang "Spectral Methods: Evolution to Complex Geometries and % Applications to Fluid Dynamics" Springer (2007) % % The Helmholtz decomposition is used to project the magnetic field onto a % divergence free subspace. Initial work on this has been done with Damian % San Roman Alerigi % % For plotting, the program uses vol3d, which is available at: % % http://www.mathworks.com/matlabcentral/fileexchange/22940-vol3d-v2 % clear all; format compact; format short; set(0,'defaultaxesfontsize',30,'defaultaxeslinewidth',.7,... 'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,... 'defaultaxesfontweight','bold') % set up grid tic Lx = 1; % period 2*pi*L Ly = 1; % period 2*pi*L Lz = 1; % period 2*pi*L Nx = 32; % number of harmonics Ny = 32; % number of harmonics Nz = 32; % number of harmonics Nt = 50; % number of time slices dt = 2/Nt; % time step t=0; % initial time Re = 100.0; % Reynolds number Rem = 100.0; % Magnetic Reynolds number tol=10^(-10); % initialise variables x = (2*pi/Nx)*(-Nx/2:Nx/2 -1)'*Lx; % x coordinate kx = 1i*[0:Nx/2-1 0 -Nx/2+1:-1]'/Lx; % wave vector y = (2*pi/Ny)*(-Ny/2:Ny/2 -1)'*Ly; % y coordinate ky = 1i*[0:Ny/2-1 0 -Ny/2+1:-1]'/Ly; % wave vector z = (2*pi/Nz)*(-Nz/2:Nz/2 -1)'*Lz; % y coordinate kz = 1i*[0:Nz/2-1 0 -Nz/2+1:-1]'/Lz; % wave vector [xx,yy,zz]=meshgrid(x,y,z); [kxm,kym,kzm]=meshgrid(kx,ky,kz); [k2xm,k2ym,k2zm]=meshgrid(kx.^2,ky.^2,kz.^2); % initial conditions for Taylor-Green vortex % theta=0; % u=(2/sqrt(3))*sin(theta+2*pi/3)*sin(xx).*cos(yy).*cos(zz); % v=(2/sqrt(3))*sin(theta-2*pi/3)*cos(xx).*sin(yy).*cos(zz); % w=(2/sqrt(3))*sin(theta)*cos(xx).*cos(yy).*sin(zz); % initial condition sl=1; sk=1; sm=1; lamlkm=sqrt(sl.^2+sk.^2+sm.^2); u=-0.5*(lamlkm*sl*cos(sk*xx).*sin(sl*yy).*sin(sm.*zz)... +sm*sk*sin(sk*xx).*cos(sl*yy).*cos(sm.*zz))... .*exp(-t*(lamlkm^2)/Re); v=0.5*(lamlkm*sk*sin(sk*xx).*cos(sl*yy).*sin(sm.*zz)... -sm*sl*cos(sk*xx).*sin(sl*yy).*cos(sm.*zz))... *exp(-t*(lamlkm^2)/Re); w=cos(sk*xx).*cos(sl*yy).*sin(sm*zz)*exp(-t*(lamlkm^2)/Re); uhat=fftn(u); vhat=fftn(v); what=fftn(w); ux=ifftn(uhat.*kxm);uy=ifftn(uhat.*kym);uz=ifftn(uhat.*kzm); vx=ifftn(vhat.*kxm);vy=ifftn(vhat.*kym);vz=ifftn(vhat.*kzm); wx=ifftn(what.*kxm);wy=ifftn(what.*kym);wz=ifftn(what.*kzm); bx=sin(sl*yy).*sin(sm.*zz); by=sin(sm.*zz); bz=cos(sk*xx).*cos(sl*yy); bxhat=fftn(bx); byhat=fftn(by); bzhat=fftn(bz); bxx=ifftn(bxhat.*kxm);bxy=ifftn(bxhat.*kym);bxz=ifftn(bxhat.*kzm); byx=ifftn(byhat.*kxm);byy=ifftn(byhat.*kym);byz=ifftn(byhat.*kzm); bzx=ifftn(bzhat.*kxm);bzy=ifftn(bzhat.*kym);bzz=ifftn(bzhat.*kzm); % calculate fluid vorticity for plotting omegax=wy-vz; omegay=uz-wx; omegaz=vx-uy; omegatot=omegax.^2+omegay.^2+omegaz.^2; % calculate magnetic vorticity for plotting omegabx=bxy-byz; omegaby=bxz-bzx; omegabz=byx-bxy; omegabtot=omegabx.^2+omegaby.^2+omegabz.^2; figure(1); clf; n=0; subplot(2,1,1); title(['|Fluid omega|^2 ',num2str(t)]); H = vol3d('CData',omegatot,'texture','3D','XData',x,'YData',y,'ZData',z); xlabel('x'); ylabel('y'); zlabel('z'); colorbar; axis equal; axis square; view(3); subplot(2,1,2); title(['|Magnetic omega|^2 ',num2str(t)]); H = vol3d('CData',omegabtot,'texture','3D','XData',x,'YData',y,'ZData',z); xlabel('x'); ylabel('y'); zlabel('z'); colorbar; axis equal; axis square; view(3); for n=1:Nt uold=u; uxold=ux; uyold=uy; uzold=uz; vold=v; vxold=vx; vyold=vy; vzold=vz; wold=w; wxold=wx; wyold=wy; wzold=wz; bxold=bx; bxxold=bxx; bxyold=bxy; bxzold=bxz; byold=by; byxold=byx; byyold=byy; byzold=byz; bzold=bz; bzxold=bzx; bzyold=bzy; bzzold=bzz; rhsuhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*uhat; rhsvhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*vhat; rhswhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*what; rhsbxhatfix=(1/dt+(0.5/Rem)*(k2xm+k2ym+k2zm)).*bxhat; rhsbyhatfix=(1/dt+(0.5/Rem)*(k2xm+k2ym+k2zm)).*byhat; rhsbzhatfix=(1/dt+(0.5/Rem)*(k2xm+k2ym+k2zm)).*bzhat; chg=1; t=t+dt; while (chg>tol) nonlinu=0.25*((u+uold).*(ux+uxold)... +(v+vold).*(uy+uyold)... +(w+wold).*(uz+uzold)... -(bx+bxold).*(bxx+bxxold)... -(by+byold).*(bxy+bxyold)... -(bz+bzold).*(bxz+bxzold)); nonlinv=0.25*((u+uold).*(vx+vxold)... +(v+vold).*(vy+vyold)... +(w+wold).*(vz+vzold)... -(bx+bxold).*(byx+byxold)... -(by+byold).*(byy+byyold)... -(bz+bzold).*(byz+byzold)); nonlinw=0.25*((u+uold).*(wx+wxold)... +(v+vold).*(wy+wyold)... +(w+wold).*(wz+wzold)... -(bx+bxold).*(bzx+bzxold)... -(by+byold).*(bzy+bzyold)... -(bz+bzold).*(bzz+bzzold)); nonlinuhat=fftn(nonlinu); nonlinvhat=fftn(nonlinv); nonlinwhat=fftn(nonlinw); phat=-1.0*(kxm.*nonlinuhat+kym.*nonlinvhat+kzm.*nonlinwhat)... ./(k2xm+k2ym+k2zm+0.1^13); uhat=(rhsuhatfix-nonlinuhat-kxm.*phat)... ./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm)); vhat=(rhsvhatfix-nonlinvhat-kym.*phat)... ./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm)); what=(rhswhatfix-nonlinwhat-kzm.*phat)... ./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm)); ux=ifftn(uhat.*kxm);uy=ifftn(uhat.*kym);uz=ifftn(uhat.*kzm); vx=ifftn(vhat.*kxm);vy=ifftn(vhat.*kym);vz=ifftn(vhat.*kzm); wx=ifftn(what.*kxm);wy=ifftn(what.*kym);wz=ifftn(what.*kzm); utemp=u; vtemp=v; wtemp=w; u=ifftn(uhat); v=ifftn(vhat); w=ifftn(what); nonlinu=0.25*((u+uold).*(bxx+bxxold)... +(v+vold).*(bxy+bxyold)... +(w+wold).*(bxz+bxzold)... -(bx+bxold).*(ux+uxold)... -(by+byold).*(uy+uyold)... -(bz+bzold).*(uz+uzold)); nonlinv=0.25*((u+uold).*(byx+byxold)... +(v+vold).*(byy+byyold)... +(w+wold).*(byz+byzold)... -(bx+bxold).*(vx+vxold)... -(by+byold).*(vy+vyold)... -(bz+bzold).*(vz+vzold)); nonlinw=0.25*((u+uold).*(bzx+bzxold)... +(v+vold).*(bzy+bzyold)... +(w+wold).*(bzz+bzzold)... -(bx+bxold).*(wx+wxold)... -(by+byold).*(wy+wyold)... -(bz+bzold).*(wz+wzold)); nonlinuhat=fftn(nonlinu); nonlinvhat=fftn(nonlinv); nonlinwhat=fftn(nonlinw); phat=-1.0*(kxm.*nonlinuhat+kym.*nonlinvhat+kzm.*nonlinwhat)... ./(k2xm+k2ym+k2zm+0.1^13); bxhat=(rhsbxhatfix-nonlinuhat-kxm.*phat)... ./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm)); byhat=(rhsbyhatfix-nonlinvhat-kym.*phat)... ./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm)); bzhat=(rhsbzhatfix-nonlinwhat-kzm.*phat)... ./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm)); bxx=ifftn(bxhat.*kxm);bxy=ifftn(bxhat.*kym);bxz=ifftn(bxhat.*kzm); byx=ifftn(byhat.*kxm);byy=ifftn(byhat.*kym);byz=ifftn(byhat.*kzm); bzx=ifftn(bzhat.*kxm);bzy=ifftn(bzhat.*kym);bzz=ifftn(bzhat.*kzm); bxtemp=bx; bytemp=by; bztemp=bz; bx=ifftn(bxhat); by=ifftn(byhat); bz=ifftn(bzhat); chg=max(abs(utemp-u)) +max(abs(vtemp-v)) +max(abs(wtemp-w))+... max(abs(bxtemp-bx)) +max(abs(bytemp-by))+max(abs(bztemp-bz)); end % calculate vorticity for plotting omegax=wy-vz; omegay=uz-wx; omegaz=vx-uy; omegatot=omegax.^2+omegay.^2+omegaz.^2; % calculate magnetic vorticity for plotting omegabx=bxy-byz; omegaby=bxz-bzx; omegabz=byx-bxy; omegabtot=omegabx.^2+omegaby.^2+omegabz.^2; figure(1); clf; subplot(2,1,1); title(['|Fluid omega|^2 ',num2str(t)]); H = vol3d('CData',omegatot,'texture','3D','XData',x,'YData',y,'ZData',z); xlabel('x'); ylabel('y'); zlabel('z'); colorbar; axis equal; axis square; view(3); subplot(2,1,2); title(['|Magnetic omega|^2 ',num2str(t)]); H = vol3d('CData',omegabtot,'texture','3D','XData',x,'YData',y,'ZData',z); xlabel('x'); ylabel('y'); zlabel('z'); colorbar; axis equal; axis square; view(3); end B1 A Python program which finds a numerical solution to the 3D magnetohydrodynamic equations Code download #!/usr/bin/env python """ A program to solve the 3D incompressible magnetohydrodynamics equations using the implicit midpoint rule The program is based on the Orszag-Patterson algorithm as documented on pg. 98 of C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zhang "Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics" Springer (2007) The Helmholtz decomposition is used to project the magnetic field onto a divergence free subspace. Initial work on this has been done with Damian San Roman Alerigi More information on visualization can be found on the Mayavi website, in particular: http://github.enthought.com/mayavi/mayavi/mlab.html which was last checked on 6 April 2012 """ import math import numpy from mayavi import mlab import matplotlib.pyplot as plt import time # Grid Lx=1.0 # Period 2*pi*Lx Ly=1.0 # Period 2*pi*Ly Lz=1.0 # Period 2*pi*Lz Nx=64 # Number of harmonics Ny=64 # Number of harmonics Nz=64 # Number of harmonics Nt=20 # Number of time slices tmax=0.2 # Maximum time dt=tmax/Nt # time step t=0.0 # initial time Re=1.0 # Reynolds number Rem=1.0 # Magnetic Reynolds number tol=0.1**(10) # tolerance for fixed point iterations x = [i*2.0*math.pi*(Lx/Nx) for i in xrange(-Nx/2,1+Nx/2)] y = [i*2.0*math.pi*(Ly/Ny) for i in xrange(-Ny/2,1+Ny/2)] z = [i*2.0*math.pi*(Lz/Nz) for i in xrange(-Nz/2,1+Nz/2)] k_x = (1.0/Lx)*numpy.array([complex(0,1)*n for n in range(0,Nx/2) \ + [0] + range(-Nx/2+1,0)]) k_y = (1.0/Ly)*numpy.array([complex(0,1)*n for n in range(0,Ny/2) \ + [0] + range(-Ny/2+1,0)]) k_z = (1.0/Lz)*numpy.array([complex(0,1)*n for n in range(0,Nz/2) \ + [0] + range(-Nz/2+1,0)]) kxm=numpy.zeros((Nx,Ny,Nz), dtype=complex) kym=numpy.zeros((Nx,Ny,Nz), dtype=complex) kzm=numpy.zeros((Nx,Ny,Nz), dtype=complex) k2xm=numpy.zeros((Nx,Ny,Nz), dtype=float) k2ym=numpy.zeros((Nx,Ny,Nz), dtype=float) k2zm=numpy.zeros((Nx,Ny,Nz), dtype=float) xx=numpy.zeros((Nx,Ny,Nz), dtype=float) yy=numpy.zeros((Nx,Ny,Nz), dtype=float) zz=numpy.zeros((Nx,Ny,Nz), dtype=float) for i in xrange(Nx): for j in xrange(Ny): for k in xrange(Nz): kxm[i,j,k] = k_x[i] kym[i,j,k] = k_y[j] kzm[i,j,k] = k_z[k] k2xm[i,j,k] = numpy.real(k_x[i]**2) k2ym[i,j,k] = numpy.real(k_y[j]**2) k2zm[i,j,k] = numpy.real(k_z[k]**2) xx[i,j,k] = x[i] yy[i,j,k] = y[j] zz[i,j,k] = z[k] # allocate arrays u=numpy.zeros((Nx,Ny,Nz), dtype=float) uold=numpy.zeros((Nx,Ny,Nz), dtype=float) v=numpy.zeros((Nx,Ny,Nz), dtype=float) vold=numpy.zeros((Nx,Ny,Nz), dtype=float) w=numpy.zeros((Nx,Ny,Nz), dtype=float) wold=numpy.zeros((Nx,Ny,Nz), dtype=float) bx=numpy.zeros((Nx,Ny,Nz), dtype=float) bxold=numpy.zeros((Nx,Ny,Nz), dtype=float) by=numpy.zeros((Nx,Ny,Nz), dtype=float) byold=numpy.zeros((Nx,Ny,Nz), dtype=float) bz=numpy.zeros((Nx,Ny,Nz), dtype=float) bzold=numpy.zeros((Nx,Ny,Nz), dtype=float) utemp=numpy.zeros((Nx,Ny,Nz), dtype=float) vtemp=numpy.zeros((Nx,Ny,Nz), dtype=float) wtemp=numpy.zeros((Nx,Ny,Nz), dtype=float) bxtemp=numpy.zeros((Nx,Ny,Nz), dtype=float) bytemp=numpy.zeros((Nx,Ny,Nz), dtype=float) bztemp=numpy.zeros((Nx,Ny,Nz), dtype=float) omegax=numpy.zeros((Nx,Ny,Nz), dtype=float) omegay=numpy.zeros((Nx,Ny,Nz), dtype=float) omegaz=numpy.zeros((Nx,Ny,Nz), dtype=float) omegatot=numpy.zeros((Nx,Ny,Nz), dtype=float) omegabx=numpy.zeros((Nx,Ny,Nz), dtype=float) omegaby=numpy.zeros((Nx,Ny,Nz), dtype=float) omegabz=numpy.zeros((Nx,Ny,Nz), dtype=float) omegabtot=numpy.zeros((Nx,Ny,Nz), dtype=float) ux=numpy.zeros((Nx,Ny,Nz), dtype=float) uy=numpy.zeros((Nx,Ny,Nz), dtype=float) uz=numpy.zeros((Nx,Ny,Nz), dtype=float) vx=numpy.zeros((Nx,Ny,Nz), dtype=float) vy=numpy.zeros((Nx,Ny,Nz), dtype=float) vz=numpy.zeros((Nx,Ny,Nz), dtype=float) wx=numpy.zeros((Nx,Ny,Nz), dtype=float) wy=numpy.zeros((Nx,Ny,Nz), dtype=float) wz=numpy.zeros((Nx,Ny,Nz), dtype=float) uxold=numpy.zeros((Nx,Ny,Nz), dtype=float) uyold=numpy.zeros((Nx,Ny,Nz), dtype=float) uzold=numpy.zeros((Nx,Ny,Nz), dtype=float) vxold=numpy.zeros((Nx,Ny,Nz), dtype=float) vyold=numpy.zeros((Nx,Ny,Nz), dtype=float) vzold=numpy.zeros((Nx,Ny,Nz), dtype=float) wxold=numpy.zeros((Nx,Ny,Nz), dtype=float) wyold=numpy.zeros((Nx,Ny,Nz), dtype=float) wzold=numpy.zeros((Nx,Ny,Nz), dtype=float) bxx=numpy.zeros((Nx,Ny,Nz), dtype=float) bxy=numpy.zeros((Nx,Ny,Nz), dtype=float) bxz=numpy.zeros((Nx,Ny,Nz), dtype=float) byx=numpy.zeros((Nx,Ny,Nz), dtype=float) byy=numpy.zeros((Nx,Ny,Nz), dtype=float) byz=numpy.zeros((Nx,Ny,Nz), dtype=float) bzx=numpy.zeros((Nx,Ny,Nz), dtype=float) bzy=numpy.zeros((Nx,Ny,Nz), dtype=float) bzz=numpy.zeros((Nx,Ny,Nz), dtype=float) bxxold=numpy.zeros((Nx,Ny,Nz), dtype=float) bxyold=numpy.zeros((Nx,Ny,Nz), dtype=float) bxzold=numpy.zeros((Nx,Ny,Nz), dtype=float) byxold=numpy.zeros((Nx,Ny,Nz), dtype=float) byyold=numpy.zeros((Nx,Ny,Nz), dtype=float) byzold=numpy.zeros((Nx,Ny,Nz), dtype=float) bzxold=numpy.zeros((Nx,Ny,Nz), dtype=float) bzyold=numpy.zeros((Nx,Ny,Nz), dtype=float) bzzold=numpy.zeros((Nx,Ny,Nz), dtype=float) nonlinu=numpy.zeros((Nx,Ny,Nz), dtype=float) nonlinv=numpy.zeros((Nx,Ny,Nz), dtype=float) nonlinw=numpy.zeros((Nx,Ny,Nz), dtype=float) uhat=numpy.zeros((Nx,Ny,Nz), dtype=complex) what=numpy.zeros((Nx,Ny,Nz), dtype=complex) vhat=numpy.zeros((Nx,Ny,Nz), dtype=complex) bxhat=numpy.zeros((Nx,Ny,Nz), dtype=complex) byhat=numpy.zeros((Nx,Ny,Nz), dtype=complex) bzhat=numpy.zeros((Nx,Ny,Nz), dtype=complex) phat=numpy.zeros((Nx,Ny,Nz), dtype=complex) temphat=numpy.zeros((Nx,Ny,Nz), dtype=complex) rhsuhatfix=numpy.zeros((Nx,Ny,Nz), dtype=complex) rhsvhatfix=numpy.zeros((Nx,Ny,Nz), dtype=complex) rhswhatfix=numpy.zeros((Nx,Ny,Nz), dtype=complex) rhsbxhatfix=numpy.zeros((Nx,Ny,Nz), dtype=complex) rhsbyhatfix=numpy.zeros((Nx,Ny,Nz), dtype=complex) rhsbzhatfix=numpy.zeros((Nx,Ny,Nz), dtype=complex) nonlinuhat=numpy.zeros((Nx,Ny,Nz), dtype=complex) nonlinvhat=numpy.zeros((Nx,Ny,Nz), dtype=complex) nonlinwhat=numpy.zeros((Nx,Ny,Nz), dtype=complex) tdata=numpy.zeros((Nt), dtype=float) # initial conditions for Taylor-Green Vortex theta=0.0 u=(2.0/(3.0**0.5))*numpy.sin(theta+2.0*math.pi/3.0)*numpy.sin(xx)*numpy.cos(yy)*numpy.cos(zz) v=(2.0/(3.0**0.5))*numpy.sin(theta-2.0*math.pi/3.0)*numpy.cos(xx)*numpy.sin(yy)*numpy.cos(zz) w=(2.0/(3.0**0.5))*numpy.sin(theta)*numpy.cos(xx)*numpy.cos(yy)*numpy.sin(zz) # Exact solution #sl=1 #sk=1 #sm=1 #lamlkm=(sl**2.0+sk**2.0+sm**2.0)**0.5 #u=-0.5*(lamlkm*sl*numpy.cos(sk*xx)*numpy.sin(sl*yy)*numpy.sin(sm*zz) \ #+sm*sk*numpy.sin(sk*xx)*numpy.cos(sl*yy)*numpy.cos(sm*zz))*numpy.exp(-t*(lamlkm**2.0)/Rey) #v= 0.5*(lamlkm*sk*numpy.sin(sk*xx)*numpy.cos(sl*yy)*numpy.sin(sm*zz) \ #-sm*sl*numpy.cos(sk*xx)*numpy.sin(sl*yy)*numpy.cos(sm*zz))*numpy.exp(-t*(lamlkm**2.0)/Rey) #w= numpy.cos(sk*xx)*numpy.cos(sl*yy)*numpy.sin(sm*zz)*numpy.exp(-t*(lamlkm**2.0)/Rey) # initial fluid field terms uhat=numpy.fft.fftn(u) vhat=numpy.fft.fftn(v) what=numpy.fft.fftn(w) temphat=kxm*uhat ux=numpy.real(numpy.fft.ifftn(temphat)) temphat=kym*uhat uy=numpy.real(numpy.fft.ifftn(temphat)) temphat=kzm*uhat uz=numpy.real(numpy.fft.ifftn(temphat)) temphat=kxm*vhat vx=numpy.real(numpy.fft.ifftn(temphat)) temphat=kym*vhat vy=numpy.real(numpy.fft.ifftn(temphat)) temphat=kzm*vhat vz=numpy.real(numpy.fft.ifftn(temphat)) temphat=kxm*what wx=numpy.real(numpy.fft.ifftn(temphat)) temphat=kym*what wy=numpy.real(numpy.fft.ifftn(temphat)) temphat=kzm*what wz=numpy.real(numpy.fft.ifftn(temphat)) # Calculate fluid vorticity for plotting omegax=wy-vz omegay=uz-wx omegaz=vx-uy omegatot=omegax**2.0 + omegay**2.0 + omegaz**2.0 # initial magnetic field terms bxhat=numpy.fft.fftn(bx) byhat=numpy.fft.fftn(by) bzhat=numpy.fft.fftn(bz) temphat=kxm*bxhat bxx=numpy.real(numpy.fft.ifftn(temphat)) temphat=kym*bxhat bxy=numpy.real(numpy.fft.ifftn(temphat)) temphat=kzm*bxhat bxz=numpy.real(numpy.fft.ifftn(temphat)) temphat=kxm*byhat byx=numpy.real(numpy.fft.ifftn(temphat)) temphat=kym*byhat byy=numpy.real(numpy.fft.ifftn(temphat)) temphat=kzm*byhat byz=numpy.real(numpy.fft.ifftn(temphat)) temphat=kxm*bzhat bzx=numpy.real(numpy.fft.ifftn(temphat)) temphat=kym*bzhat bzy=numpy.real(numpy.fft.ifftn(temphat)) temphat=kzm*bzhat bzz=numpy.real(numpy.fft.ifftn(temphat)) # Calculate magnetic vorticity for plotting omegabx=bzy-byz omegaby=bxz-bzx omegabz=byx-bxy omegabtot=omegabx**2.0 + omegaby**2.0 + omegabz**2.0 #src=mlab.contour3d(xx,yy,zz,u,colormap='jet',opacity=0.1,contours=4) src = mlab.pipeline.scalar_field(xx,yy,zz,omegatot,colormap='YlGnBu') mlab.pipeline.iso_surface(src, contours=[omegatot.min()+0.1*omegatot.ptp(), ], \ colormap='YlGnBu',opacity=0.85) mlab.pipeline.iso_surface(src, contours=[omegatot.max()-0.1*omegatot.ptp(), ], \ colormap='YlGnBu',opacity=1.0) #src = mlab.pipeline.scalar_field(xx,yy,zz,omegabtot,colormap='YlGnBu') #mlab.pipeline.iso_surface(src, contours=[omegabtot.min()+0.1*omegatot.ptp(), ], \ # colormap='YlGnBu',opacity=0.85) #mlab.pipeline.iso_surface(src, contours=[omegabtot.max()-0.1*omegatot.ptp(), ], \ # colormap='YlGnBu',opacity=1.0) mlab.pipeline.image_plane_widget(src,plane_orientation='z_axes', \ slice_index=Nz/2,colormap='YlGnBu', \ opacity=0.01) mlab.pipeline.image_plane_widget(src,plane_orientation='y_axes', \ slice_index=Ny/2,colormap='YlGnBu', \ opacity=0.01) mlab.pipeline.image_plane_widget(src,plane_orientation='x_axes', \ slice_index=Nx/2,colormap='YlGnBu', \ opacity=0.01) mlab.scalarbar() mlab.xlabel('x',object=src) mlab.ylabel('y',object=src) mlab.zlabel('z',object=src) t=0.0 tdata[0]=t #solve pde and plot results for n in xrange(Nt): uold=u uxold=ux uyold=uy uzold=uz vold=v vxold=vx vyold=vy vzold=vz wold=w wxold=wx wyold=wy wzold=wz bxold=bx bxxold=bxx bxyold=bxy bxzold=bxz byold=by byxold=byx byyold=byy byzold=byz bzold=bz bzxold=bzx bzyold=bzy bzzold=bzz rhsuhatfix=(1.0/dt + (0.5/Re)*(k2xm+k2ym+k2zm))*uhat rhsvhatfix=(1.0/dt + (0.5/Re)*(k2xm+k2ym+k2zm))*vhat rhswhatfix=(1.0/dt + (0.5/Re)*(k2xm+k2ym+k2zm))*what rhsbxhatfix=(1.0/dt + (0.5/Rem)*(k2xm+k2ym+k2zm))*bxhat rhsbyhatfix=(1.0/dt + (0.5/Rem)*(k2xm+k2ym+k2zm))*byhat rhsbzhatfix=(1.0/dt + (0.5/Rem)*(k2xm+k2ym+k2zm))*bzhat chg=1.0 t=t+dt while(chg>tol): # Fluid field nonlinu=0.25*((u+uold)*(ux+uxold)+(v+vold)*(uy+uyold)+(w+wold)*(uz+uzold)-\ (bx+bxold)*(bxx+bxxold)-(by+byold)*(bxy+bxyold)-(bz+bzold)*(bxz+bxzold)) nonlinv=0.25*((u+uold)*(vx+vxold)+(v+vold)*(vy+vyold)+(w+wold)*(vz+vzold)-\ (bx+bxold)*(byx+byxold)-(by+byold)*(byy+byyold)-(bz+bzold)*(byz+byzold)) nonlinw=0.25*((u+uold)*(wx+wxold)+(v+vold)*(wy+wyold)+(w+wold)*(wz+wzold)-\ (bx+bxold)*(bzx+bzxold)-(by+byold)*(bzy+bzyold)-(bz+bzold)*(bzz+bzzold)) nonlinuhat=numpy.fft.fftn(nonlinu) nonlinvhat=numpy.fft.fftn(nonlinv) nonlinwhat=numpy.fft.fftn(nonlinw) phat=-1.0*(kxm*nonlinuhat+kym*nonlinvhat+kzm*nonlinwhat)/(k2xm+k2ym+k2zm+0.1**13) uhat=(rhsuhatfix-nonlinuhat-kxm*phat)/(1.0/dt - (0.5/Re)*(k2xm+k2ym+k2zm)) vhat=(rhsvhatfix-nonlinvhat-kym*phat)/(1.0/dt - (0.5/Re)*(k2xm+k2ym+k2zm)) what=(rhswhatfix-nonlinwhat-kzm*phat)/(1.0/dt - (0.5/Re)*(k2xm+k2ym+k2zm)) temphat=kxm*uhat ux=numpy.real(numpy.fft.ifftn(temphat)) temphat=kym*uhat uy=numpy.real(numpy.fft.ifftn(temphat)) temphat=kzm*uhat uz=numpy.real(numpy.fft.ifftn(temphat)) temphat=kxm*vhat vx=numpy.real(numpy.fft.ifftn(temphat)) temphat=kym*vhat vy=numpy.real(numpy.fft.ifftn(temphat)) temphat=kzm*vhat vz=numpy.real(numpy.fft.ifftn(temphat)) temphat=kxm*what wx=numpy.real(numpy.fft.ifftn(temphat)) temphat=kym*what wy=numpy.real(numpy.fft.ifftn(temphat)) temphat=kzm*what wz=numpy.real(numpy.fft.ifftn(temphat)) utemp=u vtemp=v wtemp=w u=numpy.real(numpy.fft.ifftn(uhat)) v=numpy.real(numpy.fft.ifftn(vhat)) w=numpy.real(numpy.fft.ifftn(what)) # Magnetic field nonlinu=0.25*((u+uold)*(bxx+bxxold)+(v+vold)*(bxy+bxyold)+(w+wold)*(bxz+bxzold)-\ (bx+bxold)*(ux+uxold)-(by+byold)*(uy+uyold)-(bz+bzold)*(uz+uzold)) nonlinv=0.25*((u+uold)*(byx+byxold)+(v+vold)*(byy+byyold)+(w+wold)*(byz+byzold)-\ (bx+bxold)*(vx+vxold)-(by+byold)*(vy+vyold)-(bz+bzold)*(vz+vzold)) nonlinw=0.25*((u+uold)*(bzx+bzxold)+(v+vold)*(bzy+bzyold)+(w+wold)*(bzz+bzzold)-\ (bx+bxold)*(wx+wxold)-(by+byold)*(wy+wyold)-(bz+bzold)*(wz+wzold)) nonlinuhat=numpy.fft.fftn(nonlinu) nonlinvhat=numpy.fft.fftn(nonlinv) nonlinwhat=numpy.fft.fftn(nonlinw) phat=-1.0*(kxm*nonlinuhat+kym*nonlinvhat+kzm*nonlinwhat)/(k2xm+k2ym+k2zm+0.1**13) bxhat=(rhsbxhatfix-nonlinuhat-kxm*phat)/(1.0/dt - (0.5/Rem)*(k2xm+k2ym+k2zm)) byhat=(rhsbyhatfix-nonlinvhat-kym*phat)/(1.0/dt - (0.5/Rem)*(k2xm+k2ym+k2zm)) bzhat=(rhsbzhatfix-nonlinwhat-kzm*phat)/(1.0/dt - (0.5/Rem)*(k2xm+k2ym+k2zm)) temphat=kxm*bxhat bxx=numpy.real(numpy.fft.ifftn(temphat)) temphat=kym*bxhat bxy=numpy.real(numpy.fft.ifftn(temphat)) temphat=kzm*bxhat bxz=numpy.real(numpy.fft.ifftn(temphat)) temphat=kxm*byhat byx=numpy.real(numpy.fft.ifftn(temphat)) temphat=kym*byhat byy=numpy.real(numpy.fft.ifftn(temphat)) temphat=kzm*byhat byz=numpy.real(numpy.fft.ifftn(temphat)) temphat=kxm*bzhat bzx=numpy.real(numpy.fft.ifftn(temphat)) temphat=kym*bzhat bzy=numpy.real(numpy.fft.ifftn(temphat)) temphat=kzm*bzhat bzz=numpy.real(numpy.fft.ifftn(temphat)) bxtemp=bx bytemp=by bztemp=bz bx=numpy.real(numpy.fft.ifftn(bxhat)) by=numpy.real(numpy.fft.ifftn(byhat)) bz=numpy.real(numpy.fft.ifftn(bzhat)) chg=numpy.max(abs(u-utemp))+numpy.max(abs(v-vtemp))+numpy.max(abs(w-wtemp))+\ numpy.max(abs(bx-bxtemp))+numpy.max(abs(by-bytemp))+numpy.max(abs(bz-bztemp)) # calculate vorticity for plotting omegax=wy-vz omegay=uz-wx omegaz=vx-uy omegatot=omegax**2.0 + omegay**2.0 + omegaz**2.0 src.mlab_source.scalars = omegatot tdata[n]=t omegabx=bzy-byz omegaby=bxz-bzx omegabz=byx-bxy omegabtot=omegabx**2.0 + omegaby**2.0 + omegabz**2.0 #src.mlab_source.scalars = omegatot Fortran Programs The programs in this section are based on those for the Navier-Stokes equations in https://en.wikibooks.org/wiki/Parallel_Spectral_Numerical_Methods/The_Two-_and_Three-Dimensional_Navier-Stokes_Equations. C A serial Fortran program which finds a numerical solution to the 2D magnetohydrodynamic equations Code download PROGRAM main !-------------------------------------------------------------------- ! ! ! PURPOSE ! ! This program numerically solves the 2D incompressible magnetohydrodynamics ! equations on a Square Domain [0,1]x[0,1] using pseudo-spectral methods and ! implicit midpoint rule timestepping. A passive scalar is also advected ! ! Periodic free-slip boundary conditions and Initial conditions: ! u(x,y,0)=sin(2*pi*x)cos(2*pi*y) ! v(x,y,0)=-cos(2*pi*x)sin(2*pi*y) ! ! .. Parameters .. ! Nx = number of modes in x - power of 2 for FFT ! Ny = number of modes in y - power of 2 for FFT ! Nt = number of timesteps to take ! Tmax = maximum simulation time ! FFTW_IN_PLACE = value for FFTW input ! FFTW_MEASURE = value for FFTW input ! FFTW_EXHAUSTIVE = value for FFTW input ! FFTW_PATIENT = value for FFTW input ! FFTW_ESTIMATE = value for FFTW input ! FFTW_FORWARD = value for FFTW input ! FFTW_BACKWARD = value for FFTW input ! pi = 3.14159265358979323846264338327950288419716939937510d0 ! mu = viscosity ! rho = density ! .. Scalars .. ! i = loop counter in x direction ! j = loop counter in y direction ! n = loop counter for timesteps direction ! allocatestatus = error indicator during allocation ! count = keep track of information written to disk ! iol = size of array to write to disk ! start = variable to record start time of program ! finish = variable to record end time of program ! count_rate = variable for clock count rate ! planfx = Forward 1d fft plan in x ! planbx = Backward 1d fft plan in x ! planfy = Forward 1d fft plan in y ! planby = Backward 1d fft plan in y ! dt = timestep ! .. Arrays .. ! u = velocity in x direction ! uold = velocity in x direction at previous timestep ! v = velocity in y direction ! vold = velocity in y direction at previous timestep ! u_y = y derivative of velocity in x direction ! v_x = x derivative of velocity in y direction ! bx = x component of magnetic field ! by = y component of magnetic field ! bx = x component of magnetic field at previous time step ! by = y component of magnetic field at previous time step ! omeg = vorticity in real space ! omegold = vorticity in real space at previous ! iterate ! omegcheck = store of vorticity at previous iterate ! omegoldhat = 2D Fourier transform of vorticity at previous ! iterate ! omegoldhat_x = x-derivative of vorticity in Fourier space ! at previous iterate ! omegold_x = x-derivative of vorticity in real space ! at previous iterate ! omegoldhat_y = y-derivative of vorticity in Fourier space ! at previous iterate ! omegold_y = y-derivative of vorticity in real space ! at previous iterate ! nlold = nonlinear term in real space ! at previous iterate ! nloldhat = nonlinear term in Fourier space ! at previous iterate ! omeghat = 2D Fourier transform of vorticity ! at next iterate ! omeghat_x = x-derivative of vorticity in Fourier space ! at next timestep ! omeghat_y = y-derivative of vorticity in Fourier space ! at next timestep ! omeg_x = x-derivative of vorticity in real space ! at next timestep ! omeg_y = y-derivative of vorticity in real space ! at next timestep ! alpha = magnetic vorticity in real space ! alphaold = magnetic vorticity in real space at previous ! iterate ! alphacheck = store of vorticity potential at previous iterate ! alphaoldhat = 2D Fourier transform of vorticity potential at previous ! iterate ! alphaoldhat_x = x-derivative of magnetic vorticity in Fourier space ! at previous iterate ! alphaold_x = x-derivative of magnetic vorticity in real space ! at previous iterate ! alphaoldhat_y = y-derivative of magnetic vorticity in Fourier space ! at previous iterate ! alphaold_y = y-derivative of magnetic vorticity in real space ! alphahat = 2D Fourier transform of magnetic vorticity ! at next iterate ! alphahat_x = x-derivative of magnetic vorticity in Fourier space ! at next timestep ! alphahat_y = y-derivative of magnetic vorticity in Fourier space ! at next timestep ! alpha_x = x-derivative of magnetic vorticity in real space ! at next timestep ! alpha_y = y-derivative of magnetic vorticity in real space ! at next timestep ! psihat = streamfunction at next timestep ! psi_x = x-derivative of streamfunction in real space ! at next timestep ! psi_y = y-derivative of streamfunction in real space ! at next timestep ! phi = Magnetic potential at next timestep ! phihat = 2D Fourier transform of magnetic potential at next timestep ! phi_x = x-derivative of streamfunction in real space ! at next timestep ! phi_y = y-derivative of streamfunction in real space ! at next timestep ! .. Vectors .. ! kx = fourier frequencies in x direction ! ky = fourier frequencies in y direction ! kxx = square of fourier frequencies in x direction ! kyy = square of fourier frequencies in y direction ! x = x locations ! y = y locations ! time = times at which save data ! name_config = array to store filename for data to be saved ! fftfx = array to setup x Fourier transform ! fftbx = array to setup y Fourier transform ! REFERENCES ! ! ACKNOWLEDGEMENTS ! ! ACCURACY ! ! ERROR INDICATORS AND WARNINGS ! ! FURTHER COMMENTS ! This program has not been optimized to use the least amount of memory ! but is intended as an example only for which all states can be saved !-------------------------------------------------------------------- ! External routines required ! ! External libraries required ! FFTW3 -- Fast Fourier Transform in the West Library ! (http://www.fftw.org/) ! declare variables IMPLICIT NONE INTEGER(kind=4), PARAMETER :: Nx=64 INTEGER(kind=4), PARAMETER :: Ny=64 REAL(kind=8), PARAMETER :: dt=0.00125 REAL(kind=8), PARAMETER & :: pi=3.14159265358979323846264338327950288419716939937510 REAL(kind=8), PARAMETER :: Re=1.0d0 REAL(kind=8), PARAMETER :: Rem=1.0d0 REAL(kind=8), PARAMETER :: tol=0.1d0**10 REAL(kind=8) :: chg INTEGER(kind=4), PARAMETER :: nplots=50 REAL(kind=8), DIMENSION(:), ALLOCATABLE :: time COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE :: kx,kxx COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE :: ky,kyy REAL(kind=8), DIMENSION(:), ALLOCATABLE :: x REAL(kind=8), DIMENSION(:), ALLOCATABLE :: y COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE :: & u,uold,v,vold,u_y,v_x,omegold, omegcheck, omeg,& omegoldhat, omegold_x, omegold_y, & omeghat, omeg_x, omeg_y,& nl, nlhat, psihat, psi_x, psi_y,& theta, thetaold, thetax, thetay, thetacheck, thetahat, thetaoldhat, & alpha, alphaold, alpha_x, alpha_y, bx, by, bxold, byold, alphahat, alphaoldhat, alphacheck, & phi, phihat, phi_x, phi_y, temp INTEGER(kind=4) :: i,j,k,n, allocatestatus, count, iol INTEGER(kind=4) :: start, finish, count_rate INTEGER(kind=4), PARAMETER :: FFTW_IN_PLACE = 8, FFTW_MEASURE = 0, & FFTW_EXHAUSTIVE = 8, FFTW_PATIENT = 32, & FFTW_ESTIMATE = 64 INTEGER(kind=4),PARAMETER :: FFTW_FORWARD = -1, FFTW_BACKWARD=1 COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE :: fftfx,fftbx INTEGER(kind=8) :: planfxy,planbxy CHARACTER*100 :: name_config CALL system_clock(start,count_rate) ALLOCATE(time(1:nplots),kx(1:Nx),kxx(1:Nx),ky(1:Ny),kyy(1:Ny),x(1:Nx),y(1:Ny),& u(1:Nx,1:Ny),uold(1:Nx,1:Ny),v(1:Nx,1:Ny),vold(1:Nx,1:Ny),u_y(1:Nx,1:Ny),& v_x(1:Nx,1:Ny),omegold(1:Nx,1:Ny),omegcheck(1:Nx,1:Ny), omeg(1:Nx,1:Ny),& omegoldhat(1:Nx,1:Ny), omegold_x(1:Nx,1:Ny), omegold_y(1:Nx,1:Ny), & omeghat(1:Nx,1:Ny), omeg_x(1:Nx,1:Ny), omeg_y(1:Nx,1:Ny),& nl(1:Nx,1:Ny), nlhat(1:Nx,1:Ny), psihat(1:Nx,1:Ny), & psi_x(1:Nx,1:Ny), psi_y(1:Nx,1:Ny),& theta(1:Nx,1:Ny), thetaold(1:Nx,1:Ny), thetax(1:Nx,1:Ny), thetay(1:Nx,1:Ny), & thetacheck(1:Nx,1:Ny), thetahat(1:Nx,1:Ny), thetaoldhat(1:Nx,1:Ny), & alpha(1:Nx,1:Ny), alphaold(1:Nx,1:Ny), alpha_x(1:Nx,1:Ny), alpha_y(1:Nx,1:Ny),& bx(1:Nx,1:Ny), by(1:Nx,1:Ny), bxold(1:Nx,1:Ny), byold(1:Nx,1:Ny),& alphahat(1:Nx,1:Ny), alphaoldhat(1:Nx,1:Ny),alphacheck(1:Nx,1:Ny), & phi(1:Nx,1:Ny), phihat(1:Nx,1:Ny), & phi_x(1:Nx,1:Ny), phi_y(1:Nx,1:Ny), temp(1:Nx,1:Ny),& fftfx(1:Nx,1:Ny),fftbx(1:Nx,1:Ny),stat=AllocateStatus) IF (AllocateStatus .ne. 0) STOP PRINT *,'allocated space' ! set up ffts CALL dfftw_plan_dft_2d_(planfxy,Nx,Ny,fftfx(1:Nx,1:Ny),fftbx(1:Nx,1:Ny),& FFTW_FORWARD,FFTW_EXHAUSTIVE) CALL dfftw_plan_dft_2d_(planbxy,Nx,Ny,fftbx(1:Nx,1:Ny),fftfx(1:Nx,1:Ny),& FFTW_BACKWARD,FFTW_EXHAUSTIVE) ! setup fourier frequencies in x-direction DO i=1,1+Nx/2 kx(i)= 2.0d0*pi*cmplx(0.0d0,1.0d0)*REAL(i-1,kind(0d0)) END DO kx(1+Nx/2)=0.0d0 DO i = 1,Nx/2 -1 kx(i+1+Nx/2)=-kx(1-i+Nx/2) END DO DO i=1,Nx kxx(i)=kx(i)*kx(i) END DO DO i=1,Nx x(i)=REAL(i-1,kind(0d0))/REAL(Nx,kind(0d0)) END DO ! setup fourier frequencies in y-direction DO j=1,1+Ny/2 ky(j)= 2.0d0*pi*cmplx(0.0d0,1.0d0)*REAL(j-1,kind(0d0)) END DO ky(1+Ny/2)=0.0d0 DO j = 1,Ny/2 -1 ky(j+1+Ny/2)=-ky(1-j+Ny/2) END DO DO j=1,Ny kyy(j)=ky(j)*ky(j) END DO DO j=1,Ny y(j)=REAL(j-1,kind(0d0))/REAL(Ny,kind(0d0)) END DO PRINT *,'Setup grid and fourier frequencies' DO j=1,Ny DO i=1,Nx u(i,j)=sin(2.0d0*pi*x(i))*cos(2.0d0*pi*y(j)) v(i,j)=-cos(2.0d0*pi*x(i))*sin(2.0d0*pi*y(j)) u_y(i,j)=-2.0d0*pi*sin(2.0d0*pi*x(i))*sin(2.0d0*pi*y(j)) v_x(i,j)=2.0d0*pi*sin(2.0d0*pi*x(i))*sin(2.0d0*pi*y(j)) omeg(i,j)=v_x(i,j)-u_y(i,j) theta(i,j)=exp(-2.0*((4.0*pi*(x(i)-0.3))**2+(4.0*pi*(y(j)-0.3))**2)) alpha(i,j)=sin(4.0*pi*x(i))*cos(6.0*pi*y(j)) END DO END DO ! Vorticity to Fourier Space CALL dfftw_execute_dft_(planfxy,omeg(1:Nx,1:Ny),omeghat(1:Nx,1:Ny)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! time(1)=0.0d0 PRINT *,'Got initial data, starting timestepping' DO n=1,nplots chg=1 ! save old values uold(1:Nx,1:Ny)=u(1:Nx,1:Ny) vold(1:Nx,1:Ny)=v(1:Nx,1:Ny) omegold(1:Nx,1:Ny)=omeg(1:Nx,1:Ny) omegcheck(1:Nx,1:Ny)=omeg(1:Nx,1:Ny) omegoldhat(1:Nx,1:Ny)=omeghat(1:Nx,1:Ny) thetaold(1:Nx,1:Ny)=theta(1:Nx,1:Ny) thetaoldhat(1:Nx,1:Ny)=thetahat(1:Nx,1:Ny) thetacheck(1:Nx,1:Ny)=theta(1:Nx,1:Ny) bxold(1:Nx,1:Ny)=bx(1:Nx,1:Ny) byold(1:Nx,1:Ny)=by(1:Nx,1:Ny) alphaoldhat(1:Nx,1:Ny)=alphahat(1:Nx,1:Ny) alphaold(1:Nx,1:Ny)=alpha(1:Nx,1:Ny) alphacheck(1:Nx,1:Ny)=alpha(1:Nx,1:Ny) DO WHILE (chg>tol) ! Fluid flow field !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!nonlinear fixed (n,k+1)!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! obtain \hat{\omega}_x^{n+1,k} DO i=1,Nx; DO j=1,Ny; temp(i,j)=0.5d0*(omeghat(i,j)+omegoldhat(i,j))*kx(i)/REAL(Nx*Ny,kind(0d0)) END DO; END DO ! convert back to real space CALL dfftw_execute_dft_(planbxy,temp(1:Nx,1:Ny),omeg_x(1:Nx,1:Ny)) ! obtain \hat{\omega}_y^{n+1,k} DO i=1,Nx; DO j=1,Ny; temp(i,j)=0.5d0*(omeghat(i,j)+omegoldhat(i,j))*ky(j)/REAL(Nx*Ny,kind(0d0)) END DO; END DO ! convert back to real space CALL dfftw_execute_dft_(planbxy,temp(1:Nx,1:Ny),omeg_y(1:Nx,1:Ny)) ! calculate nonlinear term in real space DO i=1,Nx; DO j=1,Ny; nl(i,j)=0.5d0*( (u(i,j) + uold(i,j))*omeg_x(i,j)+& (v(i,j) + vold(i,j))*omeg_y(i,j)-& (bx(i,j) + bxold(i,j))*alpha_x(i,j)-& (by(i,j) + byold(i,j))*alpha_y(i,j) ) END DO; END DO ! convert back to fourier CALL dfftw_execute_dft_(planfxy,nl(1:Nx,1:Ny),nlhat(1:Nx,1:Ny)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! obtain \hat{\omega}^{n+1,k+1} with implicit midpoint rule timestepping DO i=1,Nx; DO j=1,Ny; omeghat(i,j)=( (1.0d0/dt+0.5d0*(1.0d0/Re)*(kxx(i)+kyy(j)))& *omegoldhat(i,j) - nlhat(i,j))/& (1.0d0/dt-0.5d0*(1.0d0/Re)*(kxx(i)+kyy(j))) END DO; END DO ! calculate \hat{\psi}^{n+1,k+1} DO i=1,Nx; DO j=1,Ny; psihat(i,j)=-omeghat(i,j)/(kxx(i)+kyy(j)) END DO; END DO psihat(1,1)=0.0d0 psihat(Nx/2+1,Ny/2+1)=0.0d0 psihat(Nx/2+1,1)=0.0d0 psihat(1,Ny/2+1)=0.0d0 ! obtain \psi_x^{n+1,k+1} and \psi_y^{n+1,k+1} DO i=1,Nx; DO j=1,Ny; temp(i,j)=psihat(i,j)*kx(i)/REAL(Nx*Ny,kind(0d0)) END DO; END DO CALL dfftw_execute_dft_(planbxy,temp(1:Nx,1:Ny),psi_x(1:Nx,1:Ny)) DO i=1,Nx; DO j=1,Ny temp(i,j)=psihat(i,j)*ky(j)/REAL(Nx*Ny,kind(0d0)) END DO; END DO CALL dfftw_execute_dft_(planbxy,temp(1:Ny,1:Ny),psi_y(1:Ny,1:Ny)) ! obtain \omega^{n+1,k+1} CALL dfftw_execute_dft_(planbxy,omeghat(1:Nx,1:Ny),omeg(1:Nx,1:Ny)) DO i=1,Nx; DO j=1,Ny omeg(i,j)=omeg(i,j)/REAL(Nx*Ny,kind(0d0)) END DO; END DO ! obtain u^{n+1,k+1} and v^{n+1,k+1} using stream function (\psi) in real space DO i=1,Nx; DO j=1,Ny u(i,j)=psi_y(i,j) v(i,j)=-psi_x(i,j) END DO; END DO ! Magnetic field !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!nonlinear fixed (n,k+1)!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! obtain \hat{\omega}_x^{n+1,k} DO i=1,Nx; DO j=1,Ny; temp(i,j)=0.5*(alphahat(i,j)+alphaoldhat(i,j))*kx(i)/REAL(Nx*Ny,kind(0d0)) END DO; END DO ! convert back to real space CALL dfftw_execute_dft_(planbxy,temp(1:Nx,1:Ny),alpha_x(1:Nx,1:Ny)) ! obtain \hat{\omega}_y^{n+1,k} DO i=1,Nx; DO j=1,Ny; temp(i,j)=0.5*(alphahat(i,j)+alphaoldhat(i,j))*ky(j)/REAL(Nx*Ny,kind(0d0)) END DO; END DO ! convert back to real space CALL dfftw_execute_dft_(planbxy,temp(1:Nx,1:Ny),omeg_y(1:Nx,1:Ny)) ! calculate nonlinear term in real space DO i=1,Nx; DO j=1,Ny nl(i,j)=0.5d0*( (u(i,j) + uold(i,j))*alpha_x(i,j)+& (v(i,j) + vold(i,j))*alpha_y(i,j)-& (bx(i,j) + bxold(i,j))*omeg_x(i,j)-& (by(i,j) + byold(i,j))*omeg_y(i,j) ) END DO; END DO ! convert back to fourier CALL dfftw_execute_dft_(planfxy,nl(1:Nx,1:Ny),nlhat(1:Nx,1:Ny)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! obtain \hat{\omega}^{n+1,k+1} with implicit midpoint rule timestepping DO i=1,Nx; DO j=1,Ny alphahat(i,j)=( (1.0d0/dt+0.5d0*(1.0d0/Rem)*(kxx(i)+kyy(j)))& *alphaoldhat(i,j) - nlhat(i,j))/& (1.0d0/dt-0.5d0*(1.0d0/Rem)*(kxx(i)+kyy(j))) END DO; END DO ! calculate \hat{\psi}^{n+1,k+1} DO i=1,Nx; DO j=1,Ny phihat(i,j)=-alphahat(i,j)/(kxx(i)+kyy(j)) END DO; END DO phihat(1,1)=0.0d0 phihat(Nx/2+1,Ny/2+1)=0.0d0 phihat(Nx/2+1,1)=0.0d0 phihat(1,Ny/2+1)=0.0d0 ! obtain \psi_x^{n+1,k+1} and \psi_y^{n+1,k+1} DO i=1,Nx; DO j=1,Ny; temp(i,j)=phihat(i,j)*kx(i)/REAL(Nx*Ny,kind(0d0)) END DO; END DO CALL dfftw_execute_dft_(planbxy,temp(1:Nx,1:Ny),phi_x(1:Nx,1:Ny)) DO i=1,Nx; DO j=1,Ny temp(i,j)=phihat(i,j)*ky(j)/REAL(Nx*Ny,kind(0d0)) END DO; END DO CALL dfftw_execute_dft_(planbxy,temp(1:Ny,1:Ny),phi_y(1:Ny,1:Ny)) ! obtain \omega^{n+1,k+1} CALL dfftw_execute_dft_(planbxy,alphahat(1:Nx,1:Ny),alpha(1:Nx,1:Ny)) DO i=1,Nx; DO j=1,Ny alpha(i,j)=alpha(i,j)/REAL(Nx*Ny,kind(0d0)) END DO; END DO ! obtain u^{n+1,k+1} and v^{n+1,k+1} using stream function (\psi) in real space DO i=1,Nx; DO j=1,Ny bx(i,j)=phi_y(i,j) END DO; END DO DO i=1,Nx; DO j=1,Ny by(i,j)=-phi_x(i,j) END DO; END DO ! check for convergence chg=maxval(abs(omeg-omegcheck)) + maxval(abs(alpha-alphacheck)) ! saves {n+1,k+1} to {n,k} for next iteration DO i=1,Nx; DO j=1,Ny omegcheck(i,j)=omeg(i,j) END DO; END DO DO i=1,Nx; DO j=1,Ny alphacheck(i,j)=alpha(i,j) END DO; END DO END DO time(n+1)=time(n)+dt PRINT *,'TIME ',time(n+1) END DO name_config = 'omegafinal.datbin' INQUIRE(iolength=iol) omeg(1,1) OPEN(unit=11,FILE=name_config,form="unformatted", access="direct",recl=iol) count = 1 DO j=1,Ny DO i=1,Nx WRITE(11,rec=count) REAL(omeg(i,j),KIND(0d0)) count=count+1 END DO END DO CLOSE(11) name_config = 'alphfinal.datbin' INQUIRE(iolength=iol) alpha(1,1) OPEN(unit=11,FILE=name_config,form="unformatted", access="direct",recl=iol) count = 1 DO j=1,Ny DO i=1,Nx WRITE(11,rec=count) REAL(alpha(i,j),KIND(0d0)) count=count+1 END DO END DO CLOSE(11) name_config = 'xcoord.dat' OPEN(unit=11,FILE=name_config,status="UNKNOWN") REWIND(11) DO i=1,Nx WRITE(11,*) x(i) END DO CLOSE(11) name_config = 'ycoord.dat' OPEN(unit=11,FILE=name_config,status="UNKNOWN") REWIND(11) DO j=1,Ny WRITE(11,*) y(j) END DO CLOSE(11) CALL dfftw_destroy_plan_(planfxy) CALL dfftw_destroy_plan_(planbxy) CALL dfftw_cleanup_() DEALLOCATE(time,kx,kxx,ky,kyy,x,y,& u,uold,v,vold,u_y,& v_x,omegold,omegcheck, omeg,& omegoldhat, omegold_x, & omegold_y, & omeghat, omeg_x,& omeg_y, nl, nlhat, psihat, & psi_x, psi_y,& theta, thetaold, thetax, thetay, & thetacheck, thetahat, thetaoldhat, & alpha, alphaold, alpha_x, alpha_y,& bx, by,& bxold, byold, alphahat, alphaoldhat, & alphacheck, phi, phihat, & phi_x, phi_y, temp,& fftfx,fftbx,stat=AllocateStatus) IF (AllocateStatus .ne. 0) STOP PRINT *,'Program execution complete' END PROGRAM main D A makefile to compile the program which finds a numerical solution to the 2D magnetohydrodynamic equations Code download COMPILER = gfortran FLAGS = -O0 LIBS = -L${FFTW_LINK} -lfftw3 -lm SOURCES = magnetohydrodynamics2D.f90 mhd2d: $(SOURCES) ${COMPILER} -o mhd2d $(FLAGS) $(SOURCES) $(LIBS) clean: rm -f *.o rm -f *.mod clobber: rm -f mhd2d Parallel Fortran Programs To make a parallel program, the library 2DECOMP&FFT is used for the domain decomposition and the Fourier transforms. 2DECOMP&FFT is available from http://www.2decomp.org/index.html. The website includes some examples which indicate how this library should be used, in particular the sample code at http://www.2decomp.org/case_study1.html is a very helpful indication of how one converts a code that uses FFTW to one that uses MPI and the aforementioned library. To get the programs to work, first compile 2DECOMP&FFT, update the link in the makefile, and then compile the parallel Fortran program. E A parallel Fortran program which finds a numerical solution to the 3D magnetohydrodynamic equations Code download PROGRAM main !----------------------------------------------------------------------------------- ! ! ! PURPOSE ! ! This program numerically solves the 3D incompressible magnetohydrodynamic equations ! on a Cubic Domain [0,2pi]x[0,2pi]x[0,2pi] using fourier pseudo-spectral methods and ! Implicit Midpoint rule timestepping. ! ! Analytical Solution: ! u(x,y,z,t)=-0.25*(cos(x)sin(y)sin(z)+sin(x)cos(y)cos(z))exp(-t/Re) ! v(x,y,z,t)= 0.25*(sin(x)cos(y)sin(z)-cos(x)sin(y)cos(z))exp(-t/Re) ! w(x,y,z,t)= 0.5*cos(x)cos(y)sin(z)exp(-t/Re) ! ! .. Parameters .. ! Nx = number of modes in x - power of 2 for FFT ! Ny = number of modes in y - power of 2 for FFT ! Nz = number of modes in z - power of 2 for FFT ! Nt = number of timesteps to take ! Tmax = maximum simulation time ! FFTW_IN_PLACE = value for FFTW input ! FFTW_MEASURE = value for FFTW input ! FFTW_EXHAUSTIVE = value for FFTW input ! FFTW_PATIENT = value for FFTW input ! FFTW_ESTIMATE = value for FFTW input ! FFTW_FORWARD = value for FFTW input ! FFTW_BACKWARD = value for FFTW input ! pi = 3.14159265358979323846264338327950288419716939937510d0 ! Re = Reynolds number ! .. Scalars .. ! i = loop counter in x direction ! j = loop counter in y direction ! k = loop counter in z direction ! n = loop counter for timesteps direction ! allocatestatus = error indicator during allocation ! count = keep track of information written to disk ! iol = size of array to write to disk ! start = variable to record start time of program ! finish = variable to record end time of program ! count_rate = variable for clock count rate ! planfxyz = Forward 3d fft plan ! planbxyz = Backward 3d fft plan ! dt = timestep ! .. Arrays .. ! u = velocity in x direction ! v = velocity in y direction ! w = velocity in z direction ! uold = velocity in x direction at previous timestep ! vold = velocity in y direction at previous timestep ! wold = velocity in z direction at previous timestep ! ux = x derivative of velocity in x direction ! uy = y derivative of velocity in x direction ! uz = z derivative of velocity in x direction ! vx = x derivative of velocity in y direction ! vy = y derivative of velocity in y direction ! vz = z derivative of velocity in y direction ! wx = x derivative of velocity in z direction ! wy = y derivative of velocity in z direction ! wz = z derivative of velocity in z direction ! uxold = x derivative of velocity in x direction ! uyold = y derivative of velocity in x direction ! uzold = z derivative of velocity in x direction ! vxold = x derivative of velocity in y direction ! vyold = y derivative of velocity in y direction ! vzold = z derivative of velocity in y direction ! wxold = x derivative of velocity in z direction ! wyold = y derivative of velocity in z direction ! wzold = z derivative of velocity in z direction ! utemp = temporary storage of u to check convergence ! vtemp = temporary storage of v to check convergence ! wtemp = temporary storage of w to check convergence ! temp_r = temporary storage for untransformed variables ! uhat = Fourier transform of u ! vhat = Fourier transform of v ! what = Fourier transform of w ! rhsuhatfix = Fourier transform of righthand side for u for timestepping ! rhsvhatfix = Fourier transform of righthand side for v for timestepping ! rhswhatfix = Fourier transform of righthand side for w for timestepping ! bx = x component of magnetic field ! by = y component of magnetic field ! bz = z component of magnetic field ! bxold = x component of magnetic field at previous timestep ! byold = y component of magnetic field at previous timestep ! bzold = z component of magnetic field at previous timestep ! bxx = x derivative of x component of magnetic field ! bxy = y derivative of x component of magnetic field ! bxz = z derivative of x component of magnetic field ! byx = x derivative of y component of magnetic field ! byy = y derivative of y component of magnetic field ! byz = z derivative of y component of magnetic field ! bzx = x derivative of z component of magnetic field ! bzy = y derivative of z component of magnetic field ! bzz = z derivative of z component of magnetic field ! bxxold = x derivative of x component of magnetic field ! bxyold = y derivative of x component of magnetic field ! bxzold = z derivative of x component of magnetic field ! byxold = x derivative of y component of magnetic field ! byyold = y derivative of y component of magnetic field ! byzold = z derivative of y component of magnetic field ! bzxold = x derivative of z component of magnetic field ! bzyold = y derivative of z component of magnetic field ! bzzold = z derivative of z component of magnetic field ! bxtemp = temporary storage of bx to check convergence ! bytemp = temporary storage of by to check convergence ! bztemp = temporary storage of bz to check convergence ! bxhat = Fourier transform of bx ! byhat = Fourier transform of by ! bzhat = Fourier transform of bz ! rhsbxhatfix = Fourier transform of righthand side for bx for timestepping ! rhsbyhatfix = Fourier transform of righthand side for by for timestepping ! rhsbzhatfix = Fourier transform of righthand side for bz for timestepping ! nonlinuhat = Fourier transform of nonlinear term for u ! nonlinvhat = Fourier transform of nonlinear term for u ! nonlinwhat = Fourier transform of nonlinear term for u ! phat = Fourier transform of nonlinear term for pressure, p ! temp_c = temporary storage for Fourier transforms ! realtemp = Real storage ! ! .. Vectors .. ! kx = fourier frequencies in x direction ! ky = fourier frequencies in y direction ! kz = fourier frequencies in z direction ! x = x locations ! y = y locations ! z = y locations ! time = times at which save data ! name_config = array to store filename for data to be saved ! ! REFERENCES ! ! ! ACKNOWLEDGEMENTS ! ! ACCURACY ! ! ERROR INDICATORS AND WARNINGS ! ! FURTHER COMMENTS ! ! This program has not been optimized to use the least amount of memory ! but is intended as an example only for which all states can be saved ! !-------------------------------------------------------------------------------- ! External routines required ! ! External libraries required ! 2DECOMP&FFT -- Fast Fourier Transform in the West Library ! (http://2decomp.org/) USE decomp_2d USE decomp_2d_fft USE decomp_2d_io USE MPI IMPLICIT NONE ! declare variables INTEGER(kind=4), PARAMETER :: Nx=64 INTEGER(kind=4), PARAMETER :: Ny=64 INTEGER(kind=4), PARAMETER :: Nz=64 INTEGER(kind=4), PARAMETER :: Lx=1 INTEGER(kind=4), PARAMETER :: Ly=1 INTEGER(kind=4), PARAMETER :: Lz=1 INTEGER(kind=4), PARAMETER :: Nt=5 REAL(kind=8), PARAMETER :: dt=0.05d0/Nt REAL(kind=8), PARAMETER :: Re=1.0d0 REAL(kind=8), PARAMETER :: Rem=1.0d0 REAL(kind=8), PARAMETER :: tol=0.1d0**10 REAL(kind=8), PARAMETER :: theta=0.0d0 REAL(kind=8), PARAMETER & :: pi=3.14159265358979323846264338327950288419716939937510d0 REAL(kind=8), PARAMETER :: ReInv=1.0d0/REAL(Re,kind(0d0)) REAL(kind=8), PARAMETER :: RemInv=1.0d0/REAL(Rem,kind(0d0)) REAL(kind=8), PARAMETER :: dtInv=1.0d0/REAL(dt,kind(0d0)) REAL(kind=8) :: scalemodes,chg,factor REAL(kind=8), DIMENSION(:), ALLOCATABLE :: x, y, z, time,mychg,allchg COMPLEX(kind=8), DIMENSION(:,:,:), ALLOCATABLE :: u, v, w,& ux, uy, uz,& vx, vy, vz,& wx, wy, wz,& uold, uxold, uyold, uzold,& vold, vxold, vyold, vzold,& wold, wxold, wyold, wzold,& utemp, vtemp, wtemp, temp_r, & bx, by, bz,& bxx, bxy, bxz,& byx, byy, byz,& bzx, bzy, bzz,& bxold, bxxold, bxyold, bxzold,& byold, byxold, byyold, byzold,& bzold, bzxold, bzyold, bzzold,& bxtemp, bytemp, bztemp COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE :: kx, ky, kz COMPLEX(kind=8), DIMENSION(:,:,:), ALLOCATABLE :: uhat, vhat, what, & bxhat, byhat, bzhat, & rhsuhatfix, rhsvhatfix, & rhswhatfix, rhsbxhatfix, & rhsbyhatfix, rhsbzhatfix, & nonlinuhat, nonlinvhat, & nonlinwhat, phat, temp_c REAL(kind=8), DIMENSION(:,:,:), ALLOCATABLE :: realtemp ! MPI and 2DECOMP variables TYPE(DECOMP_INFO) :: decomp INTEGER(kind=MPI_OFFSET_KIND) :: filesize, disp INTEGER(kind=4) :: p_row=0, p_col=0, numprocs, myid, ierr ! variables used for saving data and timing INTEGER(kind=4) :: count, iol INTEGER(kind=4) :: i,j,k,n,t,allocatestatus INTEGER(kind=4) :: ind, numberfile CHARACTER*100 :: name_config INTEGER(kind=4) :: start, finish, count_rate ! initialisation of 2DECOMP&FFT CALL MPI_INIT(ierr) CALL MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr) CALL MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr) ! do automatic domain decomposition CALL decomp_2d_init(Nx,Ny,Nz,p_row,p_col) ! get information about domain decomposition choosen CALL decomp_info_init(Nx,Ny,Nz,decomp) ! initialise FFT library CALL decomp_2d_fft_init IF (myid.eq.0) THEN PRINT *,'Grid:',Nx,'X',Ny,'Y',Nz,'Z' PRINT *,'dt:',dt END IF ALLOCATE(x(1:Nx),y(1:Ny),z(1:Nz),time(1:Nt+1),mychg(1:6),allchg(1:6),& u(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& v(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& w(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& ux(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& uy(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& uz(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& vx(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& vy(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& vz(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& wx(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& wy(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& wz(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& uold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& uxold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& uyold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& uzold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& vold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& vxold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& vyold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& vzold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& wold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& wxold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& wyold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& wzold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& utemp(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& vtemp(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& wtemp(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& temp_r(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bx(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& by(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bz(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bxx(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bxy(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bxz(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& byx(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& byy(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& byz(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bzx(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bzy(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bzz(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bxold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bxxold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bxyold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bxzold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& byold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& byxold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& byyold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& byzold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bzold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bzxold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bzyold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bzzold(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bxtemp(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bytemp(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& bztemp(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)),& kx(1:Nx),ky(1:Ny),kz(1:Nz),& uhat(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& vhat(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& what(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& bxhat(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& byhat(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& bzhat(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& rhsuhatfix(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& rhsvhatfix(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& rhswhatfix(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& rhsbxhatfix(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& rhsbyhatfix(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& rhsbzhatfix(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& nonlinuhat(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& nonlinvhat(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& nonlinwhat(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& phat(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& temp_c(decomp%zst(1):decomp%zen(1),& decomp%zst(2):decomp%zen(2),& decomp%zst(3):decomp%zen(3)),& realtemp(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)), stat=AllocateStatus) IF (AllocateStatus .ne. 0) STOP IF (myid.eq.0) THEN PRINT *,'allocated space' END IF ! setup fourier frequencies in x-direction DO i=1,Nx/2+1 kx(i)= cmplx(0.0d0,1.0d0)*REAL(i-1,kind(0d0))/Lx END DO kx(1+Nx/2)=0.0d0 DO i = 1,Nx/2 -1 kx(i+1+Nx/2)=-kx(1-i+Nx/2) END DO ind=1 DO i=-Nx/2,Nx/2-1 x(ind)=2.0d0*pi*REAL(i,kind(0d0))*Lx/REAL(Nx,kind(0d0)) ind=ind+1 END DO ! setup fourier frequencies in y-direction DO j=1,Ny/2+1 ky(j)= cmplx(0.0d0,1.0d0)*REAL(j-1,kind(0d0))/Ly END DO ky(1+Ny/2)=0.0d0 DO j = 1,Ny/2 -1 ky(j+1+Ny/2)=-ky(1-j+Ny/2) END DO ind=1 DO j=-Ny/2,Ny/2-1 y(ind)=2.0d0*pi*REAL(j,kind(0d0))*Ly/REAL(Ny,kind(0d0)) ind=ind+1 END DO ! setup fourier frequencies in z-direction DO k=1,Nz/2+1 kz(k)= cmplx(0.0d0,1.0d0)*REAL(k-1,kind(0d0))/Lz END DO kz(1+Nz/2)=0.0d0 DO k = 1,Nz/2 -1 kz(k+1+Nz/2)=-kz(1-k+Nz/2) END DO ind=1 DO k=-Nz/2,Nz/2-1 z(ind)=2.0d0*pi*REAL(k,kind(0d0))*Lz/REAL(Nz,kind(0d0)) ind=ind+1 END DO scalemodes=1.0d0/REAL(Nx*Ny*Nz,kind(0d0)) IF (myid.eq.0) THEN PRINT *,'Setup grid and fourier frequencies' END IF ! Initial conditions for fluid flow field !initial conditions for Taylor-Green vortex ! factor=2.0d0/sqrt(3.0d0) ! DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) ! u(i,j,k)=factor*sin(theta+2.0d0*pi/3.0d0)*sin(x(i))*cos(y(j))*cos(z(k)) ! END DO; END DO; END DO ! DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) ! v(i,j,k)=factor*sin(theta-2.0d0*pi/3.0d0)*cos(x(i))*sin(y(j))*cos(z(k)) ! END DO ; END DO ; END DO ! DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) ! w(i,j,k)=factor*sin(theta)*cos(x(i))*cos(y(j))*sin(z(k)) ! END DO ; END DO ; END DO time(1)=0.0d0 factor=sqrt(3.0d0) DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) u(i,j,k)=-0.5*( factor*cos(x(i))*sin(y(j))*sin(z(k))& +sin(x(i))*cos(y(j))*cos(z(k)) )*exp(-(factor**2)*time(1)/Re) END DO; END DO; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) v(i,j,k)=0.5*( factor*sin(x(i))*cos(y(j))*sin(z(k))& -cos(x(i))*sin(y(j))*cos(z(k)) )*exp(-(factor**2)*time(1)/Re) END DO ; END DO ; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) w(i,j,k)=cos(x(i))*cos(y(j))*sin(z(k))*exp(-(factor**2)*time(1)/Re) END DO ; END DO ; END DO CALL decomp_2d_fft_3d(u,uhat,DECOMP_2D_FFT_FORWARD) CALL decomp_2d_fft_3d(v,vhat,DECOMP_2D_FFT_FORWARD) CALL decomp_2d_fft_3d(w,what,DECOMP_2D_FFT_FORWARD) ! derivative of u with respect to x, y, and z DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=uhat(i,j,k)*kx(i)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,ux,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=uhat(i,j,k)*ky(j)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,uy,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=uhat(i,j,k)*kz(k)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,uz,DECOMP_2D_FFT_BACKWARD) ! derivative of v with respect to x, y, and z DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=vhat(i,j,k)*kx(i)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,vx,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=vhat(i,j,k)*ky(j)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,vy,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=vhat(i,j,k)*kz(k)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,vz,DECOMP_2D_FFT_BACKWARD) ! derivative of w with respect to x, y, and z DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=what(i,j,k)*kx(i)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,wx,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=what(i,j,k)*ky(j)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,wy,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=what(i,j,k)*kz(k)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,wz,DECOMP_2D_FFT_BACKWARD) ! save initial data n=0 DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) realtemp(i,j,k)=REAL(wy(i,j,k)-vz(i,j,k),KIND=8) END DO ; END DO ; END DO name_config='./data/omegax' CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp) !omegay DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) realtemp(i,j,k)=REAL(uz(i,j,k)-wx(i,j,k),KIND=8) END DO ; END DO ; END DO name_config='./data/omegay' CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp) !omegaz DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) realtemp(i,j,k)=REAL(vx(i,j,k)-uy(i,j,k),KIND=8) END DO ; END DO ; END DO name_config='./data/omegaz' CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp) ! Initial conditions for magnetic field DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) bx(i,j,k)=-0.5*( sin(y(j))*sin(z(k))) END DO; END DO; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) by(i,j,k)=0.5*( sin(x(i))*sin(z(k))) END DO ; END DO ; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) bz(i,j,k)=cos(x(i))*cos(y(j)) END DO ; END DO ; END DO CALL decomp_2d_fft_3d(bx,bxhat,DECOMP_2D_FFT_FORWARD) CALL decomp_2d_fft_3d(by,byhat,DECOMP_2D_FFT_FORWARD) CALL decomp_2d_fft_3d(bz,bzhat,DECOMP_2D_FFT_FORWARD) ! derivative of u with respect to x, y, and z DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=bxhat(i,j,k)*kx(i)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,bxx,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=bxhat(i,j,k)*ky(j)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,bxy,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=bxhat(i,j,k)*kz(k)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,bxz,DECOMP_2D_FFT_BACKWARD) ! derivative of v with respect to x, y, and z DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=byhat(i,j,k)*kx(i)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,byx,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=byhat(i,j,k)*ky(j)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,byy,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=byhat(i,j,k)*kz(k)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,byz,DECOMP_2D_FFT_BACKWARD) ! derivative of w with respect to x, y, and z DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=bzhat(i,j,k)*kx(i)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,wx,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=bzhat(i,j,k)*ky(j)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,wy,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=bzhat(i,j,k)*kz(k)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,bzz,DECOMP_2D_FFT_BACKWARD) ! save initial data n=0 DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) realtemp(i,j,k)=REAL(bzy(i,j,k)-byz(i,j,k),KIND=8) END DO ; END DO ; END DO name_config='./data/omegabx' CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp) !omegay DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) realtemp(i,j,k)=REAL(bxz(i,j,k)-bzx(i,j,k),KIND=8) END DO ; END DO ; END DO name_config='./data/omegaby' CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp) !omegaz DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) realtemp(i,j,k)=REAL(byx(i,j,k)-bxy(i,j,k),KIND=8) END DO ; END DO ; END DO name_config='./data/omegabz' CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp) !start timer CALL system_clock(start,count_rate) DO n=1,Nt !fixed point iteration ! store old fluid field terms DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) uold(i,j,k)=u(i,j,k) uxold(i,j,k)=ux(i,j,k) uyold(i,j,k)=uy(i,j,k) uzold(i,j,k)=uz(i,j,k) END DO ; END DO ; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) vold(i,j,k)=v(i,j,k) vxold(i,j,k)=vx(i,j,k) vyold(i,j,k)=vy(i,j,k) vzold(i,j,k)=vz(i,j,k) END DO ; END DO ; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) wold(i,j,k)=w(i,j,k) wxold(i,j,k)=wx(i,j,k) wyold(i,j,k)=wy(i,j,k) wzold(i,j,k)=wz(i,j,k) END DO ; END DO ; END DO DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) rhsuhatfix(i,j,k) = (dtInv+(0.5*ReInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)))*uhat(i,j,k) END DO ; END DO ; END DO DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) rhsvhatfix(i,j,k) = (dtInv+(0.5*ReInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)))*vhat(i,j,k) END DO ; END DO ; END DO DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) rhswhatfix(i,j,k) = (dtInv+(0.5*ReInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)))*what(i,j,k) END DO ; END DO ; END DO ! Store old Magnetic field terms DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) bxold(i,j,k)=bx(i,j,k) bxxold(i,j,k)=bxx(i,j,k) bxyold(i,j,k)=bxy(i,j,k) bxzold(i,j,k)=bxz(i,j,k) END DO ; END DO ; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) byold(i,j,k)=by(i,j,k) byxold(i,j,k)=byx(i,j,k) byyold(i,j,k)=byy(i,j,k) byzold(i,j,k)=byz(i,j,k) END DO ; END DO ; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) bzold(i,j,k)=bz(i,j,k) bzxold(i,j,k)=bzx(i,j,k) bzyold(i,j,k)=bzy(i,j,k) bzzold(i,j,k)=bzz(i,j,k) END DO ; END DO ; END DO DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) rhsbxhatfix(i,j,k) = (dtInv+(0.5d0*RemInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)))*bxhat(i,j,k) END DO ; END DO ; END DO DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) rhsbyhatfix(i,j,k) = (dtInv+(0.5*RemInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)))*byhat(i,j,k) END DO ; END DO ; END DO DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) rhsbzhatfix(i,j,k) = (dtInv+(0.5*RemInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)))*bzhat(i,j,k) END DO ; END DO ; END DO chg=1 DO WHILE (chg .gt. tol) ! Fluid field DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) temp_r(i,j,k)=0.25d0*((u(i,j,k)+uold(i,j,k))*(ux(i,j,k)+uxold(i,j,k))& +(v(i,j,k)+vold(i,j,k))*(uy(i,j,k)+uyold(i,j,k))& +(w(i,j,k)+wold(i,j,k))*(uz(i,j,k)+uzold(i,j,k))& -(bx(i,j,k)+bxold(i,j,k))*(bxx(i,j,k)+bxxold(i,j,k))& -(by(i,j,k)+byold(i,j,k))*(bxy(i,j,k)+bxyold(i,j,k))& -(bz(i,j,k)+bzold(i,j,k))*(bxz(i,j,k)+bxzold(i,j,k))) END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_r,nonlinuhat,DECOMP_2D_FFT_FORWARD) DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) temp_r(i,j,k)=0.25d0*((u(i,j,k)+uold(i,j,k))*(vx(i,j,k)+vxold(i,j,k))& +(v(i,j,k)+vold(i,j,k))*(vy(i,j,k)+vyold(i,j,k))& +(w(i,j,k)+wold(i,j,k))*(vz(i,j,k)+vzold(i,j,k))& -(bx(i,j,k)+bxold(i,j,k))*(byx(i,j,k)+byxold(i,j,k))& -(by(i,j,k)+byold(i,j,k))*(byy(i,j,k)+byyold(i,j,k))& -(bz(i,j,k)+bzold(i,j,k))*(byz(i,j,k)+byzold(i,j,k))) END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_r,nonlinvhat,DECOMP_2D_FFT_FORWARD) DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) temp_r(i,j,k)=0.25d0*((u(i,j,k)+uold(i,j,k))*(wx(i,j,k)+wxold(i,j,k))& +(v(i,j,k)+vold(i,j,k))*(wy(i,j,k)+wyold(i,j,k))& +(w(i,j,k)+wold(i,j,k))*(wz(i,j,k)+wzold(i,j,k))& -(bx(i,j,k)+bxold(i,j,k))*(bzx(i,j,k)+bzxold(i,j,k))& -(by(i,j,k)+byold(i,j,k))*(bzy(i,j,k)+bzyold(i,j,k))& -(bz(i,j,k)+bzold(i,j,k))*(bzz(i,j,k)+bzzold(i,j,k))) END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_r,nonlinwhat,DECOMP_2D_FFT_FORWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) phat(i,j,k)=-1.0d0*( kx(i)*nonlinuhat(i,j,k)& +ky(j)*nonlinvhat(i,j,k)& +kz(k)*nonlinwhat(i,j,k))& /(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)+0.1d0**13) END DO ; END DO ; END DO DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) uhat(i,j,k)=(rhsuhatfix(i,j,k)-nonlinuhat(i,j,k)-kx(i)*phat(i,j,k))/& (dtInv-(0.5d0*ReInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k))) !*scalemodes END DO ; END DO ; END DO DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) vhat(i,j,k)=(rhsvhatfix(i,j,k)-nonlinvhat(i,j,k)-ky(j)*phat(i,j,k))/& (dtInv-(0.5d0*ReInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k))) !*scalemodes END DO ; END DO ; END DO DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) what(i,j,k)=(rhswhatfix(i,j,k)-nonlinwhat(i,j,k)-kz(k)*phat(i,j,k))/& (dtInv-(0.5d0*ReInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k))) !*scalemodes END DO ; END DO ; END DO ! derivative of u with respect to x, y, and z DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=uhat(i,j,k)*kx(i)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,ux,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3); DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=uhat(i,j,k)*ky(j)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,uy,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3); DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=uhat(i,j,k)*kz(k)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,uz,DECOMP_2D_FFT_BACKWARD) ! derivative of v with respect to x, y, and z DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=vhat(i,j,k)*kx(i)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,vx,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3); DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=vhat(i,j,k)*ky(j)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,vy,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=vhat(i,j,k)*kz(k)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,vz,DECOMP_2D_FFT_BACKWARD) ! derivative of w with respect to x, y, and z DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=what(i,j,k)*kx(i)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,wx,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=what(i,j,k)*ky(j)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,wy,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=what(i,j,k)*kz(k)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,wz,DECOMP_2D_FFT_BACKWARD) DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) utemp(i,j,k)=u(i,j,k) END DO ; END DO ; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) vtemp(i,j,k)=v(i,j,k) END DO ; END DO ; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) wtemp(i,j,k)=w(i,j,k) END DO ; END DO ; END DO CALL decomp_2d_fft_3d(uhat,u,DECOMP_2D_FFT_BACKWARD) CALL decomp_2d_fft_3d(vhat,v,DECOMP_2D_FFT_BACKWARD) CALL decomp_2d_fft_3d(what,w,DECOMP_2D_FFT_BACKWARD) DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) u(i,j,k)=u(i,j,k)*scalemodes END DO ; END DO ; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) v(i,j,k)=v(i,j,k)*scalemodes END DO ; END DO ; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) w(i,j,k)=w(i,j,k)*scalemodes END DO ; END DO ; END DO ! Magnetic field DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) temp_r(i,j,k)=0.25d0*((u(i,j,k)+uold(i,j,k))*(bxx(i,j,k)+bxxold(i,j,k))& +(v(i,j,k)+vold(i,j,k))*(bxy(i,j,k)+bxyold(i,j,k))& +(w(i,j,k)+wold(i,j,k))*(bxz(i,j,k)+bxzold(i,j,k))& -(bx(i,j,k)+bxold(i,j,k))*(ux(i,j,k)+uxold(i,j,k))& -(by(i,j,k)+byold(i,j,k))*(uy(i,j,k)+uyold(i,j,k))& -(bz(i,j,k)+bzold(i,j,k))*(uz(i,j,k)+uzold(i,j,k))) END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_r,nonlinuhat,DECOMP_2D_FFT_FORWARD) DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) temp_r(i,j,k)=0.25d0*((u(i,j,k)+uold(i,j,k))*(byx(i,j,k)+byxold(i,j,k))& +(v(i,j,k)+vold(i,j,k))*(byy(i,j,k)+byyold(i,j,k))& +(w(i,j,k)+wold(i,j,k))*(byz(i,j,k)+byzold(i,j,k))& -(bx(i,j,k)+bxold(i,j,k))*(vx(i,j,k)+vxold(i,j,k))& -(by(i,j,k)+byold(i,j,k))*(vy(i,j,k)+vyold(i,j,k))& -(bz(i,j,k)+bzold(i,j,k))*(vz(i,j,k)+vzold(i,j,k))) END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_r,nonlinvhat,DECOMP_2D_FFT_FORWARD) DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) temp_r(i,j,k)=0.25d0*((u(i,j,k)+uold(i,j,k))*(bzx(i,j,k)+bzxold(i,j,k))& +(v(i,j,k)+vold(i,j,k))*(bzy(i,j,k)+bzyold(i,j,k))& +(w(i,j,k)+wold(i,j,k))*(bzz(i,j,k)+bzzold(i,j,k))& -(bx(i,j,k)+bxold(i,j,k))*(wx(i,j,k)+wxold(i,j,k))& -(by(i,j,k)+byold(i,j,k))*(wy(i,j,k)+wyold(i,j,k))& -(bz(i,j,k)+bzold(i,j,k))*(wz(i,j,k)+wzold(i,j,k))) END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_r,nonlinwhat,DECOMP_2D_FFT_FORWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) phat(i,j,k)=-1.0d0*( kx(i)*nonlinuhat(i,j,k)& +ky(j)*nonlinvhat(i,j,k)& +kz(k)*nonlinwhat(i,j,k))& /(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)+0.1d0**13) END DO ; END DO ; END DO DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) bxhat(i,j,k)=(rhsbxhatfix(i,j,k)-nonlinuhat(i,j,k)-kx(i)*phat(i,j,k))/& (dtInv-(0.5d0*RemInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k))) !*scalemodes END DO ; END DO ; END DO DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) byhat(i,j,k)=(rhsbyhatfix(i,j,k)-nonlinvhat(i,j,k)-ky(j)*phat(i,j,k))/& (dtInv-(0.5d0*RemInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k))) !*scalemodes END DO ; END DO ; END DO DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) bzhat(i,j,k)=(rhsbzhatfix(i,j,k)-nonlinwhat(i,j,k)-kz(k)*phat(i,j,k))/& (dtInv-(0.5d0*RemInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k))) !*scalemodes END DO ; END DO ; END DO ! derivative of u with respect to x, y, and z DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=bxhat(i,j,k)*kx(i)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,bxx,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3); DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=bxhat(i,j,k)*ky(j)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,bxy,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3); DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=uhat(i,j,k)*kz(k)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,bxz,DECOMP_2D_FFT_BACKWARD) ! derivative of v with respect to x, y, and z DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=byhat(i,j,k)*kx(i)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,byx,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3); DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=byhat(i,j,k)*ky(j)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,byy,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=byhat(i,j,k)*kz(k)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,byz,DECOMP_2D_FFT_BACKWARD) ! derivative of w with respect to x, y, and z DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=bzhat(i,j,k)*kx(i)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,bzx,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=bzhat(i,j,k)*ky(j)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,bzy,DECOMP_2D_FFT_BACKWARD) DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1) temp_c(i,j,k)=bzhat(i,j,k)*kz(k)*scalemodes END DO ; END DO ; END DO CALL decomp_2d_fft_3d(temp_c,bzz,DECOMP_2D_FFT_BACKWARD) DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) bxtemp(i,j,k)=bx(i,j,k) END DO ; END DO ; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) bytemp(i,j,k)=by(i,j,k) END DO ; END DO ; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) bztemp(i,j,k)=bz(i,j,k) END DO ; END DO ; END DO CALL decomp_2d_fft_3d(bxhat,bx,DECOMP_2D_FFT_BACKWARD) CALL decomp_2d_fft_3d(byhat,by,DECOMP_2D_FFT_BACKWARD) CALL decomp_2d_fft_3d(bzhat,bz,DECOMP_2D_FFT_BACKWARD) DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) bx(i,j,k)=bx(i,j,k)*scalemodes END DO ; END DO ; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) by(i,j,k)=by(i,j,k)*scalemodes END DO ; END DO ; END DO DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) bz(i,j,k)=bz(i,j,k)*scalemodes END DO ; END DO ; END DO mychg(1) =maxval(abs(utemp-u)) mychg(2) =maxval(abs(vtemp-v)) mychg(3) =maxval(abs(wtemp-w)) mychg(4) =maxval(abs(bxtemp-bx)) mychg(5) =maxval(abs(bytemp-by)) mychg(6) =maxval(abs(bztemp-bz)) CALL MPI_ALLREDUCE(mychg,allchg,6,MPI_DOUBLE_PRECISION,MPI_MAX,MPI_COMM_WORLD,ierr) chg=allchg(1)+allchg(2)+allchg(3)+allchg(4)+allchg(5)+allchg(6) IF (myid.eq.0) THEN PRINT *,'chg:',chg END IF END DO time(n+1)=n*dt IF (myid.eq.0) THEN PRINT *,'time',n*dt END IF !save omegax, omegay, omegaz, omegabx, omegaby, and omegabz !omegax DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) realtemp(i,j,k)=REAL(wy(i,j,k)-vz(i,j,k),KIND=8) END DO ; END DO ; END DO name_config='./data/omegax' CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp) !omegay DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) realtemp(i,j,k)=REAL(uz(i,j,k)-wx(i,j,k),KIND=8) END DO ; END DO ; END DO name_config='./data/omegay' CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp) !omegaz DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) realtemp(i,j,k)=REAL(vx(i,j,k)-uy(i,j,k),KIND=8) END DO ; END DO ; END DO name_config='./data/omegaz' CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp) !omegabx DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) realtemp(i,j,k)=REAL(bzy(i,j,k)-byz(i,j,k),KIND=8) END DO ; END DO ; END DO name_config='./data/omegabx' CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp) !omegaby DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) realtemp(i,j,k)=REAL(bxz(i,j,k)-bzx(i,j,k),KIND=8) END DO ; END DO ; END DO name_config='./data/omegaby' CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp) !omegabz DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1) realtemp(i,j,k)=REAL(byx(i,j,k)-bxy(i,j,k),KIND=8) END DO ; END DO ; END DO name_config='./data/omegabz' CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp) END DO CALL system_clock(finish,count_rate) IF (myid.eq.0) then PRINT *, 'Program took', REAL(finish-start)/REAL(count_rate), 'for main timestepping loop' END IF IF (myid.eq.0) THEN name_config = './data/tdata.dat' OPEN(unit=11,FILE=name_config,status="UNKNOWN") REWIND(11) DO n=1,1+Nt WRITE(11,*) time(n) END DO CLOSE(11) name_config = './data/xcoord.dat' OPEN(unit=11,FILE=name_config,status="UNKNOWN") REWIND(11) DO i=1,Nx WRITE(11,*) x(i) END DO CLOSE(11) name_config = './data/ycoord.dat' OPEN(unit=11,FILE=name_config,status="UNKNOWN") REWIND(11) DO j=1,Ny WRITE(11,*) y(j) END DO CLOSE(11) name_config = './data/zcoord.dat' OPEN(unit=11,FILE=name_config,status="UNKNOWN") REWIND(11) DO k=1,Nz WRITE(11,*) z(k) END DO CLOSE(11) PRINT *,'Saved data' END IF ! clean up CALL decomp_2d_fft_finalize CALL decomp_2d_finalize DEALLOCATE(x,y,z,time,mychg,allchg,& u,v,w,ux,uy,uz,vx,vy,vz,& wx,wy,wz,uold,uxold,& uyold,uzold,vold,vxold,& vyold,vzold,wold,wxold,& wyold,wzold,utemp,vtemp,& wtemp,temp_r,& bx,by,bz,bxx,bxy,bxz,& byx,byy,byz,bzx,bzy,bzz,& bxold,bxxold,bxyold,bxzold,& byold,byxold,byyold,byzold,& bzold,bzxold,bzyold,bzzold,& bxtemp,bytemp,bztemp,kx,ky,kz,& uhat,vhat,what,bxhat,& byhat,bzhat,rhsuhatfix,rhsvhatfix,& rhswhatfix,rhsbxhatfix,rhsbyhatfix,rhsbzhatfix,& nonlinuhat,nonlinvhat,nonlinwhat,phat,& temp_c,realtemp,stat=AllocateStatus) IF (AllocateStatus .ne. 0) STOP IF (myid.eq.0) THEN PRINT *,'Program execution complete' END IF CALL MPI_FINALIZE(ierr) END PROGRAM main F A subroutine to save real array data for the parallel MPI Fortran program to solve the 3D incompressible magnetohydrodynamics equations Code download SUBROUTINE savedata(Nx,Ny,Nz,plotnum,name_config,field,decomp) !-------------------------------------------------------------------- ! ! ! PURPOSE ! ! This subroutine saves a three dimensional real array in binary ! format ! ! INPUT ! ! .. Scalars .. ! Nx = number of modes in x - power of 2 for FFT ! Ny = number of modes in y - power of 2 for FFT ! Nz = number of modes in z - power of 2 for FFT ! plotnum = number of plot to be made ! .. Arrays .. ! field = real data to be saved ! name_config = root of filename to save to ! ! .. Output .. ! plotnum = number of plot to be saved ! .. Special Structures .. ! decomp = contains information on domain decomposition ! see http://www.2decomp.org/ for more information ! LOCAL VARIABLES ! ! .. Scalars .. ! i = loop counter in x direction ! j = loop counter in y direction ! k = loop counter in z direction ! count = counter ! iol = size of file ! .. Arrays .. ! number_file = array to hold the number of the plot ! ! REFERENCES ! ! ACKNOWLEDGEMENTS ! ! ACCURACY ! ! ERROR INDICATORS AND WARNINGS ! ! FURTHER COMMENTS !-------------------------------------------------------------------- ! External routines required ! ! External libraries required ! 2DECOMP&FFT -- Domain decomposition and Fast Fourier Library ! (http://www.2decomp.org/index.html) ! MPI library USE decomp_2d USE decomp_2d_fft USE decomp_2d_io IMPLICIT NONE INCLUDE 'mpif.h' ! Declare variables INTEGER(KIND=4), INTENT(IN) :: Nx,Ny,Nz INTEGER(KIND=4), INTENT(IN) :: plotnum TYPE(DECOMP_INFO), INTENT(IN) :: decomp REAL(KIND=8), DIMENSION(decomp%xst(1):decomp%xen(1),& decomp%xst(2):decomp%xen(2),& decomp%xst(3):decomp%xen(3)), & INTENT(IN) :: field CHARACTER*100, INTENT(IN) :: name_config INTEGER(kind=4) :: i,j,k,iol,count,ind CHARACTER*100 :: number_file ! create character array with full filename ind = index(name_config,' ') - 1 WRITE(number_file,'(i0)') 10000000+plotnum number_file = name_config(1:ind)//number_file ind = index(number_file,' ') - 1 number_file = number_file(1:ind)//'.datbin' CALL decomp_2d_write_one(1,field,number_file) END SUBROUTINE savedata G A makefile to compile the program which finds a numerical solution to the 3D magnetohydrodynamic equations Code download COMPILER = mpif90 decompdir= ../2decomp_fft FLAGS = -O0 DECOMPLIB = -I${decompdir}/include -L${decompdir}/lib -l2decomp_fft LIBS = #-L${FFTW_LINK} -lfftw3 -lm SOURCES = Magnetohydrodynamics3DfftIMR.f90 savedata.f90 ns3d: $(SOURCES) ${COMPILER} -o mhd3d $(FLAGS) $(SOURCES) $(LIBS) $(DECOMPLIB) clean: rm -f *.o rm -f *.mod clobber: rm -f mhd3d Other Publicly Available Fourier Spectral Programs for Solving the Magnetohydrodynamics Equations One other publicly available Fast Fourier Transform based program for solving the incompressible magnetohydrodynamics equations is: Tarang available from http://turbulencehub.org/index.php/codes/tarang/ Exercises The three dimensional program we have written is not the most efficient since one can use a real to complex transform to halve the work done. Implement a real to complex transform in one of the Navier-Stokes programs. The programs we have written can also introduce some aliasing errors. By reading a book on spectral methods, such as Canuto et al.Canuto et al (2006)Canuto et al (2007), find out what aliasing errors are. Explain why the strategy explained in JohnstoneJohnstone (2012) can reduce aliasing errors. See if you can find other open source Fast Fourier Transform based solvers for the incompressible magnetohydrodynamics equations and add them to the list above. There are suggestions that the incompressible magnetohydrodynamics equations can exhibit singular behavior by having one of the norms of the derivatives of the magnetic or velocity field becoming infinite. Experiment with simulations to see if you can find evidence for this. One of the first papers reporting computational experiments in this area is Orszag and TangOrszag and Tang (1979). Translate one of the programs above to another language. The following site may help, http://rosettacode.org/wiki/Rosetta_Code The full magnetohydrodynamic equations are : 12 : 13 : 14 : 15 : 16 : 17 where is the velocity, is the pressure, is the current density, is the magnetic field and is the electric field. Simulate these full equations in two and three dimensions - it is at present unknown whether solutions to these equations exist for all time in three dimensions see Germain, Ibrahim and Masmoudi Germain, Ibrahim and Masmoudi (2012), though Masmoudi Masmoudi (2010) shows that some solutions are well behaved in two dimensions. Analysis of the reduced equations is in Duvaut and Lions Duvaut and Lions (1972) and can also be found in Gerbeau, Le Bris, and Lelièvre Gerbeau, Le Bris, and Lelièvre (2006). An example program is in listing H H A Matlab program which finds a numerical solution to the 3D full magnetohydrodynamic equations Code download % A program to solve the Full 3D incompressible magnetohydrodynamic equations % with periodic boundary % conditions. The program is based on the Orszag-Patterson algorithm as % documented on pg. 98 of C. Canuto, M.Y. Hussaini, A. Quarteroni and % T.A. Zhang "Spectral Methods: Evolution to Complex Geometries and % Applications to Fluid Dynamics" Springer (2007) % % The Helmholtz decomposition is used to project the magnetic field onto a % divergence free subspace. Initial work on this has been done with Damian % San Roman Alerigi % % For plotting, the program uses vol3d, which is available at: % % http://www.mathworks.com/matlabcentral/fileexchange/22940-vol3d-v2 % clear all; format compact; format short; set(0,'defaultaxesfontsize',30,'defaultaxeslinewidth',.7,... 'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,... 'defaultaxesfontweight','bold') % set up grid tic Lx = 1; % period 2*pi*L Ly = 1; % period 2*pi*L Lz = 1; % period 2*pi*L Nx = 32; % number of harmonics Ny = 32; % number of harmonics Nz = 32; % number of harmonics Nt = 50; % number of time slices dt = 1/Nt; % time step t=0; % initial time Re = 1.0; % Reynolds number Rem = 1.0; % Magnetic Reynolds number tol=10^(-10); % initialise variables x = (2*pi/Nx)*(-Nx/2:Nx/2 -1)'*Lx; % x coordinate kx = 1i*[0:Nx/2-1 0 -Nx/2+1:-1]'/Lx; % wave vector y = (2*pi/Ny)*(-Ny/2:Ny/2 -1)'*Ly; % y coordinate ky = 1i*[0:Ny/2-1 0 -Ny/2+1:-1]'/Ly; % wave vector z = (2*pi/Nz)*(-Nz/2:Nz/2 -1)'*Lz; % y coordinate kz = 1i*[0:Nz/2-1 0 -Nz/2+1:-1]'/Lz; % wave vector [xx,yy,zz]=meshgrid(x,y,z); [kxm,kym,kzm]=meshgrid(kx,ky,kz); [k2xm,k2ym,k2zm]=meshgrid(kx.^2,ky.^2,kz.^2); % initial conditions for Taylor-Green vortex % theta=0; % u=(2/sqrt(3))*sin(theta+2*pi/3)*sin(xx).*cos(yy).*cos(zz); % v=(2/sqrt(3))*sin(theta-2*pi/3)*cos(xx).*sin(yy).*cos(zz); % w=(2/sqrt(3))*sin(theta)*cos(xx).*cos(yy).*sin(zz); % initial condition sl=1; sk=1; sm=1; lamlkm=sqrt(sl.^2+sk.^2+sm.^2); u=-0.5*(lamlkm*sl*cos(sk*xx).*sin(sl*yy).*sin(sm.*zz)... +sm*sk*sin(sk*xx).*cos(sl*yy).*cos(sm.*zz))... .*exp(-t*(lamlkm^2)/Re); v=0.5*(lamlkm*sk*sin(sk*xx).*cos(sl*yy).*sin(sm.*zz)... -sm*sl*cos(sk*xx).*sin(sl*yy).*cos(sm.*zz))... *exp(-t*(lamlkm^2)/Re); w=cos(sk*xx).*cos(sl*yy).*sin(sm*zz)*exp(-t*(lamlkm^2)/Re); uhat=fftn(u); vhat=fftn(v); what=fftn(w); ux=ifftn(uhat.*kxm);uy=ifftn(uhat.*kym);uz=ifftn(uhat.*kzm); vx=ifftn(vhat.*kxm);vy=ifftn(vhat.*kym);vz=ifftn(vhat.*kzm); wx=ifftn(what.*kxm);wy=ifftn(what.*kym);wz=ifftn(what.*kzm); % Magnetic field bx=sin(sl*yy).*sin(sm.*zz); by=sin(sm.*zz); bz=cos(sk*xx).*cos(sl*yy); bxhat=fftn(bx); byhat=fftn(by); bzhat=fftn(bz); bxx=ifftn(bxhat.*kxm);bxy=ifftn(bxhat.*kym);bxz=ifftn(bxhat.*kzm); byx=ifftn(byhat.*kxm);byy=ifftn(byhat.*kym);byz=ifftn(byhat.*kzm); bzx=ifftn(bzhat.*kxm);bzy=ifftn(bzhat.*kym);bzz=ifftn(bzhat.*kzm); % Electric field ex=sin(sl*yy).*sin(sm.*zz); ey=sin(sm.*zz); ez=cos(sk*xx).*cos(sl*yy); exhat=fftn(ex); eyhat=fftn(ey); ezhat=fftn(ez); exx=ifftn(exhat.*kxm);exy=ifftn(exhat.*kym);exz=ifftn(exhat.*kzm); eyx=ifftn(eyhat.*kxm);eyy=ifftn(eyhat.*kym);eyz=ifftn(eyhat.*kzm); ezx=ifftn(ezhat.*kxm);ezy=ifftn(ezhat.*kym);ezz=ifftn(ezhat.*kzm); % Current density jx=ex+uy.*bz-uz.*by; jy=ey+uz.*bx-ux.*bz; jz=ez+ux.*by-uy.*bx; jxhat=fftn(jx); jyhat=fftn(jy); jzhat=fftn(jz); % calculate fluid vorticity for plotting omegax=wy-vz; omegay=uz-wx; omegaz=vx-uy; omegatot=omegax.^2+omegay.^2+omegaz.^2; % calculate magnetic vorticity for plotting omegabx=bxy-byz; omegaby=bxz-bzx; omegabz=byx-bxy; omegabtot=omegabx.^2+omegaby.^2+omegabz.^2; figure(1); clf; n=0; subplot(2,1,1); title(['|Fluid omega|^2 ',num2str(t)]); H = vol3d('CData',omegatot,'texture','3D','XData',x,'YData',y,'ZData',z); xlabel('x'); ylabel('y'); zlabel('z'); colorbar; axis equal; axis square; view(3); subplot(2,1,2); title(['|Magnetic omega|^2 ',num2str(t)]); H = vol3d('CData',omegabtot,'texture','3D','XData',x,'YData',y,'ZData',z); xlabel('x'); ylabel('y'); zlabel('z'); colorbar; axis equal; axis square; view(3); for n=1:Nt uold=u; uxold=ux; uyold=uy; uzold=uz; vold=v; vxold=vx; vyold=vy; vzold=vz; wold=w; wxold=wx; wyold=wy; wzold=wz; bxold=bx; byold=by; bzold=bz; exold=ex; exxold=exx; exyold=exy; exzold=exz; eyold=ey; eyxold=eyx; eyyold=eyy; eyzold=eyz; ezold=ez; ezxold=ezx; ezyold=ezy; ezzold=ezz; jxold=jx; jyold=jy; jzold=jz; rhsuhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*uhat; rhsvhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*vhat; rhswhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*what; rhsbxhatfix=bxhat-0.5*dt*(kym.*ezhat-kzm.*eyhat); rhsbyhatfix=byhat-0.5*dt*(kzm.*exhat-kxm.*ezhat); rhsbzhatfix=bzhat-0.5*dt*(kxm.*eyhat-kym.*exhat); rhsexhatfix=exhat+0.5*dt*(kym.*bzhat-kzm.*byhat)-0.5*dt*jxhat; rhseyhatfix=eyhat+0.5*dt*(kzm.*bxhat-kxm.*bzhat)-0.5*dt*jyhat; rhsezhatfix=ezhat+0.5*dt*(kxm.*byhat-kym.*bxhat)-0.5*dt*jzhat; chg=1; t=t+dt; while (chg>tol) nonlinu=0.25*( (u+uold).*(ux+uxold)... +(v+vold).*(uy+uyold)... +(w+wold).*(uz+uzold)... -(jy+jyold).*(bz+bzold)... +(jz+jzold).*(by+byold)); nonlinv=0.25*( (u+uold).*(vx+vxold)... +(v+vold).*(vy+vyold)... +(w+wold).*(vz+vzold)... +(jx+jxold).*(bz+bzold)... -(jz+jzold).*(bx+bxold)); nonlinw=0.25*( (u+uold).*(wx+wxold)... +(v+vold).*(wy+wyold)... +(w+wold).*(wz+wzold)... -(jx+jxold).*(by+byold)... +(jy+jyold).*(bx+bxold)); nonlinuhat=fftn(nonlinu); nonlinvhat=fftn(nonlinv); nonlinwhat=fftn(nonlinw); phat=-1.0*(kxm.*nonlinuhat+kym.*nonlinvhat+kzm.*nonlinwhat)... ./(k2xm+k2ym+k2zm+0.1^13); uhat=(rhsuhatfix-nonlinuhat-kxm.*phat)... ./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm)); vhat=(rhsvhatfix-nonlinvhat-kym.*phat)... ./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm)); what=(rhswhatfix-nonlinwhat-kzm.*phat)... ./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm)); ux=ifftn(uhat.*kxm);uy=ifftn(uhat.*kym);uz=ifftn(uhat.*kzm); vx=ifftn(vhat.*kxm);vy=ifftn(vhat.*kym);vz=ifftn(vhat.*kzm); wx=ifftn(what.*kxm);wy=ifftn(what.*kym);wz=ifftn(what.*kzm); utemp=u; vtemp=v; wtemp=w; u=ifftn(uhat); v=ifftn(vhat); w=ifftn(what); bxhat=rhsbxhatfix-0.5*dt*(kym.*ezhat-kzm.*eyhat); byhat=rhsbyhatfix-0.5*dt*(kzm.*exhat-kxm.*ezhat); bzhat=rhsbzhatfix-0.5*dt*(kxm.*eyhat-kym.*exhat); bxtemp=bx; bytemp=by; bztemp=bz; bx=ifftn(bxhat); by=ifftn(byhat); bz=ifftn(bzhat); exhat=rhsexhatfix+0.5*dt*(kym.*bzhat-kzm.*byhat-jxhat); eyhat=rhseyhatfix+0.5*dt*(kzm.*bxhat-kxm.*bzhat-jyhat); ezhat=rhsezhatfix+0.5*dt*(kxm.*byhat-kym.*bxhat-jzhat); extemp=ex; eytemp=ey; eztemp=ez; ex=ifftn(exhat); ey=ifftn(eyhat); ez=ifftn(ezhat); jx=ex+uy.*bz-uz.*by; jy=ey+uz.*bx-ux.*bz; jz=ez+ux.*by-uy.*bx; jxhat=fftn(jx); jyhat=fftn(jy); jzhat=fftn(jz); chg=max(abs(utemp-u)) +max(abs(vtemp-v)) +max(abs(wtemp-w))+... max(abs(bxtemp-bx)) +max(abs(bytemp-by))+max(abs(bztemp-bz))+... max(abs(extemp-ex)) +max(abs(eytemp-ey))+max(abs(eztemp-ez)); end % calculate vorticity for plotting omegax=wy-vz; omegay=uz-wx; omegaz=vx-uy; omegatot=omegax.^2+omegay.^2+omegaz.^2; % calculate magnetic vorticity for plotting omegabx=bxy-byz; omegaby=bxz-bzx; omegabz=byx-bxy; omegabtot=omegabx.^2+omegaby.^2+omegabz.^2; figure(1); clf; subplot(2,1,1); title(['|Fluid omega|^2 ',num2str(t)]); H = vol3d('CData',omegatot,'texture','3D','XData',x,'YData',y,'ZData',z); xlabel('x'); ylabel('y'); zlabel('z'); colorbar; axis equal; axis square; view(3); subplot(2,1,2); title(['|Magnetic omega|^2 ',num2str(t)]); H = vol3d('CData',omegabtot,'texture','3D','XData',x,'YData',y,'ZData',z); xlabel('x'); ylabel('y'); zlabel('z'); colorbar; axis equal; axis square; view(3); end A compressible magnetohydrodynamic model is : 18 : 19 : 20 : 21 : 22 : 23 : 24 : 25 : 26 : 27 : 28 where is the velocity, is the pressure, is the current density, is the magnetic field, is the electric field, is the energy density, is the viscous fluid stress, is the temperature field, is the constant volume specific heat capacity, is the ideal gas constant, is the heat flux. Simulate these full equations in two and three dimensions. Notes References CanHusQuaZan06 Canuto C. Hussaini M.Y. Quarteroni A. Zang T.A. Spectral Methods: Fundamentals in Single Domains Springer 2006 CanHusQuaZan07 Canuto C. Hussaini M.Y. Quarteroni A. Zang T.A. Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics Springer 2007 CarPou09 Carbone V. Pouquet A. Turbulence in Space Plasmas :https://sites.google.com/site/pouquetannick/home An Introduction to Fluid and MHD Turbulence for Astrophysical Flows: Theory, Observational and Numerical Data, and Modeling Lecture Notes in Physics 778 71–128 Springer-Verlag 2009 10.1007/978-3-642-00210-6_3 978-3-642-00209-0 Cha61 Chandrasekhar Subrahmanyan :https://en.wikipedia.org/wiki/Subrahmanyan_Chandrasekhar Hydrodynamic and Hyrdomagnetic Stability Dover 1961 0-486-64071-X Dav01 Davidson P.A. :http://www3.eng.cam.ac.uk/~pad3/ An Introduction to Magnetohydrodynamics Cambridge University Press 2001 10.1017/CBO9780511626333 9780521794879 Dor07 Dorch Soren http://www.astro.ku.dk/~dorch/ Magnetohydrodynamic Scholarpedia 2 4 2295 2007 10.4249/scholarpedia.2295 free DuvLio72 Duvaut G. Lions J.L. Inéquations en thermoélasticité et magnétohydrodynamique Archive for Rational Mechanics and Analysis 46 4 241–279 1972 10.1007/BF00250512 120571424 GerLeBLel06 Gerbeau J.-F. :https://team.inria.fr/reo/team-members/gerbeau/ Le Bris C. :http://cermics.enpc.fr/~lebris/ Lelièvre T. :http://cermics.enpc.fr/~lelievre/ Mathematical Methods for the Magnetohydrodynamics of Liquid Metals Oxford University Press 2006 10.1093/acprof:oso/9780198566656.001.0001 978-0-19-856665-6 GerIbrMas12 Germain P. :http://www.cims.nyu.edu/~pgermain/indexenglish.html Ibrahim S. :http://www.math.uvic.ca/~ibrahim/ Masmoudi N. :http://www.math.nyu.edu/faculty/masmoudi/ Well Posedness of the Navier-Stokes-Maxwell Equations 2012 math.AP 1207.6187 Joh12 Johnstone R. Improved Scaling for Direct Numerical Simulations of Turbulence HECTOR distributed Computational Science and Engineering Report 2012 http://www.hector.ac.uk/cse/distributedcse/reports/ss3f-swt/ Mas10 Masmoudi N. :http://www.math.nyu.edu/faculty/masmoudi/ Global well posedness for the Maxwell-Navier-Stokes system in 2D J. Math. Pures. Appl. 93 559–571 2010 http://www.math.nyu.edu/faculty/masmoudi/Maxwell-NS-2D.pdf 10.1017/S002211207900210X 123354420 OrsTan79 Orszag S. Tang C.M. Small-scale structure of two-dimensional magnetohydrodynamic turbulence Journal of Fluid Mechanics 90 1 129–143 1979 10.1017/S002211207900210X 123354420 Sch09 Schnack D.D. Lectures in Magnetohydrodynamics Springer Lecture Notes in Physics 780 2009 http://www.physics.wisc.edu/grads/courses/726-f07/lecture_notes.html 10.1007/978-3-642-00688-3 978-3-642-00687-6 SerTem83 Sermange M. Temam R. :http://mypage.iu.edu/~temam/ Some mathematical questions related to the MHD equations Communications on Pure and Applied Mathematics 36 5 635–664 1983 https://hal.inria.fr/inria-00076373/file/RR-0185.pdf 10.1002/cpa.3160360506 Ver04 Verma Mahendra :http://home.iitk.ac.in/~mkv/Site/Welcome.html Statistical theory of magnetohydrodynamic turbulence: recent results Physics Reports 401 5–6 220–380 2004 10.1016/j.physrep.2004.07.007 nlin/0404043 119352240 WikipediaMHD Magnetohydrodynamics https://en.wikipedia.org/wiki/Magnetohydrodynamics Zhe12 Zheng, Linjin Topics in Magnetohydrodynamics InTech 2009 http://www.intechopen.com/books/topics-in-magnetohydrodynamics 10.5772/2080 978-953-51-0211-3