Hi everyone 
I currently try to simply a programm calculating low-Prandl number convection using derived types. For the problem I need to take the inputs variables from a external file which I do with a subroutine read_inputs and namelist. I however now unable to create an executable even though the compiler does not indicate an error. I get the following message on the terminal:
/usr/bin/ld:lowPrandt_parameters_FULL.txt: file format not recognized; treating as linker script
/usr/bin/ld:lowPrandt_parameters_FULL.txt:2: syntax error
collect2: error: ld returned 1 exit status
./run: line 4: ./a.out: No such file or directory
Where I was executing the following commands on linux:
gfortran -o ex1 -fbacktrace -g -fcheck=all -fbounds-check -finit-real=snan -ffpe-trap=invalid,zero,overflow -O0 -O -Wall -g modul_Low_Prandl.f90 exercise_11.f90 lowPrandt_parameters_FULL.txt
./a.out
Do you know how I can resolve this and what the problem is? Also are the used subroutines and functions like second_derivative in the module private? Or if not how can I make sure that they are?
I currently try to run this .f90 file:
program LowPrandl
use module_LowPrandl
use iso_fortran_env
implicit none
type(grid_variables) :: all
real(real32) :: time
INTEGER :: step
call all%read_inputs()
call all%initialise()
time = 0.
step = 0
do while (time<all%total_time)
call all%VelocitySolve()
call all%Timestep()
time = time + all%dt
step = step + 1
print*, 'Time, step, maxUx, maxUy:', time, step, maxval(abs(all%ux)), maxval(abs(all%uy))
if((mod(step,100)==1) .or.( time>=all%total_time)) call all%write(step, time)
end do
end program LowPrandl
using this module:
module module_LowPrandl
use iso_fortran_env
implicit none
private
public :: grid_variables
real(real64) :: Pi = atan(1.)*4.
type :: grid_variables
!input parameters
real(real64) :: total_time = 0.1 !total integration time
character(len=10) :: Tinit='cosine' !Temperature initial condition
integer :: nx, ny !grid spacing
real(real64) :: a_dif, a_adv !time step constant
real(real64) :: err = 1e-6 !maximum error in Poisson solves = resiude_rms/rhs_rms<err
real(real64) :: Ra !Rayleigh number
real(real64) :: Pr !Prandt number
!character(len=100) :: fname = 'parameter.txt' !read command-line arguments
!Variables in type
real(real64), allocatable :: T(:,:),dTdx(:,:), d2Tdxdy(:,:),d2Wdxdy(:,:)
real(real64), allocatable :: ux(:,:), uy(:,:)
real(real64), allocatable :: S(:,:), ugradT(:,:), ugradW(:,:), W(:,:)
real(real64) :: dx, dy, h, dt = 0.1, dt_dif = 0.1, dt_adv = 0.1
contains
procedure,public,pass(this) :: read_inputs
procedure,public,pass(this) :: initialise
procedure,public,pass(this) :: VelocitySolve
procedure,public,pass(this) :: Timestep
procedure,public,pass(this) :: write
end type grid_variables
contains
!the maine subroutines to handle the translation of generic functions for the specific need -> should be public
subroutine read_inputs(this)
class(grid_variables), INTENT(inout) :: this
integer :: tmp_nx, tmp_ny
real(real64) :: tmp_a_dif, tmp_a_adv, tmp_err,tmp_Ra, tmp_Pr, tmp_total_time
character(len=100) :: tmp_fname
character(len=10) :: tmp_Tinit
namelist / inputs / tmp_nx,tmp_ny,tmp_total_time,tmp_Ra,tmp_err,tmp_a_dif,tmp_a_adv,tmp_Tinit,tmp_Pr
!read command-line arguments
if ( command_argument_count() > 0 ) then
call get_command_argument(1,tmp_fname)
end if
!inputs is the name of the list
OPEN(1,file=tmp_fname, status ='old')
READ(1, inputs)
CLOSE(1)
this%nx = tmp_nx
this%ny = tmp_ny
this%total_time = tmp_total_time
this%Ra = tmp_Ra
this%err = tmp_err
this%a_dif = tmp_a_dif
this%a_adv = tmp_a_adv
this%Tinit = tmp_Tinit
this%Pr = tmp_Pr
end subroutine read_inputs
subroutine initialise(this)
class(grid_variables), INTENT(inout) :: this
integer :: counter
allocate(this%T(this%nx,this%ny), this%d2Tdxdy(this%nx, this%ny), this%ugradT(this%nx, this%ny))
allocate(this%S(this%nx, this%ny), this%ux(this%nx,this%ny), this%uy(this%nx, this%ny),this% W(this%nx,this%ny))
allocate(this%ugradW(this%nx,this%ny), this%d2Wdxdy(this%nx,this%ny), this%dTdx(this%nx,this%ny))
!initalize stuff
this%S = 0.
this%W = 0.
this%ux = 0.
this%uy = 0.
this%ugradT = 0.
this%d2Tdxdy = 0.
this%d2Wdxdy = 0.
!spatial mesh distance
this%dy = 1./(this%ny-1)
this%dx = this%dy !h=dy
this%h = min(this%dx,this%dy)
!numerical computation goes here
if (this%Tinit=='random') then
call random_number(this%T)
call random_number(this%W)
this%T(1,:) = 0
this%T(this%nx,:) = 0
this%T(:,1) = 0
this%T(:,this%ny) = 0
else if (this%Tinit == 'spike') then
this%T = 0.
this%T(this%nx/2,this%ny/2) = 1
else if ( this%Tinit=='cosine' ) then
do counter = 1, this%nx
this%T(counter,:) = 0.5*(1+cos(3*Pi*(counter-1)/(this%nx-1)))
end do
else
print*, 'ERROR! Tinit must be either ^spike^, ^cosine^ or ^random^!'
stop
end if
!diffusive time step
this%dt_dif = abs(this%a_dif*min(this%dx, this%dy)**2/max(this%Pr,1.))
end subroutine initialise
subroutine VelocitySolve(this)
class(grid_variables), INTENT(inout) :: this
this%S = solve_Poisson_Eq_multigrid(this%S, this%W, This%h,this%err)
call S_to_vel(this%S, this%ux, this%uy, this%nx, this%ny, this%dx, this%dy)
end subroutine VelocitySolve
subroutine Timestep(this)
class(grid_variables), INTENT(INOUT) :: this
!calculate advective timestep and min(adv,diff)
!take advection and diffusion time step
if ( (maxval(abs(this%ux(:,:))) <= 0.0001).and.(maxval(abs(this%uy(:,:))) <= 0.0001) ) then !check if velocity = 0
this%dt = this%dt_dif
else
!Calculate dt_adv from vx, vy
this%dt_adv = this%a_adv*min(this%dx/(maxval(abs(this%ux(:,:)))), this%dy/(maxval(abs(this%uy(:,:)))))
this%dt = min(this%dt_dif, this%dt_adv)
end if
if ( this%dt > 0.6 * this%dx .OR. this%dt > 0.6 * this%dy) then
this%dt = min(this%dx, this%dy) * 0.6
print*, 'dt was corrected!'
end if
!set boundaries
this%T(1,:) = this%T(2,:)
this%T(this%nx,:) = this%T(this%nx-1,:)
this%T(:,1) = 1.
this%T(:,this%ny) = 0.
this%W(:,1) = 0.
this%W(:,this%ny) = 0.
this%W(1,:) = 0.
this%W(this%nx,:) = 0.
!Calculate del2(T)
this%d2Tdxdy = second_derivative(this%T,this%dx,this%dy)
!Calculate del2(W)
this%d2Wdxdy = second_derivative(this%W,this%dx,this%dy)
!Calculate u.grad(T)
this%ugradT = u_grad(this%ux, this%uy, this%T, this%dx, this%dy)
!Calculate u.grad(W)
this%ugradW = u_grad(this%ux, this%uy, this%W, this%dx, this%dy)
!Calculate dTdx
this%dTdx = first_derivative_centred(this%T,this%dx, 'x')
!update variables according to time step
this%T(2:this%nx-1,2:this%ny-1) = this%T(2:this%nx-1,2:this%ny-1) + this%dt * &
(this%d2Tdxdy(2:this%nx-1,2:this%ny-1) - this%ugradT(2:this%nx-1,2:this%ny-1))
this%W(2:this%nx-1,2:this%ny-1) = this%W(2:this%nx-1,2:this%ny-1) + this%dt * &
(this%Pr*this%d2Wdxdy(2:this%nx-1,2:this%ny-1)-this%ugradW(2:this%nx-1,2:this%ny-1) - &
this%Ra*this%Pr*this%dTdx(2:this%nx-1,2:this%ny-1))
end subroutine Timestep
subroutine write(time_step, time_iteration, this)
class(grid_variables), INTENT(inout) :: this
INTEGER, INTENT(IN) :: time_step
real(real32),INTENT(IN) :: time_iteration
character(len=60) :: outputfile_timestep_T = ''
character(len=60) :: outputfile_timestep_S = ''
character(len=60) :: outputfile_timestep_W = ''
character(len=60) :: outputfile_uxmax_vs_time = ''
character(len=60) :: outputfile_uymax_vs_time = ''
character(len=10) :: Pr_number = ''
!write Prandl number to file
write(Pr_number,'(F8.3)') this%Pr
call StripSpaces (Pr_number)
!write time step to file
write(Pr_number,'(F8.3)') time_step
call StripSpaces (Pr_number)
if ( time_step == 0 ) then
outputfile_uxmax_vs_time = trim(Pr_number)//'_uxmax_vs_time.data'
open(11,file=outputfile_uxmax_vs_time,action='write',position='append')
outputfile_uymax_vs_time = trim(Pr_number)//'_uymax_vs_time.data'
open(12,file=outputfile_uymax_vs_time,action='write',position='append')
end if
WRITE(11,*) maxval(this%ux)
WRITE(12,*) maxval(this%uy)
if((mod(time_step,100)==1) .or.( time_iteration>=this%total_time)) then
outputfile_timestep_T =trim(Pr_number)//'_T_timestep_' //trim(Pr_number)//'.bin'
outputfile_timestep_S =trim(Pr_number)//'_S_timestep_' //trim(Pr_number)//'.bin'
outputfile_timestep_W =trim(Pr_number)//'_W_timestep_' //trim(Pr_number)//'.bin'
print*, 'Calculating.... ', (time_iteration/this%total_time)*100, '%'
open(1,file=outputfile_timestep_T,form='unformatted',access='stream')
WRITE(1) this%nx, this%ny, this%T
close(1)
open(1,file=outputfile_timestep_S,form='unformatted',access='stream')
WRITE(1) this%nx, this%ny, this%S
close(1)
open(1,file=outputfile_timestep_W,form='unformatted',access='stream')
WRITE(1) this%nx, this%ny, this%W
close(1)
end if
if(time_iteration>=this%total_time) then
!write results for plot
WRITE(11,*) maxval(this%ux)
CLOSE(11)
WRITE(12,*) maxval(this%uy)
CLOSE(12)
end if
end subroutine write
!generic functions used -> shuld be private
function second_derivative(func,hx,hy)
real(real64), INTENT(IN) :: hx, hy
real(real64) :: func(:,:)
INTEGER :: np_x, np_y
real(real64), ALLOCATABLE :: func_second_deriv(:,:)
real(real64) :: second_derivative(size(func,1),size(func,2))
INTEGER :: i,j !local variable
ALLOCATE(func_second_deriv(size(func,1),size(func,2)))
np_x = size(func,1)
np_y = size(func,2)
func_second_deriv(1,:) = 0
func_second_deriv(np_x,:) = 0
func_second_deriv(:,1) = 0
func_second_deriv(:,np_y) = 0
do i = 2, np_x - 1
do j = 2, np_y - 1
func_second_deriv(i,j) =(func(i-1,j)+func(i+1,j) - 2*func(i,j))/hx**2 + &
(func(i,j-1) + func(i,j+1) - 2*func(i,j))/hy**2
end do
end do
!Assume derivatives at the end points
second_derivative = func_second_deriv
DEALLOCATE(func_second_deriv)
end function second_derivative
function first_derivative_centred(func,h,variable)
character(len=1) :: variable
real(real64), INTENT(IN) :: h
real(real64) :: func(:,:)
INTEGER :: np_x, np_y
real(real64), ALLOCATABLE :: func_first_deriv(:,:), first_derivative_centred(:,:)
INTEGER :: i,j !local variable
ALLOCATE(func_first_deriv(size(func,1),size(func,2)),first_derivative_centred(size(func,1),size(func,2)))
np_x = size(func,1)
np_y = size(func,2)
func_first_deriv(1,:) = 0
func_first_deriv(np_x,:) = 0
func_first_deriv(:,1) = 0
func_first_deriv(:,np_y) = 0
if ( variable == 'x' ) then
do i = 2, np_x - 1
do j = 2, np_y - 1
func_first_deriv(i,j) = (func(i+1,j) - func(i-1,j)) / (2*h)
end do
end do
else if ( variable == 'y' ) then
do i = 2, np_x - 1
do j = 2, np_y - 1
func_first_deriv(i,j) = (func(i,j+1) - func(i,j-1)) / (2*h)
end do
end do
else
print*,'Error!!! The variable is not correct!!'
stop
end if
!Assume derivatives at the end points
first_derivative_centred = func_first_deriv
DEALLOCATE(func_first_deriv)
end function first_derivative_centred
function u_grad(velocity_x, velocity_y, T_field, hx, hy)
implicit none
INTEGER :: i, j
real(real64), INTENT(IN) :: velocity_x(:,:), velocity_y(:,:), T_field(:,:), hx, hy
real(real64), ALLOCATABLE :: u_grad_cal_x(:,:), u_grad_cal_y(:,:), u_grad(:,:)
ALLOCATE(u_grad_cal_x(size(T_field,1), size(T_field,2)))
ALLOCATE(u_grad_cal_y(size(T_field,1), size(T_field,2)))
ALLOCATE(u_grad(size(T_field,1), size(T_field,2)))
if ( (maxval(abs(velocity_x(:,:))) <= 0.0001).and.(maxval(abs(velocity_y(:,:))) <= 0.0001) ) then
u_grad = 0.
else
!boundaries
u_grad_cal_x(:,size(T_field,2)) = 0
u_grad_cal_x(size(T_field,1),:) = 0
u_grad_cal_x(:,1) = 0
u_grad_cal_x(1,:) = 0
u_grad_cal_y(:,size(T_field,2)) = 0
u_grad_cal_y(size(T_field,1),:) = 0
u_grad_cal_y(:,1) = 0
u_grad_cal_y(1,:) = 0
do i = 2, size(T_field,1) - 1
do j = 2, size(T_field,2) - 1
!ux*grad(T)
if ( velocity_x(i,j) > 0 ) then
u_grad_cal_x(i,j) = velocity_x(i,j) * ( T_field(i,j) - T_field(i-1,j) ) / hx
else if ( velocity_x(i,j) < 0 ) then
u_grad_cal_x(i,j) = velocity_x(i,j) * ( T_field(i+1,j) - T_field(i,j) ) / hx
end if
!uy*grad(T)
if ( velocity_y(i,j) > 0 ) then
u_grad_cal_y(i,j) = velocity_y(i,j) * ( T_field(i,j) - T_field(i,j-1) ) / hy
else if ( velocity_y(i,j) < 0 ) then
u_grad_cal_y(i,j) = velocity_y(i,j) * ( T_field(i,j+1) - T_field(i,j) ) / hy
end if
end do
end do
u_grad = u_grad_cal_x + u_grad_cal_y
end if
DEALLOCATE(u_grad_cal_x,u_grad_cal_y)
end function u_grad
recursive function Vcycle_2DPoisson(u_f,rhs,h) result (resV)
!resV result to programm
!u_f is the fine grid approximate solution
!rhs = f, is the right hand side
!h grid spacing
implicit none
real(real64) :: resV
real(real64),intent(inout) :: u_f(:,:) ! arguments
real(real64),intent(in) :: rhs(:,:), h
integer :: nx,ny,nxc,nyc, i ! local variables
real(real64),allocatable :: res_c(:,:),corr_c(:,:),res_f(:,:),corr_f(:,:)
real(real64) :: alpha=0.7, res_rms
nx=size(u_f,1); ny=size(u_f,2) ! must be power of 2 plus 1
if( nx-1/=2*((nx-1)/2) .or. ny-1/=2*((ny-1)/2) ) stop 'ERROR:not a power of 2'
nxc=1+(nx-1)/2; nyc=1+(ny-1)/2 ! coarse grid size
if (min(nx,ny)>5) then ! not the coarsest level
allocate(res_f(nx,ny),corr_f(nx,ny), &
corr_c(nxc,nyc),res_c(nxc,nyc))
!---------- take 2 iterations on the fine grid--------------
res_rms = iteration_2DPoisson(u_f,rhs,h,alpha)
res_rms = iteration_2DPoisson(u_f,rhs,h,alpha)
!---------- restrict the residue to the coarse grid --------
call residue_2DPoisson(u_f,rhs,h,res_f)
call restrict(res_f,res_c)
!---------- solve for the coarse grid correction -----------
corr_c = 0.
res_rms = Vcycle_2DPoisson(corr_c,res_c,h*2) ! *RECURSIVE CALL*
!---- prolongate (interpolate) the correction to the fine grid
call prolongate(corr_c,corr_f)
!---------- correct the fine-grid solution -----------------
u_f = u_f - corr_f
!---------- two more smoothing iterations on the fine grid---
res_rms = iteration_2DPoisson(u_f,rhs,h,alpha)
res_rms = iteration_2DPoisson(u_f,rhs,h,alpha)
deallocate(res_f,corr_f,res_c,corr_c)
else
!----- coarsest level (ny=5): iterate to get 'exact' solution
do i = 1,100
res_rms = iteration_2DPoisson(u_f,rhs,h,alpha)
end do
end if
resV = res_rms ! returns the rms. residue
end function Vcycle_2DPoisson
function iteration_2DPoisson(u,f,h,a)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Gauss-Seidel iteration !
! is a function that returns the root mean square of the residual !
! u is the approximated solution !
! f is what is on the right hand side of the equal sine !
! h is the grid spacing !
! a is alpha which is the relaxation parameter !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer :: i, j!local variable
real(real64) :: h, a, f(:,:)
real(real64) :: u(:,:)
real(real64), allocatable :: residue (:,:)
real(real64) :: iteration_2DPoisson
allocate(residue(size(u,1),size(u,2)))
residue = 0.
do concurrent (i = 2:size(u,1)-1, j=2:size(u,2)-1)
residue(i,j) = (u(i,j+1) + u(i,j-1) + u(i+1,j) + u(i-1,j) - 4*u(i,j)) / h**2 - f(i,j)
u(i,j) = u(i,j) + a * (h**2)/4 * residue(i,j)
end do
iteration_2DPoisson = RMS(residue)
DEALLOCATE(residue)
end function iteration_2DPoisson
subroutine residue_2DPoisson(u,f,h,res)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! calculates the residue in array res !
! u is the approximated solution !
! f is what is on the right hand side of the equal sine !
! h is the grid spacing !
! res is residual array !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer :: i, j !local variable
real(real64) :: h, f(:,:)
real(real64) :: u(:,:)
real(real64), INTENT(OUT), DIMENSION(size(u,1),size(u,2)) :: res
!real(real64), INTENT(OUT),ALLOCATABLE :: res(:,:) !output residual
!ALLOCATE(res(size(u,1),size(u,2)))
res = 0.
do concurrent(i = 2:size(u,1)-1, j=2:size(u,2)-1)
res(i,j) = (u(i,j+1) + u(i,j-1) + u(i+1,j) + u(i-1,j) - 4*u(i,j)) / h**2 - f(i,j)
end do
end subroutine residue_2DPoisson
subroutine restrict(fine, coarse)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! copies every other point from fine into coarse !
! fine is the fine grid !
! coarse is the coarse grid !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer :: i_fine, j_fine !local variable
real(real64), INTENT(IN) :: fine(:,:) !input grid
!real(real64), INTENT(OUT),ALLOCATABLE :: coarse(:,:) !output grid
real(real64), INTENT(OUT) :: coarse(size(fine,1)/2+1,size(fine,2)/2+1) !output grid
!ALLOCATE(coarse(size(fine,1)/2+1,size(fine,2)/2+1))
do concurrent(i_fine=1:size(fine,1), j_fine=1:size(fine,2))
if ( (mod(i_fine,2) > 0).and.(mod(j_fine,2) > 0)) then
coarse(i_fine/2+1,j_fine/2+1) = fine(i_fine,j_fine)
end if
end do
end subroutine restrict
subroutine prolongate(coarse, fine)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! copies coarse into every other point in fine !
! does linear interpolation to fill the other points !
! fine is the fine grid !
! coarse is the coarse grid !
!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer :: i, j !local variable
real(real64), INTENT(OUT) :: fine(:,:) !input grid
real(real64), DIMENSION(size(fine,1)/2+1,size(fine,2)/2+1) :: coarse !output grid
!real(real64), ALLOCATABLE :: coarse(:,:) !output grid
!ALLOCATE(coarse(size(fine,1)/2+1,size(fine,2)/2+1))
do concurrent(i=1:size(fine,1), j=1:size(fine,2))
if ( (mod(i,2) > 0).and.(mod(j,2) > 0)) then
fine(i,j) = coarse(i/2+1,j/2+1)
else if ((mod(i,2) > 0).and.(mod(j,2) == 0)) then
fine(i,j) = 0.5*(coarse(i/2+1,j/2) + coarse(i/2+1,j/2+1))
else if ((mod(i,2) == 0).and.(mod(j,2) > 0)) then
fine(i,j) = 0.5*(coarse(i/2,j/2+1) + coarse(i/2+1,j/2+1))
end if
end do
do concurrent(i=1:size(fine,1), j=1:size(fine,2))
if ((mod(i,2) == 0).and.(mod(j,2) == 0)) then
fine(i,j) = 0.25*(fine(i-1,j)+fine(i+1,j)+fine(i,j-1)+fine(i,j+1))
end if
end do
!DEALLOCATE(coarse)
end subroutine prolongate
real(real64) function RMS(A)
implicit none
real(real64), INTENT(IN) :: A(:,:)
INTEGER :: i,j
real(real64) :: sum_squared!, RMS
sum_squared = 0.0
do concurrent (i = 1:size(A,1), j = 1:size(A,2))
sum_squared = sum_squared + A(i,j)**2
end do
RMS = sqrt(sum_squared/(size(A,1)*size(A,2)))
end function RMS
function solve_Poisson_Eq_multigrid(u,f,h, crit_given)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Solve Poisson Equation using multigrid !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer :: iteration, iteration_max
real(real64) :: h, crit_given, crit_current, f(:,:)
real(real64) :: u(:,:), residual_func
real(real64), allocatable :: solve_Poisson_Eq_multigrid(:,:)
allocate(solve_Poisson_Eq_multigrid(size(u,1),size(u,2)))
!set values back
residual_func = 0. !changed frim residual to residual_func
iteration_max = 100
iteration = 0
crit_current = 1.
if ( (RMS(f) <= 0.0001).and.(RMS(f) >= - 0.0001) ) then !check if rms(f)==0
solve_Poisson_Eq_multigrid = 0.
else
do while ((crit_current >= crit_given).and.(iteration < iteration_max) )
residual_func = Vcycle_2DPoisson(u,f,h)
iteration = iteration + 1
crit_current = residual_func/rms(f)
end do
solve_Poisson_Eq_multigrid = u
end if
end function solve_Poisson_Eq_multigrid
subroutine S_to_vel(S,ux, uy, nx, ny, dx, dy)
implicit none
real(real64), allocatable,intent(inout) :: S(:,:)
integer, INTENT(IN) :: nx, ny
real(real64),intent(in) :: dx, dy
real(real64), allocatable,intent(out) :: ux(:,:), uy(:,:)
S(1,:) = 0.
S(nx,:) = 0.
S(:,1) = 0.
S(:,ny) = 0.
!Calculate ux, uy from S
ux = first_derivative_centred(S(:,:),dy,'y')
uy = -1*first_derivative_centred(S(:,:),dx,'x')
end subroutine S_to_vel
subroutine StripSpaces(string)
character(len=*) :: string
integer :: stringLen
integer :: last, actual
stringLen = len (string)
last = 1
actual = 1
do while (actual < stringLen)
if ( string(last:last) == '.' ) then
string(last:last) = '_'
end if
if (string(last:last) == ' ') then
actual = actual + 1
string(last:last) = string(actual:actual)
string(actual:actual) = ' '
else
last = last + 1
if (actual < last) &
actual = last
endif
end do
end subroutine
end module module_LowPrandl
and this input file:
&inputs
Pr=0.01 ! <---- CHANGE ONLY THIS
nx=257 ny=65
total_time=0.1
Ra=1.e6
err=1.e-3
a_dif=0.15 a_adv=0.4
Tinit='cosine'
/