From 52c5a0021ff733f5b90d423a8294392f7e3f599b Mon Sep 17 00:00:00 2001 From: tobias Date: Thu, 1 Apr 1999 23:24:39 +0000 Subject: These changes should allow IDLinearWaves to compile... git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDLinearWaves/trunk@9 5c0f84ea-6048-4d6e-bfa4-55cd5f2e0dd7 --- param.ccl | 8 +- schedule.ccl | 2 +- src/planewaves.F | 12 ++- src/teukwaves.F | 272 ++++++++++++++++++++++++++++--------------------------- 4 files changed, 149 insertions(+), 145 deletions(-) diff --git a/param.ccl b/param.ccl index 3acb784..5a8b69f 100644 --- a/param.ccl +++ b/param.ccl @@ -3,8 +3,8 @@ friend:einstein EXTENDS KEYWORD initial_data "" -"teukwaves" : "linear waves initialdata: Teukolsky waves" -"planewave" : "linear waves initialdata: plane waves" +"teukwaves" :: "linear waves initialdata: Teukolsky waves" +"planewave" :: "linear waves initialdata: plane waves" } private: @@ -66,6 +66,6 @@ KEYWORD wavesgoing "in and outgoing waves..." KEYWORD teuk_no_vee "Initialize Teuk. waves with V=0?" { - "no" : " Bona Masso setting" - "yes": " Bona Masso setting" + "no" :: "Bona Masso setting" + "yes":: "Bona Masso setting" } "no" diff --git a/schedule.ccl b/schedule.ccl index b2f0d5d..f7b878f 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -1,6 +1,6 @@ # Schedule definitions for thorn IDLinearWaves -if (CCTK_Equals(initaldata,"linearwaves")) +if (CCTK_Equals(initial_data,"linearwaves")) { schedule LinearWaves at CCTK_INITIAL { diff --git a/src/planewaves.F b/src/planewaves.F index 17b8431..683de4c 100644 --- a/src/planewaves.F +++ b/src/planewaves.F @@ -46,7 +46,7 @@ INTEGER iin,iout REAL pi - REAL amplitude,ra,the,phi + REAL ra,the,phi REAL wave,wavep REAL kx,ky,kz,w REAL kxp,kyp,kzp,wp @@ -54,11 +54,13 @@ c local arrays (scalars! Man, what a sweat!) REAL plus,minus,plusp,minusp,ain,aout - INTEGERE i,j,k + INTEGER i,j,k + + INTEGER CCTK_Equals pi = 3.14159265358979d0 - c Wave characteristics +c Wave characteristics c amplitude is declared by DECLARE_PARAMETER, for all other parameters c we stick to the naming in the cactus3.2.0 version of cactus @@ -106,7 +108,7 @@ c precalc kzp = 2*pi*cos(the)/wavep wp = (kxp*kxp+kyp*kyp+kzp*kzp)**0.5 - c *************** plane waves ******************** +c *************** plane waves ******************** do k=1,sh(3) do j=1,sh(2) do i=1,sh(1) @@ -141,7 +143,7 @@ c the metric functions c At this pint, the CAC3.2.0 version has some code for the dxgxx,.... c defined, how does CAC4.0 handles this without #ifdef THORN_BONAMASSO ? - c and finally the extrinsic curvatures +c and finally the extrinsic curvatures c Joan: I removed the division by alp as the initialization c of the lapse may take place later. This data is consistent c for initial lapse equal one. diff --git a/src/teukwaves.F b/src/teukwaves.F index c7db62a..a2f1ca4 100644 --- a/src/teukwaves.F +++ b/src/teukwaves.F @@ -30,23 +30,22 @@ subroutine teukwaves(CCTK_FARGUMENTS) implicit none - + DECLARE_CCTK_FARGUMENTS DECLARE_PARAMETERS - + REAL amp,m,ra,pi REAL wave,wave2,wave3,wave4,wave5,wave6,wave7,wave8 - REAL wavepulse INTEGER ipacket,iparity INTEGER ingoing,outgoing REAL kappa c point values of x.y z,r REAL xp,yp,zp,rp - + c spherical metric - REAL grr,grt,grp,gtt,gtp,gpp - REAL hrr,hrt,hrp,htt,htp,hpp + REAL teuk_grr,teuk_grt,teuk_grp,teuk_gtt,teuk_gtp,teuk_gpp + REAL teuk_hrr,teuk_hrt,teuk_hrp,teuk_htt,teuk_htp,teuk_hpp c special coefficeints needed for teukolsky waves c All these were 3d arrays in old versions of Newage,G,etc.. @@ -61,14 +60,17 @@ c from old spheretocart REAL drx,dry,drz,dtx,dty,dtz,dpx,dpy,dpz INTEGER i,j,k - + + INTEGER CCTK_Equals + cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc pi = 3.14159265358979 kappa = sqrt(143.0d0/pi)/12288.0d0 - + if (wavecenter.ne.0.) then - CCTK_Warn(0,"Teukwaves: wavecenter must be zero for time symmetry") + call CCTK_Warn(0,"Teukwaves: wavecenter must be + & zero for time symmetry") endif c CACTUS 4.0: all REAL/INTEGER parameter are available as variables and @@ -76,24 +78,24 @@ c need not be derived from the database. Nevertheless I will assign c the parameter variables to the variables of CAC3.0, to allow for c easy portation (wavepulse stays wavepulse btw). (GERD) - am = amplitude + amp = amplitude m = mvalue wave= wavelength - + c to be consistent with Evans & Abrahams, rescale the amplitude c by kappa, because they use that in the packet and c 4/3 because they use I not Teuk. F - if (CCTK_Equals(packet,'eppley')) then + if (CCTK_Equals(packet,'eppley') == 1) then ipacket = 1 write(*,*)'Teuk.Packet= eppley' - elseif (CCTK_Equals(packet,'evans')) then + elseif (CCTK_Equals(packet,'evans') == 1) then ipacket = 2 write(*,*)'Teuk.Packet= evans' c to be consistent with Evans & Abrahams, rescale the amplitude c by kappa, because they use that in the packet and c 4/3 because they use I not Teuk. F amp = amp*(4.0d0/3.0d0)*kappa - elseif (CCTK_Equals(packet,'square')) then + elseif (CCTK_Equals(packet,'square') == 1) then ipacket = 3 write(*,*)'Teuk.Packet= square' else @@ -107,25 +109,25 @@ c 4/3 because they use I not Teuk. F write(*,*)'wavecenter = ',ra write(*,*)'wavelength = ',wave - if (CCTK_Equals(parity,'even')) then + if (CCTK_Equals(parity,'even') == 1) then iparity = 0 - else (CONTAINS(parity,'odd')) then + elseif (CCTK_Equals(parity,'odd') == 1) then iparity = 1 endif - - if(CCTK_Equals(wavegoing,'in')) then - ingoing = 1 - endif - - if(CCTK_Equals(wavegoing,'out')) then + + if(CCTK_Equals(wavesgoing,'in') == 1) then + ingoing = 1 + endif + + if(CCTK_Equals(wavesgoing,'out') == 1) then outgoing = 1 endif - - if(CCTK_Equals(wavegoing,'both')) then + + if(CCTK_Equals(wavesgoing,'both') == 1) then ingoing = 1 outgoing = 1 endif - + if (ipacket.eq.1) then if (ingoing.eq.0.or.outgoing.eq.0) then write(*,*)'WARNING: Epply packet is only non-singular at' @@ -183,9 +185,9 @@ c to keep the analytic solution valid for the eppley packet: teuk_tn = (time-rp-ra) tn2 = teuk_tn**2 endif - + if (ingoing.eq.1) then - + tmp = exp(-tp2) fp = amp * teuk_tp*tmp fpa = amp * (1 - 2*tp2)*tmp @@ -194,23 +196,23 @@ c to keep the analytic solution valid for the eppley packet: fpd = amp * teuk_tp*(60 - 80*tp2 + 16*tp2*tp2)*tmp fpe = amp * ( & (60 - 360*tp2 + 240*tp2*tp2 - 32*tp2*tp2*tp2)*tmp ) - + if (outgoing.eq.0.and.(teuk_tp+wavepulse).le.0) & then - + fp = 0.0d0 fpa = 0.0d0 fpb = 0.0d0 fpc = 0.0d0 fpd = 0.0d0 fpe = 0.0d0 - + endif - + endif - + if (outgoing.eq.1) then - + tmp = exp(-tn2) fn = amp * teuk_tn*tmp fna = amp * (1 - 2*tn2)*tmp @@ -219,36 +221,36 @@ c to keep the analytic solution valid for the eppley packet: fnd = amp * teuk_tn*(60 - 80*tn2 + 16*tn2*tn2)*tmp fne = amp * ( & (60 - 360*tn2 + 240*tn2*tn2 - 32*tn2*tn2*tn2)*tmp ) - + endif endif - - + + c evans package --> w^4(1-x^2/w^2)^6 if(ipacket.eq.2) then - + if (ingoing.eq.1) then - + tmp = 1.0d0 - tp2/wave2 - + if (abs(teuk_tp).lt.wave) then - + fp = amp*wave5*tmp**6 fpa = -12.0d0*amp*wave3*teuk_tp*tmp**5 fpb = -12.0d0*amp*wave3*( & tmp**5 - 10.0d0*tp2/wave2*tmp**4 ) fpc = 120.0d0*amp*wave3*( - & 3.0d0*teuk_tp/wave2*tmp**4 + & 3.0d0*teuk_tp/wave2*tmp**4 & -8.0d0*teuk_tp*tp2/wave4*tmp**3 ) fpd = 360.0d0*amp*wave3*( - & 1.0d0/wave2*tmp**4 + & 1.0d0/wave2*tmp**4 & -16.0d0*tp2/wave4*tmp**3 & +16.0d0*tp2*tp2/wave6*tmp**2 ) fpe = 2880.0d0*amp*wave3*( & -5.0d0*teuk_tp/wave4*tmp**3 & +20.0d0*teuk_tp*tp2/wave6*tmp**2 & -8.0d0*teuk_tp*tp2*tp2/wave8*tmp ) - + else fp=0. fpa=0. @@ -257,31 +259,31 @@ c evans package --> w^4(1-x^2/w^2)^6 fpd=0. fpe=0. endif - + endif - + if (outgoing.eq.1) then - + tmp = 1.0d0 - tn2/wave2 - + if (abs(teuk_tn).lt.wave) then - + fn = amp*wave5*tmp**6 fna = -12.0d0*amp*wave3*teuk_tn*tmp**5 fnb = -12.0d0*amp*wave3*( & tmp**5 - 10.0d0*tn2/wave2*tmp**4 ) fnc = 120.0d0*amp*wave3*( - & 3.0d0*teuk_tn/wave2*tmp**4 + & 3.0d0*teuk_tn/wave2*tmp**4 & -8.0d0*teuk_tn*tn2/wave4*tmp**3 ) fnd = 360.0d0*amp*wave3*( - & 1.0d0/wave2*tmp**4 + & 1.0d0/wave2*tmp**4 & -16.0d0*tn2/wave4*tmp**3 & +16.0d0*tn2*tn2/wave6*tmp**2 ) fne = 2880.0d0*amp*wave3*( & -5.0d0*teuk_tn/wave4*tmp**3 & +20.0d0*teuk_tn*tn2/wave6*tmp**2 & -8.0d0*teuk_tn*tn2*tn2/wave8*tmp ) - + else fn=0. fna=0. @@ -290,17 +292,17 @@ c evans package --> w^4(1-x^2/w^2)^6 fnd=0. fne=0. endif - + endif - + endif c square package --> (1-x^2/w^2)^2 if(ipacket.eq.3) then - + write(*,*)'WARNING: need to calculate fpe,fne' write(*,*)'WARNING: need conditionals for in out waves' - + if (abs(teuk_tp).lt.wave) then fp = amp * (1-(tp2/wave2))**2 fpa = amp * (-4*teuk_tp)*(1-tp2/wave2)/wave2 @@ -329,7 +331,7 @@ c square package --> (1-x^2/w^2)^2 fnd=0. endif endif - + c c coefficients c @@ -347,7 +349,7 @@ c of the above coefficients. That is because the above coefficients c assume an ingoing and outgoing wave packet of the form: c wave = ingoing - outgoing. For a single wave, we want ingoing c and outgoing to have the same sign. - + if(ingoing.eq.0.and.outgoing.eq.1) then ca=-ca cb=-cb @@ -438,23 +440,23 @@ c mvalue c parity if (iparity.eq.0) then - grr = 1 + ca*frr - grt = cb*frt*rp - grp = cb*frp*(xp**2+yp**2)**0.5 - gtt = rp**2*(1 + cc*ftt1 + ca*ftt2) - gtp = (ca - 2*cc)*ftp*rp*(xp**2+yp**2)**0.5 - gpp = (xp**2+yp**2)*(1 + cc*fpp1 + ca*fpp2) + teuk_grr = 1 + ca*frr + teuk_grt = cb*frt*rp + teuk_grp = cb*frp*(xp**2+yp**2)**0.5 + teuk_gtt = rp**2*(1 + cc*ftt1 + ca*ftt2) + teuk_gtp = (ca - 2*cc)*ftp*rp*(xp**2+yp**2)**0.5 + teuk_gpp = (xp**2+yp**2)*(1 + cc*fpp1 + ca*fpp2) else - grr = 1. - grt = ck*drt*rp - grp = ck*drp*rp*sint - gtt = (1+cl*dtt)*rp**2 - gtp = cl*dtp*rp**2*sint - gpp = (1+cl*dpp)*rp**2*sint**2 + teuk_grr = 1. + teuk_grt = ck*drt*rp + teuk_grp = ck*drp*rp*sint + teuk_gtt = (1+cl*dtt)*rp**2 + teuk_gtp = cl*dtp*rp**2*sint + teuk_gpp = (1+cl*dpp)*rp**2*sint**2 endif - + c Finally, convert the spherical components to cartesian - + c define derivatives drx = (dr/dx) drx = xp/rp dry = yp/rp @@ -467,31 +469,34 @@ c define derivatives drx = (dr/dx) dpx = -yp/(xp**2+yp**2) dpy = xp/(xp**2+yp**2) dpz = 0. - + c cartesian g - gxx(i,j,k) = drx**2*grr+2.*drx*dtx*grt - & +2.*drx*dpx*grp+dtx**2*gtt - & +2.*dtx*dpx*gtp+dpx**2*gpp + gxx(i,j,k) = drx**2*teuk_grr+2.*drx*dtx*teuk_grt + & +2.*drx*dpx*teuk_grp+dtx**2*teuk_gtt + & +2.*dtx*dpx*teuk_gtp+dpx**2*teuk_gpp - gyy(i,j,k) = dry**2*grr+2.*dry*dty*grt - & +2.*dry*dpy*grp+dty**2*gtt - & +2.*dty*dpy*gtp+dpy**2*gpp + gyy(i,j,k) = dry**2*teuk_grr+2.*dry*dty*teuk_grt + & +2.*dry*dpy*teuk_grp+dty**2*teuk_gtt + & +2.*dty*dpy*teuk_gtp+dpy**2*teuk_gpp - gzz(i,j,k) = drz**2*grr+2.*drz*dtz*grt - & +2.*drz*dpz*grp+dtz**2*gtt - & +2.*dtz*dpz*gtp+dpz**2*gpp + gzz(i,j,k) = drz**2*teuk_grr+2.*drz*dtz*teuk_grt + & +2.*drz*dpz*teuk_grp+dtz**2*teuk_gtt + & +2.*dtz*dpz*teuk_gtp+dpz**2*teuk_gpp - gxy(i,j,k) = drx*dry*grr+(drx*dty+dtx*dry)*grt - & +(drx*dpy+dpx*dry)*grp - & +dtx*dty*gtt+(dtx*dpy+dpx*dty)*gtp+dpx*dpy*gpp + gxy(i,j,k) = drx*dry*teuk_grr+(drx*dty+dtx*dry)*teuk_grt + & +(drx*dpy+dpx*dry)*teuk_grp + & +dtx*dty*teuk_gtt+(dtx*dpy+dpx*dty)*teuk_gtp+ + & dpx*dpy*teuk_gpp - gxz(i,j,k) = drx*drz*grr+(drx*dtz+dtx*drz)*grt - & +(drx*dpz+dpx*drz)*grp - & +dtx*dtz*gtt+(dtx*dpz+dpx*dtz)*gtp+dpx*dpz*gpp + gxz(i,j,k) = drx*drz*teuk_grr+(drx*dtz+dtx*drz)*teuk_grt + & +(drx*dpz+dpx*drz)*teuk_grp + & +dtx*dtz*teuk_gtt+(dtx*dpz+dpx*dtz)*teuk_gtp+ + & dpx*dpz*teuk_gpp - gyz(i,j,k) = dry*drz*grr+(dry*dtz+dty*drz)*grt - & +(dry*dpz+dpy*drz)*grp - & +dty*dtz*gtt+(dty*dpz+dpy*dtz)*gtp+dpy*dpz*gpp + gyz(i,j,k) = dry*drz*teuk_grr+(dry*dtz+dty*drz)*teuk_grt + & +(dry*dpz+dpy*drz)*teuk_grp + & +dty*dtz*teuk_gtt+(dty*dpz+dpy*dtz)*teuk_gtp+ + & dpy*dpz*teuk_gpp @@ -506,64 +511,61 @@ c for the ext. curv. well need... ckdot = (fpc-fnc)/rp**2-3*(fpb+fnb)/rp**3+3*(fpa-fna)/rp**4 cldot = -(fpd+fnd)/rp+2*(fpc-fnc)/rp**2-3*(fpb+fnb)/rp**3+ & 3*(fpa-fna)/rp**4 - + c parity if (iparity.eq.0) then - hrr = -0.5*(cadot*frr) - hrt = -0.5*(cbdot*frt*rp) - hrp = -0.5*(cbdot*frp*(xp**2+yp**2)**0.5) - htt = -0.5*(rp**2*(ccdot*ftt1 + cadot*ftt2)) - htp = -0.5*((cadot - 2*ccdot)*ftp*rp*(xp**2+yp**2)**0.5) - hpp = -0.5*((xp**2+yp**2)*(ccdot*fpp1 + cadot*fpp2)) + teuk_hrr = -0.5*(cadot*frr) + teuk_hrt = -0.5*(cbdot*frt*rp) + teuk_hrp = -0.5*(cbdot*frp*(xp**2+yp**2)**0.5) + teuk_htt = -0.5*(rp**2*(ccdot*ftt1 + cadot*ftt2)) + teuk_htp = + & -0.5*((cadot - 2*ccdot)*ftp*rp*(xp**2+yp**2)**0.5) + teuk_hpp = -0.5*((xp**2+yp**2)*(ccdot*fpp1 + cadot*fpp2)) else - hrr = 0. - hrt = -0.5*(ckdot*drt*rp) - hrp = -0.5*(ckdot*drp*rp*sint) - htt = -0.5*((cldot*dtt)*rp**2) - htp = -0.5*(cldot*dtp*rp**2*sint) - hpp = -0.5*((cldot*dpp)*rp**2*sint**2) + teuk_hrr = 0. + teuk_hrt = -0.5*(ckdot*drt*rp) + teuk_hrp = -0.5*(ckdot*drp*rp*sint) + teuk_htt = -0.5*((cldot*dtt)*rp**2) + teuk_htp = -0.5*(cldot*dtp*rp**2*sint) + teuk_hpp = -0.5*((cldot*dpp)*rp**2*sint**2) endif c time symmetry - hxx(i,j,k) = drx**2*hrr+2.*drx*dtx*hrt - & +2.*drx*dpx*hrp+dtx**2*htt - & +2.*dtx*dpx*htp+dpx**2*hpp - - hyy(i,j,k) = dry**2*hrr+2.*dry*dty*hrt - & +2.*dry*dpy*hrp+dty**2*htt - & +2.*dty*dpy*htp+dpy**2*hpp + hxx(i,j,k) = drx**2*teuk_hrr+2.*drx*dtx*teuk_hrt + & +2.*drx*dpx*teuk_hrp+dtx**2*teuk_htt + & +2.*dtx*dpx*teuk_htp+dpx**2*teuk_hpp - hzz(i,j,k) = drz**2*hrr+2.*drz*dtz*hrt - & +2.*drz*dpz*hrp+dtz**2*htt - & +2.*dtz*dpz*htp+dpz**2*hpp + hyy(i,j,k) = dry**2*teuk_hrr+2.*dry*dty*teuk_hrt + & +2.*dry*dpy*teuk_hrp+dty**2*teuk_htt + & +2.*dty*dpy*teuk_htp+dpy**2*teuk_hpp - hxy(i,j,k) = drx*dry*hrr+(drx*dty+dtx*dry)*hrt - & +(drx*dpy+dpx*dry)*hrp - & +dtx*dty*htt+(dtx*dpy+dpx*dty)*htp+dpx*dpy*hpp + hzz(i,j,k) = drz**2*teuk_hrr+2.*drz*dtz*teuk_hrt + & +2.*drz*dpz*teuk_hrp+dtz**2*teuk_htt + & +2.*dtz*dpz*teuk_htp+dpz**2*teuk_hpp + + hxy(i,j,k) = drx*dry*teuk_hrr+(drx*dty+dtx*dry)*teuk_hrt + & +(drx*dpy+dpx*dry)*teuk_hrp + & +dtx*dty*teuk_htt+(dtx*dpy+dpx*dty)*teuk_htp+ + & dpx*dpy*teuk_hpp + + hxz(i,j,k) = drx*drz*teuk_hrr+(drx*dtz+dtx*drz)*teuk_hrt + & +(drx*dpz+dpx*drz)*teuk_hrp + & +dtx*dtz*teuk_htt+(dtx*dpz+dpx*dtz)*teuk_htp+ + & dpx*dpz*teuk_hpp + + hyz(i,j,k) = dry*drz*teuk_hrr+(dry*dtz+dty*drz)*teuk_hrt + & +(dry*dpz+dpy*drz)*teuk_hrp + & +dty*dtz*teuk_htt+(dty*dpz+dpy*dtz)*teuk_htp+ + & dpy*dpz*teuk_hpp - hxz(i,j,k) = drx*drz*hrr+(drx*dtz+dtx*drz)*hrt - & +(drx*dpz+dpx*drz)*hrp - & +dtx*dtz*htt+(dtx*dpz+dpx*dtz)*htp+dpx*dpz*hpp - hyz(i,j,k) = dry*drz*hrr+(dry*dtz+dty*drz)*hrt - & +(dry*dpz+dpy*drz)*hrp - & +dty*dtz*htt+(dty*dpz+dpy*dtz)*htp+dpy*dpz*hpp - - -c$$$ hxx(i,j,k) = 0. -c$$$ hxy(i,j,k) = 0. -c$$$ hxz(i,j,k) = 0. -c$$$ hyy(i,j,k) = 0. -c$$$ hyz(i,j,k) = 0. -c$$$ hzz(i,j,k) = 0. - - enddo enddo enddo - - + return - end + end + + -- cgit v1.2.3