forums.silverfrost.com

[yuvadon]

c
do k=1,NDIM
do i=1,N_BC_GIVEN(k)
V(k,I_BC_GIVEN(k,i)) = V_BC_GIVEN(k,i) ! ‹�ŠE��Œ�‚�“K—p
enddo
enddo
c
c==== output
c### call fileoutput(NUMP,X)
call fileoutput2(NUMP,X,DMG,T,V)
call fileoutput3(NUMP,X,SIGM)
c### stop
c
c==== initial neighboring table
call set_neigh
H = H0 ! ƒJ�[ƒlƒ‹ƒTƒCƒY�i‰e‹�”�Œa�j
c
c==== calculate relative distance between particles
call update_neigh
call rijcal(1)
c
c
c==== setting
DT = 1.0D-8 ! ŽžŠԍ�‚�•�ƒ�t�i�‰Š�’l�j
TIMEDIV = 1.0D0
c
c output_dt = TOTAL_TIME / 1000 ! ƒtƒ@ƒCƒ‹�o—�‚�‚�ŽžŠԍ�‚�
output_dt = TOTAL_TIME / 100 ! ƒtƒ@ƒCƒ‹�o—�‚�‚�ŽžŠԍ�‚�
noutput = 0 ! ƒtƒ@ƒCƒ‹�o—�‚�Š�Šu‚�—˜—p
c
c
c#### time increment ####
time_sim = 0.0D0
it = 0
c
do while(time_sim .lt. TOTAL_TIME)
it = it + 1
time_sim = time_sim + DT
c
c==== out for monitor
write(2,'(A,I7,4E13.5)') 'step',it,time_sim,DT,H,enrg
if(mod(it,100) .eq. 1) write(*,'(A,I7,3(A,E12.5))')
& 'step',it,' time',time_sim,' dt',DT,' H',H
c
do k=1,NDIM
do i=1,N_BC_GIVEN(k)
V(k,I_BC_GIVEN(k,i)) = V_BC_GIVEN(k,i) ! ‹�ŠE��Œ�‚�“K—p
enddo
enddo
c
c==== update neighboring list at every 10 steps
call update_neigh
c
c==== calculate artificial viscosity
call viscal
c
c
c<<<< Runge-Kutta method >>>>
do irk=1,4
c==== calculate pressure
call calc_press1 ! Mie-Gruneisen�F–�“x•�‰��C‰�“xŒ�‰�‘�
c### call calc_press2 ! Murnaghan�F–�“x•�‰��C‰�“xŒ�‰ʏ�
c
c==== calculate gradient of kernel function
call rijcal(0)
call calc_gradf
c
c==== calculate stress
call calc_stress(irk)
c
c==== calculate time-derivative of variables
call timederiv
c
c==== update quantities by Runge-Kutta integration
call rugkut_gill(irk)
c
do k=1,NDIM
do i=1,N_BC_GIVEN(k)
V(k,I_BC_GIVEN(k,i)) = V_BC_GIVEN(k,i) ! ‹�ŠE��Œ�‚�“K—p
enddo
enddo
c
c==== smoothing
call smoothing_CSA ! conservative smoothing
c### call smoothing_XSPH ! smoothing in XSPH
enddo
c<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>
c
c==== calculate relative distance and judge damage between particles
write(4,*) 'simulation time,',time_sim
call rijcal(2)
call judge_damage
c
c==== calculate yield stress by Johnson-Cook model
call calc_yield_JC(DT)
c
c==== update size of the kernel function
call kernel_size
c
c==== update time step
PDT = DT
call update_dt
c
c==== output
nout = int(time_sim / output_dt)
if(nout .ne. noutput) then
noutput = nout
c### call fileoutput(NUMP,X)
call fileoutput2(NUMP,X,DMG,T,V)
call fileoutput3(NUMP,X,SIGM)
endif
c
c==== store variables
enrg = 0.0D0
do ip=1,NUMP
HENSASOLD(1:3,ip) = HENSAS(1:3,ip)
enrg = 0.5D0 * PMASS(ip) * dot_product(V(:,ip),V(:,ip))
& + E(ip) * PMASS(ip)
enddo
c
c
c###############
c### if(mod(it,10) .eq. 1) call fileoutput2(NUMP,X,DMG,T,V)
c### if(it .eq. 200) stop
if(it .eq. MAXLOOP) stop
c###############
enddo ! time increment