CITA |
Gravitational Lensing |
NASA Images |
Search Archieves | My Research | My Photos | Contact Info | Home |
Liuxuan's HomepageLast Update: 30 June 2004 |

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

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)

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 changedcomm_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

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)

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

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

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

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

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?

CITA |
Gravitational Lensing |
NASA Images |
Search Archieves | My Research | My Photos | Contact Info | Home |
Liuxuan's Homepage© 2002 by Lucy Liuxuan ZhangLast Update: 30 May 2004 |