MHD turbulence

This page includes:

Past Solved Problems with Solutions (for reference):

2004 - May 9-15 | May 16-22 | May 23-29 | May 30-June 5 | June 6-12 | June 13-19 |
2003 - [2003]

and link to MHD turbulence | MHD related papers

[May 9-15]

1. infinitely long runtime using compiler intel-7.0 (time step about
   E-07):
   abandoned and use compiler intel-8.0 (time step about E-2)
2. crashed using compiler intel-8.0 and mpirun,
   died in fluidx backward (mhdflux) - floating point error:
   trapifc.c (which deals with floating point exceptions)
   -- turn the second line of body into a comment
3. array bound exceeded in minus_momentum, normalize_energy, addv:
   mod(floor(x)-1,ndv)+1, mod(ceiling(x)-1,ndv)+1
4. when ndv=1,
   turbulence injected power is zero;
   normalized injected power is infinity:
   It is correct though unphysical, should not use ndv=1.
5. array bound exceeded in sweep:
   should use nxg,nyg,nzg greater than a certain minimum
   depending on xbuf.
6. too much comments at run time: turn off -C compiler option
7. syntax error when attempting to print in comm_module.f:
   In .f file, must have at least 6 blanks before the command
8. When debugger refuses to tell me any information or crashes
   with the program, try switching (by source) to the corresponding
   compiler
   (eg. idb with intel-7.0).
9. In CURL, using mod vs. modulo vs. component explicit has an effect on
   shape of measured power spectrum (may be compiler dependent):
   Analysis:
     Some properties using intel-8.0 compiler
     - mod -- max(div) about 0.5, smooth profile
     - modulo -- max(div) about E-10, not smooth
     - explicit -- max(div) about E-02, smooth
       (see plots 5-9)
     However, mod and modulo are not supposed to yield different
     results as the arguments are positive integers.  This deviation
     may be due to compiler bugs dealing with combined use of
     cshift and mod/modulo.
   Solution:
     Use component explicit method instead of mod or modulo.
10. When testing CURL, looked at different power spectrums:
    Let's fix ps2v=k**6 * exp(-4.03332*(k/kpk)**2) until specified
    otherwise.
11. Checked power spectrum for ndv.ne.nxg
    - smooth, similar to ndv=nxg, max(div) about E-02
      (see plot 10)
12. failed to run on two cpu using intel-7.0 compiler
    (died in comm_bufferupdate or later):
    Analysis:
      try running mhdIIa
        - mhd.f90: added comments (version on BOB)
        - died after transpose2 in buffer initialization using
          intel-8.0 compiler
        - exited normally using intel-7.0 compiler
      (a) The isothermal version with turbulence driving seems to run
          fine using intel-7.0 compiler.
      (b) the time step looks reasonable as far as the test goes,
          about E-02.  Previous problem seems to have disappeared.
      (c) max(div) is about E-10 while it is E-02 using intel-8.0,
          the measured 1D power spectrum is less smooth.
          (see plot 11)
    Solution:
      use intel-7.0 compiler from now on
13. random seeds are now synchronized between different processes.
    (checked using 'diff [output files]' and plot comparison)

[May 16-22]

1.  momentum is not conserved during advection:
    changed momentum calculation (also apply to energy)
       -- changed summation bounds to exclude buffers
       -- sensible only when outside of sweep or after completion of
          forward/backward sweep due to summation bounds error
    changed comm_module.f
       -- to output nxmn,nymn,nzmn,npercpu
       -- changed domain decomposition to cut along y-direction
          first, to allow mod(ncpu,8)=2 while still disliking =4
	  CANCELLED! (because it does not affect advection)
       Note: the constraint that nxmn=nzmn seems to have been removed.
2.  done testing subroutine outputenergy.
3.  check that all three components of dv are nonzero and reasonable.
4.  momentum after PERTURB .neq. momentum before PERTURB + momentum(dv)
    (in ADDV) even when ncpu=1, ndv=nxg and advection is off:
    the deviation is due to rounding error
    tried - use real(kind=8) dummy:=u
               -- helped improve precision
	       -- from now on, keep dummy for debugging
	          but not for large runs
          - double precision mpcom for momentum calculationg
	       -- wrong since MPI_ALLREDUCE only takes single
	          precision real arguments
    Solution:
      from now on, keep dummy for debugging, but not for large runs.
      from now on, use only single precision reals (eg. en and mpcom)
        in MPI_ALLREDUCE.
5.  energy after PERTURB .neq. energy before PERTURB + energy(dv)
    (in ADDV) even when (ncpu=1,) ndv=nxg and advection is on:
    Analysis:
      - when ndv=nxg, ncpu=1, advection off
	energy in ADDV is correct
      - when ndv.neq.nxg, ncpu=2, advection off
	energy in ADDV is correct
      - when ndv=nxg, ncpu=1, advection on
	energy in ADDV is incorrect
      - change dummy to real(kind=16)
        does not affect energy or momentum calculation
    Tests and Solution:
      - test mhd2 for for momentum conservation -- satisfactory
      - test static inhomogeneous initial density profile
        with perturbation and without advection -- energy not additive
      - added printing dv**2 and dv*v terms in normalize_energy
        -- cross term equals zero (since originally 'en' real(kind=8))
	   FIXED this problem by making 'en' simply real
        -- energy is not normalized correctly
	   FIXED this problem by
	     - removing the extra rho in normalized_energy which
	       affects the cross term, enloc(2)
	     - also removing the extra rho in addv perspective energy
	-- removed the added print
    Post-testing:
      - when ndv=nxg, ncpu=1, advection on (use SOG, homog. density)
	energy in ADDV is now correct
      - recover ndv.neq.nxg
        energy is correct
        momentum is negligible
      - recover cpu=2
        energy is correct
        momentum is negligible
    Note: momentum is still not quite additive (about E-03 to E-05)
          -- perturbation-introduced error is slightly greater
	     than advection-introduced error
          -- may be okay since it is very tiny in percentage error

[May 23-29]

1.  simplified record output
2.  time step problems - for nxg=32, feels like it takes too long
   and perhaps too many time steps:
   - expect 351 perturbations
   - can reduce coarse grid resolution to save time
   - for (ndv=128, kpk=4),
     time step limited by perturbation time interval
     and not by advection
   - try (ndv=64, kpk=4) -- faster
   - try (ndv=16, kpk=0.5)
     -- small dt (about 1E-07) after iter=351 (see 2a)
     -- small dt (about 1E-07) after iter=127 (see 2b)
2a. time step greater than perturbation interval:
    - when determining dt, changed from
      dt=min(0.7*cfl, 0.5*(tf-time)) to
      dt=min(0.7*cfl, 0.5*(tperturb-time))
    - accompanying the above change is the uncrucial change from
      time = min(time + 2.*dt,tf) back to
      time = time + 2.*dt
2b. 'time' never quite gets to 'tperturb' while 'cfl' is large enough:
    - print iter/time/dt/tf and delt/tperturb to 9 decimal places
    - limited by single precision 'time'
    - perturb when time.ge.(tperturb-time*1E-07)

[May 30-June 5]

1.  coherence time driving
    - created new directory SOG1-coherence/
    - turned off advection
    - first used a fixed coherence time
    - for now, disregard seed synchronization
    - files modified:
        parameter_profile.f90
        mhd.f90
        perturb.f90 - added 2 new arguments
        normalize_energy.f90
        mixdv.f90 - NEW!!
        outputdv.f90 - NEW!!
1a. implementing coherence time driving:
    fluid output of x-component at (i,j,k)=(1,1,1) as a function of time
    seems to alternate between 3 different values.
    - strange correlation is indeed created by mixdv
    - coefficient computation correct
    - due to lack of normalization before mixdv, old dv dominates
    - by mistake, I 'call curl(dv)' instead of 'curl(newdv)',
      so the dv dominated by old dv gets into a cycle of period 3.
    Solution:
    - in perturb.f90, 'call curl(dv)' is changed to 'call curl(newdv)'
      -- dv is still dominated by old dv
      -- fluid output changes from multiple horizonal lines to
         a single diagonal line decreasing with time
      -- linear increment (as expected)
    - now normalize dvnew before mixdv
      -- now zero momentum and normalize dvnew before mixdv
      -- normalize again after since kinetic energy does not obey
         superposition principle
      -- now transition of dv between all new and all old is smooth
         depending tcoh:
         - dvnew u1 ranges from -25 to 15
         - dvold u1 ranges from 0 to -200 (diagonal)
1b. testing coherence time driving:
    - constant coefficient to obtain certain expected injected energy
      (ndv=n=2)
      - slight asymmetry of gaussian distribution - okay now
        (see ConstNorm2/randomSOG.dat)
      - a little more symmetric about zero, still slightly negative
        (compare ConstNorm2/u1_0.002.dat and Dynorm/u1_0.002.dat)
      - correlation (see ConstNorm2/timecor.eps)
      - energy (see ConstNorm2/energy_dr_002_002_00.eps)
      - power spectrum
        -- ConstNorm2/psV_002_002_00.eps
        -- currently half-body-diag-long 1D power spectrum
        -- not worth seeing for n=2
    - whole body diagonal powerspectrum - disgarded
      (ndv=n=2)
      - random number -- okay
      - u1 -- indistinguishable from ConstNorm2/u1_0.002.dat
      - correlation -- should be similar to ConstNorm2/timecor.eps
      - energy
        -- BodyDiag3a/energy_dr_002_002_00.eps
        -- indistinguishable from ConstNorm2/energy_dr_002_002_00.eps
      - power spectrum
        -- BodyDiag3a/psV_002_002_00.eps
        -- indistinguishable from ConstNorm2/psV_002_002_00.eps
        -- added print statements in convolve, check, perturb
           c1d_to_c3d and c3d_to_c1d
    - change back to half body diagonal
      - minor modification of c1d_to_c3d and c3d_to_c1d
        -- 'if(ix.ge.one*n/2) ix1=ix1-n' becomes
           'if(ix.gt.one*n/2) ix1=ix1-n' etc.
      - linear interpolation coefficients look correct
    - larger iterf=1000,10000 (ie. let it run for longer)
      -- tcoh=0.002
         - u1 -- largely negative (-120 to 20) ??
         - correlation -- not sharp enough
         - energy -- looks good (saturation about 1E5)
      -- tcoh=1
         - u1 -- even more largely negative (-180 to 0) ??
	 - correlation -- a little less sharp
 	 - energy -- looks good (saturation around 1E6)
    - now measure the time correlation of dv instead
      - added outputdv.f90
      - use iterf=10000
      - use kpk=0.25
        - tcoh=1: okay
        - tcoh=5: okay
        - tcoh=10: okay
        - dv(t) length not being a multiple of 2 doesn't seem to matter
      - use kpk=4.0 (better, include comparison with theoretical)
        - tcoh=5: good (see Plot 15)
        - tcoh=10: good (see Plot 16)
        - tcoh=100: good (see Plot 17)
        - dv(t) length is a multiple of 2
3.  dv power spectrum
    - use iterf=1
    - use kpk=0.25: ndv=32, 128 (see Plot 13)
      -- should not use small kpk, such as this, because
         actual power spectrum cannot achieve the theoretical sharpness
    - use kpk=4: ndv=128 (see Plot 14)
      -- looks quit good
4.  random_seed
    - necessary when ncpu.neq.1
    - re-implemented random_seed synchronization
      -- set seed(1:47)=0 and then put seed
      -- seeds happen to be synchronized on devel1 and devel2

[June 6-12]

1.  when n=ndv=32 (not when n=ndv=4,16), the cross section (z=1,nzg/2+1)
    of velocity field looks amplified near y=1 and y=nyg.
    Analysis:
      - independent of ncpu
      - independent of advection
      - not due to time correlation
      - due to random seed (incidentally, for ndv=32 only)
    Solution:
      - set seed=(1,..,47) instead of (seed=0 or seed=1)
      - looks good w/o advection, w/o time correlation, and iterf=5
      - try iterf=501, mixdv with tcoh=0.01
        -- note: strong buildup of velocity controlled by edot (okay)
      - try adding advection
        -- note: buildup of v-field is supressed by advection

[June 13-19]

1.  simulation tuning: observe the dependence of avg(de)
    - while setting rdvnew=1, and iterf large enough
    - on ndv, nxg, tcoh, kpk (and seed)
    - no dependence on ncpu
      eg. ncpu=1 and ncpu=8 gives the same avg(de) for nxg=ndv=32
2.  simulation for analysis: output for fluid power spectrum computation
    - added new subroutine output3d(34)
    - called once at the end
3.  post-analysis: preliminary (small) fluid power spectrum
    - no dependence on ncpu (eg. n=32, ncpu=1 vs. ncpu=8)
    - created programs:
      -- source files in pet_fat/Projects/MHD/Isoth/PowerSpectrum/
      -- data in mammoth:/cita/d/scratch-month/pet_fat/PowerSpectrum/
         dumped in doug:/mnt/raid/hoser1/pet_fat/
	 (takes an hour to tar, 2 minutes to move the tar file,
	 and a few minutes to extract)
    - analyzed:
      -- ncpu=512, nxg=32, iterf=100
4.  boundary condition
    - create new directory SOG3-boundary/ (future MPI version)
      - added new subroutine boundary_b (35) !! - need MPI fftw
      - the only file that's somewhat worth keeping: boundary_b(35)
    - created serial codes to test outer-BC solving scheme:
      - made new programs in ~/Projects/MHD/Isoth/Boundary/
5.  u power spectrum on the fly (possible MPI version)
    - can be replaced by post-analysis on eg. octopus
    - installed fftw
    - added new subroutine ps_ub (36) !! - need MPI fftw
      - still need to use, code (plan), compile
    - may need to modify c3d_to_c1d

[July 4-10]

1.  put in readfile
    - for time limit and decaying turbulence

[2003]

Tests:
0. adiabatic - mhd.f90 problems caused by non-constant cs
1. test1 (MPI) - graph looks exactly like previous 'minmod' under 'smooth minmod'
2. test2 (MPI) - looks unsmooth, has kinks
3. test2 (MPI) - error calculation
4. test2 (MPI) - rho sm: analytic curve first half looks especially ugly

SOG:
0. when one uses transposeb
1. vperturb (MPI) - static array more efficient?
2. SOG - implement subroutines added in mhd.f90 and makefile

Now fixed:
1. makefile - forget to add new *.f90 files, effects?
2. mhd.f90 - call trapfpe()
3. mhd.f90 - 0.8*cfl should be 0.7*cfl!
4. adiabatic - energy as initial condition?