Program driver ! -------------------------------------------------------------------------- ! MPI Version of the Cartesian Acoustic Sun Simulator. ! Copyright 2006, Shravan Hanasoge ! Hansen Experimental Physics Laboratory ! 455 Via Palou way, Stanford ! CA 94305, USA ! Email: shravan@stanford.edu ! The driver routine. It utilizes the other subroutines to perform time- ! advancements. The time advancement is performed with the explicit Low ! Dissipation, Low Dispersion Runge Kutta 4th order method, described in ! Hu et al. (1996). Radial derivatives are computed using compact finite ! difference (Lele 1992). Horizontal variations are extracted by projecting ! the desired quantity onto Fourier harmonic space followed by an ! appropriate inverse fourier transform. ! ! BASIC CHECKS: ! 1. If init is all right. ! 2. If the time-length of the simulation matches ! 3. If the directories are OK ! 4. If the forcing function is being read from the right directory ! 5. If the variables are initialized OK ! 6. If the flow profiles are allright ! 7. If output and status are deleted or appropriately moved ! ! --------------------------------------------------------------------------- use initialize use all_modules use mp_physics use time_step implicit none integer i,j, init, ierr, t, k, index(1000,2), randinit(2), aind real*8 start_time,end_time,tempf,t1,T00,nor1,nor2,nor3,nor4,nor5,e, mp_time real*8 data(nz,6),total_time, start_mp_time, avg_time,zz, tempxy, rand1, rand2 logical saved, first_iteration character*5 tempc real*8, allocatable, dimension(:,:,:) :: tempo,check,outpu real*8, allocatable, dimension(:,:,:,:) :: bg 20 format(8f12.4) ! Initializing the computation ! -------------------------------------------- call Initialize_all() start_mp_time = MPI_WTIME() start_time = MPI_WTIME() ! Read in rho,P,c,gamma,g ! Determine nondim factors ! calculate c^2 call read_in_background() if (rank==0) then print *, 'dimc = ',dimc print *, 'diml = ',diml print *, 'dimrho = ',dimrho endif if (flows) then call read_in_flows() endif if (magnetic) then call read_in_field() endif ! STABILIZATION ALGORITHM FOR CREATING CONVECTIVELY STABILIZED MODELS (CSM's) ! from Hamed Moradi, 1-10-10 ! See Schunker et al., 2010 for details if (CSM_A) then if (rank==0) print *, 'RUNNING SOUND SPEED MODIFICATION FOR CSM A' do k=1,nz c2(:,:,k) = c2(:,:,k) *(1. + 0.1 * exp(-((Rsun-Rsun*z(k))*10.0**(-8.) )**2. ) ) c_speed(:,:,k) = c2(:,:,k)**0.5 c2rho0(:,:,k) = c2(:,:,k)*rho0(:,:,k) cq(k) = c_speed(1,1,k) enddo endif if (CSM_B) then if (rank==0) print *, 'RUNNING SOUND SPEED MODIFICATION FOR CSM B' do k=1,nz c2a(:,:,k) = c2(:,:,k) *(1. + 0.1 * exp(-((Rsun-Rsun*z(k))*10.0**(-8.) )**2. ) ) c2(:,:,k) = c2a(:,:,k)* (1. - 0.03 * exp(-((Rsun-Rsun*z(k))/(5*10.0**(8.)))**2.)) c_speed(:,:,k) = c2(:,:,k)**0.5 c2rho0(:,:,k) = c2(:,:,k)*rho0(:,:,k) cq(k) = c_speed(1,1,k) enddo endif ! END SOUND SPEED MODIFICATION FOR STABILIZATION ALGORITHM ! FILTER BACKGROUND MODEL ! Note: currently no option for filtering ! background flows if (FILTERING_BG) then allocate(bg(nx,dim2(rank),nz,dimfive)) bg(:,:,:,1) = rho0 bg(:,:,:,2) = gamma bg(:,:,:,3) = c_speed bg(:,:,:,4) = g bg(:,:,:,5) = p0 if (magnetic) then bg(:,:,:,6) = box bg(:,:,:,7) = boy bg(:,:,:,8) = boz endif if (TEST_IN_2D) then print *, '2D FILTERING BACKGROUND MODEL' bg = filter_vars_2d(bg) elseif (.not. TEST_IN_2D) then print *, '3D FILTERING BACKGROUND MODEL' bg = filter_vars(bg) endif rho0 = bg(:,:,:,1) gamma = bg(:,:,:,2) c_speed = bg(:,:,:,3) g = bg(:,:,:,4) p0 = bg(:,:,:,5) if (magnetic) then box = bg(:,:,:,6) boy = bg(:,:,:,7) boz = bg(:,:,:,8) endif c2 = c_speed**2. c2rho0 = c2*rho0 cq = c_speed(1,1,:) endif !!!!!!! ADDITIONAL SOUND-SPEED AND FLOW PERTURBATIONS HERE !!!!!! ! Birch et al. 2009 sound speed perturbation algorithm was located here ! !!!!!!!!! READ IN INITIAL CONDITIONS !!!!!!!!!!!! ! set init_variables = .true. in params.i if (init_variables) then !print *, 'Initial Conditions Not Currently Implemented... Stopping!' !stop call readfits(dirbackmag//'P_init.fits',a(:,:,:,5),nz) a(:,:,:,5) = a(:,:,:,5)/(dimrho*dimc**2) endif if (.not. TEST_IN_2D) call mp_initialize_RHS() if (TEST_IN_2D) call mp_initialize_RHS_2d () end_time = MPI_WTIME() if (rank == 0) & print *, 'Initialization took: ', (end_time -start_time) !--------------------------------------------- ! LOAD INITIAL CONDITIONS FOR RESTART STATE ! vars are nondimensionalized after loading if (time0 .ne.0) call read_in_initial_condition ! Non-dimensional form, for the calculations deltat = timestep*dimc/diml if (rank == 0) & print *,'With cadence of 60 s, the timestep is', deltat*diml/dimc ! Length of simulation maxtime = floor(solartime*3600.0/timestep) + 1 maxtime = 2*(maxtime/2) total_time = wall_time * 3600. steps = floor(outputcad/timestep + 0.1) ! HORIZONTAL FREQUENCY FILTERING ! freq_filtering = floor(60./timestep + 0.1) ! FILTERING HORIZONTALLY EVERY TIMESTEP ! freq_filtering = 1 if (rank ==0) then print *,'FREQUENCY OF OUTPUT AT', steps, 'TIMESTEPS' print *, 'NUMBER OF TIMESTEPS =', maxtime print *, 'TIMESTEP =', timestep endif ! Excitation series computation ----------------------------- num_steps = FLOOR(DBLE(maxtime)*timestep/cadforcing) + 2 num_steps = 2*(num_steps/2)+2 saved = .false. allocate(vr(nx,dim2(rank),num_steps)) vr(:,:,:) = 0 ! set forcing to zero if (toggle_force) then if (.not. TEST_IN_2D) then if (rank==0) print *, 'NUMBER OF FORCING STEPS = ', num_steps call readfits(forcingfunc,vr,num_steps) elseif(TEST_IN_2D) then allocate(tempin(nx,nx,num_steps)) tempin = 0.0 call readfits_1proc(forcingfunc,tempin,nx,nx,num_steps) vr(:,1,:) = tempin(:,nx/2,:) deallocate(tempin) endif ! SOURCE APODIZATION (TEMPORAL AND SPATIAL) ! Note if periodic horizontal boundaries are not, there will be no ! horizontal spatial apodization do k=1,num_steps do j=1,dim2(rank) do i=1,nx !tempxy=(((x(i)-0.5)*200)**2.0+((y(j+ystart(rank)-1)-0.5)*200)**2.0)**0.5 vr(i,j,k) = vr(i,j,k) * (1.-spongex(i)) * (1. - spongey(j)) /(1.0+ exp((5.-k))) !* 1.0/(1.0+exp((10.0-tempxy)*1.5)) enddo enddo enddo endif time_old = 0 !!!!! SOURCE STRENGTH PERTURBATIONS HERE (i.e., in the above loop)!!!!!! !----------------Driver Kernel---------------------------- ! RMS_HIST STORES THE RMS TIME EVOLUTION OF v_z ! A VERY USEFUL INDICATOR OF STABILITY if (rank == 0) & open(19,file='rms_hist',status='replace',form='formatted',position='append') if (time0 == 0) init = 0 if (time0 .NE. 0) init = time0+1 call Initialize_step() do i = init,maxtime time = i ! code to be timed call cpu_time(t1) T00 = t1 start_time = MPI_WTIME() call lagrange_interp () call step() end_time = MPI_WTIME() call cpu_time(t1) nor1 = norm(a(:,:,:,4)) !* dimc mp_time = MPI_WTIME() ! Timing and 'convergence' statistics ! The criterion used to see if the results are 'convergent' is the ! L2 norm of the radial velocity. Quite silly, really. if (rank ==0) then print *, '----------------------------------------------------------' print *, ' Iteration #, Cpu Time and Real Time:' print *, time, (t1 - T00), (time+1)*timestep print * print *,' Wall clock time', (end_time - start_time) print * print *,'L2 norm of vz = ',nor1 print *, '----------------------------------------------------------' ! DOMBROSKI ! write(19,*) nor1 if (nor1 .GT. 10.0**(3)) then ! Defualt threshold: 10.0**(-3.0) write(*,*) "Warning: nor1 > threshold; stopping..." stop endif !!!!!!!!!!! endif if (mod(time,steps) == 0) then if (rank == 0) print *, 'Writing output, time = ', time call convert_to_string(time,tempc,5) ! Magnetogram and thermo_vars were subroutines in old code !if (WRITE_MAG) call magnetogram() !if (WRITE_THERMO) call thermo_vars() ! THIS IS NEW OUTPUT CODE FOR 2D IMPLEMENTATION: ! files are ouput in cgs units call write_out_full_state !if (.not. test_in_2D) call write_out_slice(o_Rad) !if (test_in_2D) call write_out_full_state if (rank ==0 ) print *,'Full output done, index:',time if ((mp_time - start_mp_time) > (total_time - 150.0)) then write(*,*) "Stopping at L270 driver.f90" stop endif endif ! SHRAVAN'S STUFF: !if (( ((mp_time - start_mp_time) > (total_time - 150.0)) .OR. (mod(time+1,500) == 0) )) then ! ! call write_out_full_state ! ! if (rank ==0 ) print *,'Full output done, index:',time ! if ((mp_time - start_mp_time) > (total_time - 150.0)) stop ! !endif end do ! END OF TIME ITERATION call MPI_BARRIER(MPI_COMM_WORLD, ierr) call dfftw_destroy_plan(fftw_plan_fwd_x) call dfftw_destroy_plan(fftw_plan_inv_x) call dfftw_destroy_plan(fftw_plan_fwd_y) call dfftw_destroy_plan(fftw_plan_inv_y) if (rank == 0) then close(19) endif call MPI_BARRIER(MPI_COMM_WORLD, ierr) call MPI_FINALIZE(ierr) end program driver