C#######################################################################C C C C Method of lines for solving the graph variant of the elastic C C flow in the form: C C C C alpha v_Gamma_Phi = Delta_Gamma_Phi ( Kappa_Gamma_Phi + F ) C C C C Version 0.2.2.beta C C Playing with the i.c. shape C C-----------------------------------------------------------------------C C Comes from the isotropic versions of: C C Novebre 6, 2001, April 18, 2002, April 26, 2002, October 28, 2002 C C May 29, 2003, August 26, 2005, September 19-22, 2005 C C Date: June 6, 2006, September 4, 2006 C C C C#######################################################################C C C C Dr.Ing. Michal Benes C C Dept of Maths, FJFI CVUT C C (C) 2001 - 2006 C C C C#######################################################################C C f90 -K -O3 +U77 -o elast22_f_hp elast22_f.f C C ifc -pc80 -FI -o elast22_f_i elast22_f.f C C ifc -pc80 -FI -w95 -Vaxlib -o elast22_f_i elast22_f.f C C ifc -prec_div -FI -w95 -Vaxlib -o elast22_f_i elast22_f.f C C ifc -mp -FI -w95 -tpp7 -save -Vaxlib -o elast22_f_i elast22_f.f C C f90 -K -O3 +Oall +Ovectorize +DD64 +U77 -o C C elast22_f_ani_hp elast22_f_ani.f C PROGRAM Anisotropic_Elastic_by_Lines C time measurement CHARACTER*8 time_string INTEGER*2 hour, minute, second, hundrd C date and time CHARACTER*24 FDATE REAL*4 ETIME, ELAPSED, SYS(2) INTEGER date(8) CHARACTER (LEN=10) BEN(3) C simulation settings DOUBLE PRECISION T_STEP(500000), check_sum INTEGER*4 N_TST, TST_N DOUBLE PRECISION H1, H2, DT, F, ZETA, THETA_0, R, EPS, L1, L2, / Phi(-1:401,-1:401), DELTA, STEP, ALPHA, BLOW, / SCALE, a, b, c INTEGER FRQ, FIN_ID, WRITE_ID, CUR_ID, FLD_ID, CHK_ID, / SAVE_ID, LOAD_ID, SYS_ID LOGICAL FIRST COMMON /NUMER/ H1, H2, DT, N1, N2, EPS /MERSN/ STEP /CHECK/ BLOW / /PHYS/ ALPHA, F, ZETA, THETA_0, M_FOLD, R / /PATTERN/ SCALE /ID/ WRITE_ID /STORE/ SAVE_ID, LOAD_ID / /SYS_DATA/ SYS_ID, N_TST, T_STEP, TST_N, check_sum / /ANIS/ a, b, c C-----------------------------------------------------------------------C WRITE_ID = 1 WRITE (*,*) 'Method-of-lines scheme for solving' WRITE (*,*) 'the graph formulation of the elastic flow' WRITE (*,*) 'Version 0.2.2.beta' WRITE (*,*) 'Reading data' OPEN (96,FILE='elast22.dta',STATUS='OLD') READ (96,'(////////////2x,1F12.5,I12,5F12.5)') / F,M_FOLD,ZETA,THETA_0,R,SCALE READ (96,'(///2x,6F12.5)') L1, L2, ALPHA, a, b, c READ (96,'(//////2x,2I12,F12.5,I12,F12.5)') N1,N2,DT,NT,EPS READ (96,'(///2x,3F12.5)') DELTA, STEP, BLOW READ (96,'(//////2x,6I12///2X,3I12)') FRQ, FIN_ID, INI_ID, CUR_ID, / FLD_ID, CHK_ID, SAVE_ID, LOAD_ID, SYS_ID CLOSE(96) WRITE (*,*) 'Data loaded.' IF ((N1.GT.800).OR.(N2.GT.800)) THEN WRITE (*,*) 'Grid size exceeded!' STOP ENDIF H1 = L1/N1 H2 = L2/N2 PI = 3.14159D00 if ((STEP.gt.DT).or.(STEP.le.0.0D00)) STEP = DT C-----------------------------------------------------------------------C C timing block C system data initialization C N_TST = 0 TST_N = 0 check_sum = 0.D00 OPEN (96,FILE='elast22.dta',STATUS='OLD',ACCESS='APPEND') C PC call GetTim (hour, minute, second, hundrd) C Lin CALL DATE_AND_TIME (BEN(1),BEN(2),BEN(3),DATE) WRITE (96,*) 'Process started on:' C Lin ,DATE(3),'.',DATE(2),'.',DATE(1) ELAPSED = ETIME(SYS) WRITE (96,*)'Initial time: ' C Lin,DATE(5),':',DATE(6),':',DATE(7),'.', C Lin / DATE(8) WRITE (96,*) FDATE() write (96,*) 'Elapsed time = ', ELAPSED, ' User time = ', / SYS(1), ' System time = ', SYS(2) CLOSE (96) C initial condition WRITE (*,*) 'Initialization in progress...' IF (LOAD_ID.NE.1) THEN CALL INITIAL_LEVEL (Phi) ELSE CALL STORAGE (Phi) ENDIF K = 0 IF (INI_ID.EQ.1) THEN C OPEN (77,FILE='field.ini',STATUS='unknown') CALL SAVE_FIELD (Phi, K) C CLOSE (77) ENDIF WRITE (*,*) 'Initial level set' IF (CUR_ID.EQ.1) CALL EXTRACT_LEVEL_SET (Phi, 0) WRITE (*,*) 'Initialization terminated.' C----------------------------------------------------------------------C C time loop FIRST = .TRUE. DO K=1,NT CALL UPDATE_TIME_LEVEL((K-1)*DT,Phi,K*DT,N1,N2,DELTA,FIRST) IF (CHK_ID.EQ.1) CALL CHECK_BLOWUP(Phi, K) IF (MOD(K,FRQ).EQ.0) THEN WRITE (*,*) 'Time level No.: ', K IF (CUR_ID.EQ.1) CALL EXTRACT_LEVEL_SET (Phi, K) IF (FLD_ID.EQ.1) then CALL SAVE_FIELD (Phi, K) ENDIF ENDIF IF (SYS_ID.EQ.1) THEN OPEN(65,FILE='elast22.tst',STATUS='UNKNOWN',ACCESS='APPEND') WRITE (65,'(19H# Time level number,I10)') K DO I=1,N_TST check_sum = check_sum + T_STEP(I) WRITE(65,'(1X,I6,2X,F20.10)') TST_N+I, T_STEP(I) ENDDO CLOSE (65) WRITE (*,*) 'Partial number of time steps: ', N_TST TST_N = TST_N + N_TST N_TST = 0 ENDIF ENDDO C-----------------------------------------------------------------------C IF (FIN_ID.EQ.1) THEN WRITE (*,*) 'Saving the final form of the field Phi' C OPEN (77,FILE='field.fin',STATUS='unknown') CALL SAVE_FIELD (Phi, NT) C CLOSE (77) ENDIF IF (SAVE_ID.EQ.1) THEN CALL STORAGE (Phi) ENDIF C-----------------------------------------------------------------------C OPEN (96,FILE='elast22.dta',STATUS='OLD',ACCESS='APPEND') C PC call GetTim (hour, minute, second, hundrd) C Lin CALL DATE_AND_TIME (BEN(1),BEN(2),BEN(3),DATE) C PC WRITE (96,*)'Final time: ',hour,':',minute,':',second,'.',hundrd WRITE (96,*) 'Process terminated on:' C Lin ,DATE(3),'.',DATE(2),'.', C Lin / DATE(1) ELAPSED = ETIME(SYS) WRITE (96,*) FDATE() WRITE (96,*) 'Elapsed time = ', ELAPSED, ' User time = ', / SYS(1), ' System time = ', SYS(2) C Lin WRITE (96,*)'Final time: ' C Lin,DATE(5),':',DATE(6),':',DATE(7),'.', C Lin / DATE(8) if (SYS_ID.EQ.1) then WRITE (96,*) 'Number of time steps: ', TST_N WRITE (96,*) 'Number spatial nodes: ', / (N1+1)*(N2+1) WRITE (96,'(37H Total number of degrees of freedom: ,F25.0)') / DBLE(N1+1)*DBLE(N2+1)*DBLE(TST_N) Write (96,*) 'Control sum of time steps: ', check_sum endif CLOSE (96) WRITE (*,*) 'Program is successfully going to terminate.' WRITE (*,*) 'Time-step size at the end of computation: ', STEP STOP END C-----------------------------------------------------------------------C SUBROUTINE UPDATE_TIME_LEVEL (TB, Y, TE, N1, N2, DELTA, FIRST) C Runge-Kutta fourth-order scheme with automatic time stepping C DOUBLE PRECISION Y(-1:401,-1:401),Y0(0:400,0:400),K1(0:400,0:400), /K3(0:400,0:400),K4(0:400,0:400),K5(0:400,0:400),EPS,DELTA,Q DOUBLE PRECISION STEP,T,TB,TE,T0,H,H03, T_STEP(500000) INTEGER SYS_ID, N_TST LOGICAL LAST, FIRST COMMON /MERSN/ STEP COMMON /STACK/ Y0,K1,K3,K4,K5 COMMON /SYS_DATA/ SYS_ID, N_TST, T_STEP T = TB IF (FIRST) GO TO 101 H = STEP 102 IF (DABS(T-TE).LE.DABS(H)) GO TO 101 LAST = .FALSE. GO TO 103 101 LAST = .TRUE. H = TE - T 103 T0 = T DO 104 I=0,N1 DO 104 J=0,N2 104 Y0(I,J) = Y(I,J) H03 = H/3.D00 C write(*,*) 'Writing position: T=', T, ' H=', H CALL RIGHT_HAND_SIDE(T, Y, K1) DO 105 I=0,N1 DO 105 J=0,N2 105 Y(I,J) = K1(I,J)*H03 + Y0(I,J) T = T0 + H03 CALL RIGHT_HAND_SIDE(T, Y, K3) DO 106 I=0,N1 DO 106 J=0,N2 106 Y(I,J) = (K1(I,J) + K3(I,J))/2.D00*H03 + Y0(I,J) CALL RIGHT_HAND_SIDE(T, Y, K3) DO 107 I=0,N1 DO 107 J=0,N2 107 Y(I,J) = (K1(I,J)*0.375D00+K3(I,J)*1.125D00)*H03+Y0(I,J) T = T0 + H/2.D00 CALL RIGHT_HAND_SIDE(T, Y, K4) DO 108 I=0,N1 DO 108 J=0,N2 108 Y(I,J) = (K1(I,J)*1.5D00-K3(I,J)*4.5D00+K4(I,J)*6.D00)*H03+Y0(I,J) T = T0 + H CALL RIGHT_HAND_SIDE(T, Y, K5) DO 109 I=0,N1 DO 109 J=0,N2 109 Y(I,J) = ((K1(I,J)+K5(I,J))/2.D00+K4(I,J)*2.D00)*H03+Y0(I,J) EPS = 0.D00 C WRITE (*,*) 'EPS=0 ', EPS DO 110 I=0,N1 DO 110 J=0,N2 Q = DABS(K1(I,J)*0.2D00-K3(I,J)*0.9D00+K4(I,J)*0.8D00- / K5(I,J)*0.1D00) 110 IF (Q.GT.EPS) EPS = Q C write (*,*) 'Eps = ', EPS, ' Q= ', Q C write (*,*) 'Writing EPS= ', EPS, ' at T,H = ', T, H EPS = EPS*DABS(H03) C write (*,'(1x,E20.12)') Q, EPS IF (EPS.LE.DELTA) GO TO 111 DO 112 I=0,N1 DO 112 J=0,N2 112 Y(I,J) = Y0(I,J) T = T0 GO TO 113 C here we insert the time step measurement 111 if (SYS_ID.EQ.1) then N_TST = N_TST + 1 T_STEP(N_TST) = H endif IF (.NOT.LAST) GO TO 114 IF (FIRST) STEP = H RETURN 114 IF (EPS.EQ.0.0) GO TO 101 113 H = (DELTA/EPS)**0.2D00*0.8D00*H C write(*,*) 'H changed to ', H STEP = H FIRST = .FALSE. GO TO 102 END C-----------------------------------------------------------------------C SUBROUTINE RIGHT_HAND_SIDE (T, Phi, GG) DOUBLE PRECISION Phi(-1:401,-1:401), GG(0:400,0:400), T, / Kn(-1:401,-1:401), DT, F, H1, H2, ZETA, THETA_0, E, / R, EPS, ALPHA, vv1, vv2, vv3, Psi, T0_1, T0_2, T0_3, a, b, c DOUBLE PRECISION QQ(0:401,0:401), v1(0:401,0:401),v2(0:401,0:401), / NN1(0:401,0:401), NN2(0:401,0:401), NN3(0:401,0:401), / NU1(0:401,0:401), NU2(0:401,0:401), NU3(0:401,0:401), / TH1(0:401,0:401), TH2(0:401,0:401), TH3(0:401,0:401) INTEGER WRITE_ID COMMON /NUMER/ H1, H2, DT, N1, N2, EPS /ID/ WRITE_ID / /PHYS/ ALPHA, F, ZETA, THETA_0, M_FOLD, R / /ANIS/ a, b, c C functions of anisotropy C the anisotropy rotated 45 deg on Sep 4, 2006 c Psi (vv1,vv2,vv3) = DSQRT((1.D00+a)*vv1**2 + b*vv2**2 + c*vv3**2) c / +DSQRT(a*vv1**2 + (1.D00+b)*vv2**2 + c*vv3**2) c / +DSQRT(a*vv1**2 + b*vv2**2 + (1.D00+c)*vv3**2) c T0_1(vv1,vv2,vv3) = Psi (vv1,vv2,vv3)*vv1*( c / (1.D00+a)/DSQRT((1.D00+a)*vv1**2 + b*vv2**2 + c*vv3**2+E) c / +a/DSQRT(a*vv1**2 + (1.D00+b)*vv2**2 + c*vv3**2+E) c / +a/DSQRT(a*vv1**2 + b*vv2**2 + (1.D00+c)*vv3**2+E)) c T0_2(vv1,vv2,vv3) = Psi (vv1,vv2,vv3)*vv2*( c / b/DSQRT((1.D00+a)*vv1**2 + b*vv2**2 + c*vv3**2+E) c / +(1.D00+b)/DSQRT(a*vv1**2 + (1.D00+b)*vv2**2 + c*vv3**2+E) c / +b/DSQRT(a*vv1**2 + b*vv2**2 + (1.D00+c)*vv3**2+E)) c T0_3(vv1,vv2,vv3) = Psi (vv1,vv2,vv3)*vv3*( c / c/DSQRT((1.D00+a)*vv1**2 + b*vv2**2 + c*vv3**2+E) c / +c/DSQRT(a*vv1**2 + (1.D00+b)*vv2**2 + c*vv3**2+E) c / +(1.D00+c)/DSQRT(a*vv1**2 + b*vv2**2 + (1.D00+c)*vv3**2+E)) Psi (vv1,vv2,vv3) = / DSQRT((1.D00+a)*(vv1+vv2)**2 + b*(vv1-vv2)**2 + c*vv3**2) / +DSQRT(a*(vv1+vv2)**2 + (1.D00+b)*(vv1-vv2)**2 + c*vv3**2) / +DSQRT(a*(vv1+vv2)**2 + b*(vv1-vv2)**2 + (1.D00+c)*vv3**2) T0_1(vv1,vv2,vv3) = Psi (vv1,vv2,vv3)*( / ((1.D00+a)*(vv1+vv2)+b*(vv1-vv2)) / /DSQRT((1.D00+a)*(vv1+vv2)**2 + b*(vv1-vv2)**2 + c*vv3**2+E) / +(a*(vv1+vv2)+(1.D00+b)*(vv1-vv2)) / /DSQRT(a*(vv1+vv2)**2 + (1.D00+b)*(vv1-vv2)**2 + c*vv3**2+E) / +(a*(vv1+vv2)+b*(vv1-vv2)) / /DSQRT(a*(vv1+vv2)**2 + b*(vv1-vv2)**2 + (1.D00+c)*vv3**2+E)) T0_2(vv1,vv2,vv3) = Psi (vv1,vv2,vv3)*( / ((1.D00+a)*(vv1+vv2)-b*(vv1-vv2)) / /DSQRT((1.D00+a)*(vv1+vv2)**2 + b*(vv1-vv2)**2 + c*vv3**2+E) / +(a*(vv1+vv2)-(1.D00+b)*(vv1-vv2)) / /DSQRT(a*(vv1+vv2)**2 + (1.D00+b)*(vv1-vv2)**2 + c*vv3**2+E) / +(a*(vv1+vv2)-b*(vv1-vv2)) / /DSQRT(a*(vv1+vv2)**2 + b*(vv1-vv2)**2 + (1.D00+c)*vv3**2+E)) T0_3(vv1,vv2,vv3) = Psi (vv1,vv2,vv3)*vv3*( / c/DSQRT((1.D00+a)*(vv1+vv2)**2 + b*(vv1-vv2)**2 + c*vv3**2+E) / +c/DSQRT(a*(vv1+vv2)**2 + (1.D00+b)*(vv1-vv2)**2 + c*vv3**2+E) / +(1.D00+c)/ / DSQRT(a*(vv1+vv2)**2+b*(vv1-vv2)**2+(1.D00+c)*vv3**2+E)) E = 1.D-9 C reflexion of Phi C edges and corners DO I=0,N1 Phi(I,-1) = Phi(I,0) Phi(I,N2+1) = Phi(I,N2) ENDDO DO J=0,N2 Phi(-1,J) = Phi(0,J) Phi(N1+1,J) = Phi(N1,J) ENDDO Phi(-1,-1) = Phi(0,0) Phi(N1+1,-1) = Phi(N1,0) Phi(-1,N2+1) = Phi(0,N2) Phi(N1+1,N2+1)= Phi(N1,N2) C preparing gradients and normals DO i1=0,N1+1 DO i2=0,N2+1 v1(i1,i2) = (Phi(i1,i2)-Phi(i1-1,i2))/h1 v2(i1,i2) = (Phi(i1,i2)-Phi(i1,i2-1))/h2 C QQ(i1,i2) = DSQRT(eps**2+v1(i1,i2)**2+v2(i1,i2)**2) QQ(i1,i2) = Psi(-v1(i1,i2),-v2(i1,i2),1.D00) NU1(i1,i2) = -v1(i1,i2)/QQ(i1,i2) NU2(i1,i2) = -v2(i1,i2)/QQ(i1,i2) NU3(i1,i2) = 1.D00/QQ(i1,i2) NN1(i1,i2) = T0_1(NU1(i1,i2),NU2(i1,i2),NU3(i1,i2)) NN2(i1,i2) = T0_2(NU1(i1,i2),NU2(i1,i2),NU3(i1,i2)) NN3(i1,i2) = T0_3(NU1(i1,i2),NU2(i1,i2),NU3(i1,i2)) ENDDO ENDDO C computing curvature DO i1=0,N1 DO i2=0,N2 Kn(i1,i2) = / -(NN1(i1+1,i2)-NN1(i1,i2))/h1-(NN2(i1,i2+1)-NN2(i1,i2))/h2 ENDDO ENDDO C reflexion of curvature DO I=0,N1 Kn(I,-1) = Kn(I,0) Kn(I,N2+1) = Kn(I,N2) ENDDO DO J=0,N2 Kn(-1,J) = Kn(0,J) Kn(N1+1,J) = Kn(N1,J) ENDDO Kn(-1,-1) = Kn(0,0) Kn(N1+1,-1) = Kn(N1,0) Kn(-1,N2+1) = Kn(0,N2) Kn(N1+1,N2+1) = Kn(N1,N2) C preparing gradients C adding force according to Eberhard and previous codes in 1D DO i1=0,N1+1 DO i2=0,N2+1 v1(i1,i2) = (Kn(i1,i2)-Kn(i1-1,i2))/h1 / +F*(1.D00/Phi(i1,i2)-1.D00/Phi(i1-1,i2))/h1 v2(i1,i2) = (Kn(i1,i2)-Kn(i1,i2-1))/h2 / +F*(1.D00/Phi(i1,i2)-1.D00/Phi(i1,i2-1))/h2 TH1(i1,i2) = T0_1(v1(i1,i2),v2(i1,i2),0.0D00) TH2(i1,i2) = T0_2(v1(i1,i2),v2(i1,i2),0.0D00) TH3(i1,i2) = T0_3(v1(i1,i2),v2(i1,i2),0.0D00) ENDDO ENDDO C time evolution DO I=0,N1 DO J=0,N2 GG(I,J) = -(QQ(I+1,J)*(TH1(I+1,J)- / (TH1(I+1,J)*NU1(I+1,J)+TH2(I+1,J)*NU2(I+1,J) / +TH3(I+1,J)*NU3(I+1,J))*NN1(I+1,J))- / QQ(I,J)*(TH1(I,J)- / (TH1(I,J)*NU1(I,J)+TH2(I,J)*NU2(I,J) / +TH3(I,J)*NU3(I,J))*NN1(I,J)))/H1 / -(QQ(I,J+1)*(TH2(I,J+1)- / (TH1(I,J+1)*NU1(I,J+1)+TH2(I,J+1)*NU2(I,J+1) / +TH3(I,J+1)*NU3(I,J+1))*NN2(I,J+1))- / QQ(I,J)*(TH2(I,J)- / (TH1(I,J)*NU1(I,J)+TH2(I,J)*NU2(I,J) / +TH3(I,J)*NU3(I,J))*NN2(I,J)))/H2 ENDDO ENDDO RETURN END C-----------------------------------------------------------------------C SUBROUTINE CHECK_BLOWUP(Phi, IT) DOUBLE PRECISION Phi(-1:401,-1:401), H1, H2, DT, EPS, SIZE, BLOW COMMON /NUMER/ H1, H2, DT, N1, N2, EPS /CHECK/ BLOW SIZE = 0.D00 DO I=0,N1 DO J=0,N2 SIZE = DMAX1(SIZE, DABS(Phi(I,J))) ENDDO ENDDO IF (SIZE.GE.BLOW) THEN WRITE (*,*) 'There is a blow up observed' CALL SAVE_FIELD (Phi,IT) STOP ENDIF RETURN END C-----------------------------------------------------------------------C SUBROUTINE SAVE_FIELD (Phi, IT) DOUBLE PRECISION Phi(-1:401,-1:401), H1, H2, DT, EPS, / Kn(-1:401,-1:401), GRADX(0:400,0:400), GRADY(0:400,0:400), / GR_NORM(0:400,0:400), QQ(0:400,0:400), QQ_m(0:400,0:400), / Dx_0, Dy_0, Dx_p, Dy_p, Dx_m, Dy_m, Dxx, Dyy, Dxy, NABLA_p, / NABLA_m DOUBLE PRECISION F, DF, DDF, R, K, CURV(0:400,0:400), S, CURV_AUX CHARACTER*12 FNAME_NUM, FNAME_EXA CHARACTER*4 TLEV_INDX COMMON /NUMER/ H1, H2, DT, N1, N2, EPS C-----------------------------------------------------------------------C C entering the even extrapolation of Phi C edges and corners C This part is out of sense here DO I=0,N1 Phi(I,-1) = 2.D00*Phi(I,0)-Phi(I,1) Phi(I,N2+1) = 2.D00*Phi(I,N2)-Phi(I,N2-1) ENDDO DO J=0,N2 Phi(-1,J) = 2.D00*Phi(0,J)-Phi(1,J) Phi(N1+1,J) = 2.D00*Phi(N1,J)-Phi(N1-1,J) ENDDO Phi(-1,-1) = Phi(0,0) Phi(N1+1,-1) = Phi(N1,0) Phi(-1,N2+1) = Phi(0,N2) Phi(N1+1,N2+1)= Phi(N1,N2) C This maybe is useless DO I=0,N1 DO J=0,N2 C some auxiliary notations Dx_0 = (Phi(I+1,J)-Phi(I-1,J))/H1/2.D00 Dy_0 = (Phi(I,J+1)-Phi(I,J-1))/H2/2.D00 Dx_p = (Phi(I+1,J)-Phi(I,J))/H1 Dy_p = (Phi(I,J+1)-Phi(I,J))/H2 Dx_m = (Phi(I,J)-Phi(I-1,J))/H1 Dy_m = (Phi(I,J)-Phi(I,J-1))/H2 Dxx = (Phi(I+1,J)-2.D00*Phi(I,J)+Phi(I-1,J))/H1/H1 Dyy = (Phi(I,J+1)-2.D00*Phi(I,J)+Phi(I,J-1))/H2/H2 Dxy = ((Phi(I+1,J+1)-Phi(I-1,J+1))/H1/2.D00 / -(Phi(I+1,J-1)-Phi(I-1,J-1))/H1/2.D00)/H2/2.D00 C now saving what I need GRADX(I,J) = Dx_0 GRADY(I,J) = Dy_0 GR_NORM(I,J) = DSQRT(Dx_0*Dx_0+Dy_0*Dy_0) QQ(I,J) = DSQRT(Dx_0*Dx_0+Dy_0*Dy_0+EPS**2) C temporarily modify - desperately wanted to do smthng QQ_m(I,J) = (DSQRT(Dx_p*Dx_p+Dy_p*Dy_p+EPS**2)+ / DSQRT(Dx_m*Dx_m+Dy_m*Dy_m+EPS**2))/2.D00 C Done April 25, 2002 c IF (Dx_0**2+Dy_0**2.LE.1.D-20) THEN c Kn(I,J) = 0.0D00 c ELSE Kn(I,J) = / -((Dyy/QQ(I,J))*(Dx_0/QQ(I,J))*(Dx_0/QQ(I,J))+ / (Dxx/QQ(I,J))*(Dy_0/QQ(I,J))*(Dy_0/QQ(I,J)) / -2.D00*(Dx_0/QQ(I,J))*(Dy_0/QQ(I,J))*(Dxy/QQ(I,J))+ / (EPS/QQ(I,J))*(EPS/QQ(I,J))*((Dxx+Dyy)/QQ(I,J))) c ENDIF ENDDO ENDDO C entering the even extrapolation for curvature Kn C edges and corners DO I=0,N1 Kn(I,-1) = -Kn(I,1) Kn(I,N2+1) = -Kn(I,N2-1) ENDDO DO J=0,N2 Kn(-1,J) = -Kn(1,J) Kn(N1+1,J) = -Kn(N1-1,J) ENDDO Kn(-1,-1) = Kn(0,0) Kn(N1+1,-1) = Kn(N1,0) Kn(-1,N2+1) = Kn(0,N2) Kn(N1+1,N2+1)= Kn(N1,N2) C----------------------------------------------------------------------C encode (4,'(i4.4)', TLEV_INDX) IT FNAME_NUM = 'f_'//TLEV_INDX//'.dat' WRITE(*,*) 'Writing computed temperature field at time level = ', / IT,' to the file ', FNAME_NUM OPEN (77,FILE=FNAME_NUM,STATUS = 'UNKNOWN') 10 FORMAT (17H# Time level T = , F8.5, 2I10) 12 FORMAT (2x,20F5.1) WRITE (77,10) IT*DT DO I=0,N1 DO J=0,N2 WRITE (77,'(1X,2F15.5,2F25.5)') I*H1, J*H2, Phi(I,J), Kn(I,J) ENDDO WRITE (77,*) ENDDO CLOSE (77) RETURN END C-----------------------------------------------------------------------C C Engine shape SUBROUTINE INITIAL_LEVEL_engine (Phi) DOUBLE PRECISION Phi(-1:401,-1:401), H1, H2, DT, R, X, Y, NORM, / F, ZETA, THETA_0, EPS, ALPHA, SCALE, Pi, K COMMON /NUMER/ H1, H2, DT, N1, N2, EPS / /PHYS/ ALPHA, F, ZETA, THETA_0, M_FOLD, R / /PATTERN/ SCALE NORM(X,Y) = DSQRT(X*X+Y*Y) Pi = 3.1415926536D00 K = 2.D00 C Signed distance function smoothed by cut-off C-----------------------------------------------------------------------C C initial convex perturbation of the phase field C C Curve definition C C C C x_1 = R*cos(t), C C x_2 = R*sin(t). C C C C domain described by the condition [x_1-NUCL_1,x_2-NUCL_2] < curve C C H_crit_1, H_crit_2 have no meaning here C C sinusoidal perturbation C theta_0 is the frequency, pattern localized by tanh I0 = N1/2 J0 = N2/2 DO I=0,N1 DO J=0,N2 RR = NORM(0.0D00,DBLE(J-J0)*H2) Phi(I,J) = 1.D00 + / dsin(THETA_0*Pi*I/DBLE(N1))**M_fold* / (1.d00-4.d00*((J-J0)/DBLE(N2))**2) / *(1.D00-DTANH(ALPHA*(RR-R)))/2.D00 ENDDO ENDDO C-----------------------------------------------------------------------C C edges and corners DO I=0,N1 Phi(I,-1) = Phi(I,1) Phi(I,N2+1) = Phi(I,N2-1) ENDDO DO J=0,N2 Phi(-1,J) = Phi(1,J) Phi(N1+1,J) = Phi(N1-1,J) ENDDO Phi(-1,-1) = Phi(1,1) Phi(N1+1,-1) = Phi(N1-1,1) Phi(-1,N2+1) = Phi(1,N2-1) Phi(N1+1,N2+1)= Phi(N1-1,N2-1) RETURN END C-----------------------------------------------------------------------C C Brigitte shape SUBROUTINE INITIAL_LEVEL_Brigitte (Phi) DOUBLE PRECISION Phi(-1:401,-1:401), H1, H2, DT, R, X, Y, NORM, / F, ZETA, THETA_0, EPS, ALPHA, SCALE, Pi, K COMMON /NUMER/ H1, H2, DT, N1, N2, EPS / /PHYS/ ALPHA, F, ZETA, THETA_0, M_FOLD, R / /PATTERN/ SCALE NORM(X,Y) = DSQRT(X*X+Y*Y) Pi = 3.1415926536D00 K = 2.D00 C Signed distance function smoothed by cut-off C-----------------------------------------------------------------------C C initial convex perturbation of the phase field C C Curve definition C C C C x_1 = R*cos(t), C C x_2 = R*sin(t). C C C C domain described by the condition [x_1-NUCL_1,x_2-NUCL_2] < curve C C H_crit_1, H_crit_2 have no meaning here C C sinusoidal perturbation C theta_0 is the frequency, pattern localized by tanh I0 = N1/2 J0 = N2/2 DO I=0,N1 DO J=0,N2 RR = NORM(DBLE(I-I0)*H1,DBLE(J-J0)*H2) Phi(I,J) = 1.D00 + / dsin(THETA_0*Pi*I/DBLE(N1))**M_fold* / (1.d00-4.d00*((J-J0)/DBLE(N2))**2) / *(1.D00-DTANH(ALPHA*(RR-R)))/2.D00 ENDDO ENDDO C-----------------------------------------------------------------------C C edges and corners DO I=0,N1 Phi(I,-1) = Phi(I,1) Phi(I,N2+1) = Phi(I,N2-1) ENDDO DO J=0,N2 Phi(-1,J) = Phi(1,J) Phi(N1+1,J) = Phi(N1-1,J) ENDDO Phi(-1,-1) = Phi(1,1) Phi(N1+1,-1) = Phi(N1-1,1) Phi(-1,N2+1) = Phi(1,N2-1) Phi(N1+1,N2+1)= Phi(N1-1,N2-1) RETURN END C-----------------------------------------------------------------------C SUBROUTINE EXTRACT_LEVEL_SET (Phi, IT) DOUBLE PRECISION Phi(-1:401,-1:401), H1, H2, DT, EPS COMMON /NUMER/ H1, H2, DT, N1, N2, EPS OPEN (90,FILE='elast22.crv',STATUS='UNKNOWN',ACCESS='APPEND') 10 FORMAT (17H# Time level T = , F8.5, 2I10) 12 FORMAT (2x,20F5.1) C EPS = 0.005 WRITE (90,10) IT*DT call Identify_Curve (Phi) C write (*,*) 'Control LEVEL' write(90,'(/)') CLOSE (90) RETURN END C-----------------------------------------------------------------------C Subroutine Identify_Curve (field) DOUBLE PRECISION field(-1:401,-1:401), H1, H2, DT, EPS, Pi double precision f_max,f_min,cur_1(20000),cur_2(20000), / p_1(2),p_2(2) COMMON /NUMER/ H1, H2, DT, M1, M2, EPS integer M1,M2,k,l COMMON /LU/ Pi p_level = 0.D00 k = 0 DO I=1,M1 DO J=1,M2 f_max = dmax1(field(i,j),field(i-1,j),field(i,j-1), / field(i-1,j-1)) f_min = dmin1(field(i,j),field(i-1,j),field(i,j-1), / field(i-1,j-1)) if ((f_min-p_level)*(f_max-p_level).le.0.d00) then l=0 if ((field(i,j)-p_level)*(field(i-1,j)-p_level).lt.0.d00) / then l=l+1 p_1(l)=(i-1+(p_level-field(i-1,j))/ / (field(i,j)-field(i-1,j)))*h1 p_2(l)=j*h2 endif if (field(i,j)-p_level.eq.0.d00) then l=l+1 p_1(l)=i*h1 p_2(l)=j*h2 endif if (field(i-1,j)-p_level.eq.0.d00) then l=l+1 p_1(l)=(i-1)*h1 p_2(l)=j*h2 endif if ((field(i,j)-p_level)*(field(i,j-1)-p_level).lt.0.d00) / then l=l+1 if (l.eq.3) / write(*,*) 'Warning: curve many times in element' p_2(l)=(j-1+(p_level-field(i,j-1))/ / (field(i,j)-field(i,j-1)))*h2 p_1(l)=i*h1 endif if (field(i,j-1)-p_level.eq.0.d00) then l=l+1 if (l.eq.3) / write(*,*) 'Warning: curve many times in element' p_1(l)=i*h1 p_2(l)=(j-1)*h2 endif if ((field(i-1,j)-p_level)*(field(i-1,j-1)-p_level).lt.0.d00) / then l=l+1 if (l.eq.3) / write(*,*) 'Warning: curve many times in element' p_2(l)=(j-1+(p_level-field(i-1,j-1))/ / (field(i-1,j)-field(i-1,j-1)))*h2 p_1(l)=(i-1)*h1 endif if (field(i-1,j-1)-p_level.eq.0.d00) then l=l+1 if (l.eq.3) / write(*,*) 'Warning: curve many times in element' p_1(l)=(i-1)*h1 p_2(l)=(j-1)*h2 endif if ((field(i,j-1)-p_level)*(field(i-1,j-1)-p_level).lt.0.d00) / then l=l+1 if (l.eq.3) / write(*,*) 'Warning: curve many times in element' p_1(l)=(i-1+(p_level-field(i-1,j-1))/ / (field(i,j-1)-field(i-1,j-1)))*h1 p_2(l)=(j-1)*h2 endif k = k+1 if (l.eq.1) THEN cur_1(k) = p_1(1) cur_2(k) = p_2(1) else cur_1(k) = (p_1(1)+p_1(2))/2.0d00 cur_2(k) = (p_2(1)+p_2(2))/2.0d00 endif endif enddo enddo write (*,*) 'Total ', k, ' points identified.' call Parametrize_Curve (cur_1,cur_2,k) do i=1,k write (90,'(1x,2f10.5)') cur_1(i), cur_2(i) enddo return end C-----------------------------------------------------------------------C Subroutine Parametrize_Curve (pts_1,pts_2,number) integer number, i, i0, num_in, num_out double precision pts_1(20000),pts_2(20000),dist,outpts_1(20000), / outpts_2(20000) num_in = number-1 num_out = 1 outpts_1(1) = pts_1(number) outpts_2(1) = pts_2(number) do while (num_in.ge.1) dist = 10000.d00 do i=1,num_in if (dsqrt((outpts_1(num_out)-pts_1(i))**2+(outpts_2(num_out) / -pts_2(i))**2).lt.dist) then dist = dsqrt((outpts_1(num_out)-pts_1(i))**2+ / (outpts_2(num_out)-pts_2(i))**2) i0 = i endif enddo num_out = num_out+1 outpts_1(num_out) = pts_1(i0) outpts_2(num_out) = pts_2(i0) pts_1(i0) = pts_1(num_in) pts_2(i0) = pts_2(num_in) num_in = num_in-1 enddo do i=1,number pts_1(i) = outpts_1(i) pts_2(i) = outpts_2(i) enddo pts_1(number+1) = pts_1(1) pts_2(number+1) = pts_2(1) number = number+1 return end C-----------------------------------------------------------------------C C Hills shape SUBROUTINE INITIAL_LEVEL (Phi) DOUBLE PRECISION Phi(-1:401,-1:401), H1, H2, DT, R, X, Y, NORM, / F, ZETA, THETA_0, EPS, ALPHA, SCALE, Pi, K, R_ax, / THETA COMMON /NUMER/ H1, H2, DT, N1, N2, EPS / /PHYS/ ALPHA, F, ZETA, THETA_0, M_FOLD, R / /PATTERN/ SCALE /AX/ R_ax NORM(X,Y) = DSQRT(X*X+Y*Y) Pi = 3.1415926536D00 K = 2.D00 C Signed distance function smoothed by cut-off C-----------------------------------------------------------------------C C initial convex perturbation of the phase field C C Curve definition C C C C x_1 = R*cos(t), C C x_2 = R*sin(t). C C C C domain described by the condition [x_1-NUCL_1,x_2-NUCL_2] < curve C C H_crit_1, H_crit_2 have no meaning here C C sinusoidal perturbation I0 = N1/2 J0 = N2/2 DO I=-1,N1 DO J=-1,N2 THETA = DATAN2((J-J0)*H2,(I-I0)*H1+1.D-9) RR = NORM(DBLE(I-I0)*H1,DBLE(J-J0)*H2) C RR = NORM(DBLE(I-I0)*H1,0.D00) C RR = NORM(DBLE(I)*H1,0.D00)+SCALE*DCOS(2.*Pi*J/N2) C Phi(I,J) = 1.D00 + ZETA*DSIN(ALPHA*(RR-R))* C / DEXP(-ALPHA**2*(RR-R)**2)*(THETA_0+DSIGN(1.D00,RR-R)) C theta_0 is the mountain frequency, localized in a circle Phi(I,J) = 1.D00 + / K*DSIN(2.D00*PI*THETA_0*I/DBLE(N1))**M_FOLD* / DSIN(2.D00*PI*THETA_0*J/DBLE(N2))**M_FOLD* / (1.D00-DTANH(ALPHA*(RR-R)))/2.D00 ENDDO ENDDO C-----------------------------------------------------------------------C C edges and corners c DO I=0,N1 c Phi(I,-1) = Phi(I,1) c Phi(I,N2+1) = Phi(I,N2-1) c ENDDO c DO J=0,N2 c Phi(-1,J) = Phi(N1-1,J) c Phi(N1+1,J) = Phi(1,J) c ENDDO c Phi(-1,-1) = Phi(-1,1) c Phi(N1+1,-1) = Phi(N1+1,1) c Phi(-1,N2+1) = Phi(-1,N2-1) c Phi(N1+1,N2+1)= Phi(N1+1,N2-1) RETURN END C-----------------------------------------------------------------------C SUBROUTINE INITIAL_LEVEL_2D (Phi) DOUBLE PRECISION Phi(-1:401,-1:401), H1, H2, DT, R, X, Y, NORM, / F, ZETA, THETA_0, EPS, ALPHA, SCALE, Pi, K, R_ax, / THETA COMMON /NUMER/ H1, H2, DT, N1, N2, EPS / /PHYS/ ALPHA, F, ZETA, THETA_0, M_FOLD, R / /PATTERN/ SCALE /AX/ R_ax NORM(X,Y) = DSQRT(X*X+Y*Y) Pi = 3.1415926536D00 K = 2.D00 C Signed distance function smoothed by cut-off C-----------------------------------------------------------------------C C initial convex perturbation of the phase field C C Curve definition C C C C x_1 = R*cos(t), C C x_2 = R*sin(t). C C C C domain described by the condition [x_1-NUCL_1,x_2-NUCL_2] < curve C C H_crit_1, H_crit_2 have no meaning here C C sinusoidal perturbation I0 = N1/2 J0 = N2/2 DO I=-1,N1 DO J=-1,N2 THETA = DATAN2((J-J0)*H2,(I-I0)*H1+1.D-9) C RR = NORM(DBLE(I-I0)*H1,DBLE(J-J0)*H2) C RR = NORM(DBLE(I-I0)*H1,0.D00) RR = NORM(DBLE(I)*H1,0.D00)+SCALE*DCOS(2.*Pi*J/N2) Phi(I,J) = 1.D00 + ZETA*DSIN(ALPHA*(RR-R))* / DEXP(-ALPHA**2*(RR-R)**2)*(THETA_0+DSIGN(1.D00,RR-R)) ENDDO ENDDO C-----------------------------------------------------------------------C C edges and corners c DO I=0,N1 c Phi(I,-1) = Phi(I,1) c Phi(I,N2+1) = Phi(I,N2-1) c ENDDO c DO J=0,N2 c Phi(-1,J) = Phi(N1-1,J) c Phi(N1+1,J) = Phi(1,J) c ENDDO c Phi(-1,-1) = Phi(-1,1) c Phi(N1+1,-1) = Phi(N1+1,1) c Phi(-1,N2+1) = Phi(-1,N2-1) c Phi(N1+1,N2+1)= Phi(N1+1,N2-1) RETURN END C-----------------------------------------------------------------------C C Suitable for E0, E1 - standing and expanding circle SUBROUTINE INITIAL_LEVEL_wave (Phi) DOUBLE PRECISION Phi(-1:401,-1:401), H1, H2, DT, R, X, Y, NORM, / F, ZETA, THETA_0, EPS, ALPHA, SCALE, Pi, K, RR, / THETA COMMON /NUMER/ H1, H2, DT, N1, N2, EPS / /PHYS/ ALPHA, F, ZETA, THETA_0, M_FOLD, R / /PATTERN/ SCALE NORM(X,Y) = DSQRT(X*X+Y*Y) Pi = 3.1415926536D00 K = 2.D00 C Signed distance function smoothed by cut-off C-----------------------------------------------------------------------C C initial convex perturbation of the phase field C C Curve definition C C C C x_1 = R*cos(t), C C x_2 = R*sin(t). C C C C domain described by the condition [x_1-NUCL_1,x_2-NUCL_2] < curve C C H_crit_1, H_crit_2 have no meaning here C C sinusoidal perturbation I0 = N1/2 J0 = N2/2 DO I=0,N1 DO J=0,N2 c THETA = DATAN2((J-J0)*H2,(I-I0)*H1+1.D-9) c RR = NORM(DBLE(I-I0)*H1,DBLE(0)*H2)+R/5.D00* c / DSIN(M_FOLD*THETA) c Phi(I,J) = c / (0.5D00-DEXP(-K*(RR/R)**4)) C got from II_04/version20_f THETA = DATAN2((J-J0)*H2,(I-I0)*H1+1.D-9) RR = NORM(DBLE(I)*H1,0.D00) Phi(I,J) = 1.D00 + ZETA*DSIN(ALPHA*(RR-R))* / DEXP(-ALPHA**2*(RR-R)**2)*(THETA_0+DSIGN(1.D00,RR-R)) ENDDO ENDDO C-----------------------------------------------------------------------C C edges and corners DO I=0,N1 Phi(I,-1) = Phi(I,1) Phi(I,N2+1) = Phi(I,N2-1) ENDDO DO J=0,N2 Phi(-1,J) = Phi(1,J) Phi(N1+1,J) = Phi(N1-1,J) ENDDO Phi(-1,-1) = Phi(1,1) Phi(N1+1,-1) = Phi(N1-1,1) Phi(-1,N2+1) = Phi(1,N2-1) Phi(N1+1,N2+1)= Phi(N1-1,N2-1) RETURN END C-----------------------------------------------------------------------C SUBROUTINE TEST (Phi, IT) DOUBLE PRECISION Phi(-1:401,-1:401), H1, H2, DT, EPS, / Kn(-1:401,-1:401), GRADX(0:400,0:400), GRADY(0:400,0:400), / GR_NORM(0:400,0:400), QQ(0:400,0:400), QQ_m(0:400,0:400), / Dx_0, Dy_0, Dx_p, Dy_p, Dx_m, Dy_m, Dxx, Dyy, Dxy, NABLA_p, / NABLA_m DOUBLE PRECISION F, DF, DDF, R, K, CURV(0:400,0:400), S, CURV_AUX CHARACTER*12 FNAME_NUM, FNAME_EXA CHARACTER*3 TLEV_INDX COMMON /NUMER/ H1, H2, DT, N1, N2, EPS F(S) = 0.5D00-DEXP(-K*(S/R)**4) DF(S) = 4.0D00*K/R**4*S**2*DEXP(-K*(S/R)**4) DDF(S) = -DEXP(-K*(S/R)**4)*((4.0D00*K/R**4*S**3)**2- / 12.0D00*K/R**4*S**2) K = 1.D00 R = 0.2D00 C-----------------------------------------------------------------------C C entering the even extrapolation of Phi C edges and corners C This part is out the sense C redefining Phi as F DO I=0,N1 DO J=0,N2 S = DSQRT(((I-N1/2)*H1)**2+((J-N2/2)*H2)**2) Phi(I,J) = F(S) ENDDO ENDDO DO I=0,N1 Phi(I,-1) = 2.D00*Phi(I,0)-Phi(I,1) Phi(I,N2+1) = 2.D00*Phi(I,N2)-Phi(I,N2-1) ENDDO DO J=0,N2 Phi(-1,J) = 2.D00*Phi(0,J)-Phi(1,J) Phi(N1+1,J) = 2.D00*Phi(N1,J)-Phi(N1-1,J) ENDDO Phi(-1,-1) = Phi(0,0) Phi(N1+1,-1) = Phi(N1,0) Phi(-1,N2+1) = Phi(0,N2) Phi(N1+1,N2+1)= Phi(N1,N2) C This maybe is useless DO I=0,N1 DO J=0,N2 C some auxiliary notations Dx_0 = (Phi(I+1,J)-Phi(I-1,J))/H1/2.D00 Dy_0 = (Phi(I,J+1)-Phi(I,J-1))/H2/2.D00 Dx_p = (Phi(I+1,J)-Phi(I,J))/H1 Dy_p = (Phi(I,J+1)-Phi(I,J))/H2 Dx_m = (Phi(I,J)-Phi(I-1,J))/H1 Dy_m = (Phi(I,J)-Phi(I,J-1))/H2 Dxx = (Phi(I+1,J)-2.D00*Phi(I,J)+Phi(I-1,J))/H1/H1 Dyy = (Phi(I,J+1)-2.D00*Phi(I,J)+Phi(I,J-1))/H2/H2 Dxy = ((Phi(I+1,J+1)-Phi(I-1,J+1))/H1/2.D00 / -(Phi(I+1,J-1)-Phi(I-1,J-1))/H1/2.D00)/H2/2.D00 C now saving what I need GRADX(I,J) = Dx_0 GRADY(I,J) = Dy_0 GR_NORM(I,J) = DSQRT(Dx_0*Dx_0+Dy_0*Dy_0) QQ(I,J) = DSQRT(Dx_0*Dx_0+Dy_0*Dy_0+EPS**2) C temporarily modify - desperately wanted to do smthng QQ_m(I,J) = (DSQRT(Dx_p*Dx_p+Dy_p*Dy_p+EPS**2)+ / DSQRT(Dx_m*Dx_m+Dy_m*Dy_m+EPS**2))/2.D00 C Done April 25, 2002 c IF (Dx_0**2+Dy_0**2.LE.1.D-20) THEN c Kn(I,J) = 0.0D00 c ELSE Kn(I,J) = / -((Dyy/QQ(I,J))*(Dx_0/QQ(I,J))*(Dx_0/QQ(I,J))+ / (Dxx/QQ(I,J))*(Dy_0/QQ(I,J))*(Dy_0/QQ(I,J)) / -2.D00*(Dx_0/QQ(I,J))*(Dy_0/QQ(I,J))*(Dxy/QQ(I,J))+ / (EPS/QQ(I,J))*(EPS/QQ(I,J))*((Dxx+Dyy)/QQ(I,J))) c ENDIF S = DSQRT(((I-N1/2)*H1)**2+((J-N2/2)*H2)**2) CURV(I,J) = (-EPS**2*DDF(S)-DF(S)**3*S**2-EPS**2*DF(S))/ / (DSQRT(EPS**2+(DF(S)*S)**2))**3 ENDDO ENDDO C entering the even extrapolation for curvature Kn C edges and corners DO I=0,N1 Kn(I,-1) = -Kn(I,1) Kn(I,N2+1) = -Kn(I,N2-1) ENDDO DO J=0,N2 Kn(-1,J) = -Kn(1,J) Kn(N1+1,J) = -Kn(N1-1,J) ENDDO Kn(-1,-1) = Kn(0,0) Kn(N1+1,-1) = Kn(N1,0) Kn(-1,N2+1) = Kn(0,N2) Kn(N1+1,N2+1)= Kn(N1,N2) C----------------------------------------------------------------------C encode (3,'(i3.3)', TLEV_INDX) IT FNAME_NUM = 'f_'//TLEV_INDX//'.dat' WRITE(*,*) 'Writing computed temperature field at time level = ', / IT,' to the file ', FNAME_NUM OPEN (77,FILE=FNAME_NUM,STATUS = 'UNKNOWN') 10 FORMAT (17H# Time level T = , F8.5, 2I10) 12 FORMAT (2x,20F5.1) WRITE (77,10) IT*DT DO I=0,N1 DO J=0,N2 WRITE (77,'(1X,2F15.5,3F25.5)') I*H1, J*H2, Phi(I,J), / Kn(I,J), CURV(I,J) ENDDO WRITE (77,*) ENDDO CLOSE (77) OPEN (65,FILE='curv.dat',STATUS = 'UNKNOWN') DO I=0,N1 CURV_AUX = (-EPS**2*DDF(I*H1)-DF(I*H1)**3*(I*H1)**2- / EPS**2*DF(I*H1))/(DSQRT(EPS**2+(DF(I*H1)*I*H1)**2))**3 WRITE (65,'(1X,2F15.5)') I*H1,CURV_AUX ENDDO CLOSE (65) RETURN END C-----------------------------------------------------------------------C C loading and saving data SUBROUTINE STORAGE (Phi) DOUBLE PRECISION Phi(-1:401,-1:401), H1, H2, DT, R, X, Y, NORM, / F, ZETA, THETA_0, EPS, ALPHA, SCALE, Pi, K, / STEP INTEGER SAVE_ID, LOAD_ID COMMON /NUMER/ H1, H2, DT, N1, N2, EPS /MERSN/ STEP / /PHYS/ ALPHA, F, ZETA, THETA_0, M_FOLD, R / /PATTERN/ SCALE /STORE/ SAVE_ID, LOAD_ID IF (LOAD_ID.EQ.1) THEN LOAD_ID = 0 OPEN (77,FILE='elast22.ini',STATUS = 'UNKNOWN') READ (77,'(1X,F30.15)') STEP DO I=0,N1 READ (77,'(1X,20F30.15)') (Phi(I,J), J=0,N2) ENDDO CLOSE (77) WRITE (*,*) 'Data loaded from elast22.ini' ELSE OPEN (77,FILE='elast22.ini',STATUS = 'UNKNOWN') WRITE (77,'(1X,F30.15)') STEP DO I=0,N1 WRITE (77,'(1X,20F30.15)') (Phi(I,J), J=0,N2) ENDDO CLOSE (77) WRITE (*,*) 'Data saved into elast22.ini' ENDIF RETURN END C-----------------------------------------------------------------------C C New example with topology changes SUBROUTINE INITIAL_LEVEL_TOP (Phi) DOUBLE PRECISION Phi(-1:401,-1:401), H1, H2, DT, R, X, Y, NORM, / F, ZETA, THETA_0, EPS, ALPHA, SCALE, Pi, K, RR, / THETA COMMON /NUMER/ H1, H2, DT, N1, N2, EPS / /PHYS/ ALPHA, F, ZETA, THETA_0, M_FOLD, R / /PATTERN/ SCALE NORM(X,Y) = DSQRT(X*X+Y*Y) Pi = 3.1415926536D00 K = 2.D00 C Signed distance function smoothed by cut-off C-----------------------------------------------------------------------C C initial convex perturbation of the phase field C C Curve definition C C C C x_1 = R*cos(t), C C x_2 = R*sin(t). C C C C domain described by the condition [x_1-NUCL_1,x_2-NUCL_2] < curve C C H_crit_1, H_crit_2 have no meaning here C C sinusoidal perturbation I0 = N1/2 J0 = N2/2 DO I=0,N1 DO J=0,N2 THETA = DBLE(J)/DBLE(N2)*2.D00*Pi RR = NORM(DBLE(I-I0)*H1,0.D00)+ZETA* / DSIN(M_FOLD*THETA) Phi(I,J) = (0.5D00-DEXP(-K*(RR/R)**4)) ENDDO ENDDO C-----------------------------------------------------------------------C C edges and corners DO I=0,N1 Phi(I,-1) = Phi(I,1) Phi(I,N2+1) = Phi(I,N2-1) ENDDO DO J=0,N2 Phi(-1,J) = Phi(1,J) Phi(N1+1,J) = Phi(N1-1,J) ENDDO Phi(-1,-1) = Phi(1,1) Phi(N1+1,-1) = Phi(N1-1,1) Phi(-1,N2+1) = Phi(1,N2-1) Phi(N1+1,N2+1)= Phi(N1-1,N2-1) RETURN END C-----------------------------------------------------------------------C