diff options
Diffstat (limited to 'src/bhbrill3d.m')
-rw-r--r-- | src/bhbrill3d.m | 132 |
1 files changed, 132 insertions, 0 deletions
diff --git a/src/bhbrill3d.m b/src/bhbrill3d.m new file mode 100644 index 0000000..ef7bccc --- /dev/null +++ b/src/bhbrill3d.m @@ -0,0 +1,132 @@ +$Path = Union[$Path,{"~/SetTensor"}]; +Needs["SetTensor`"]; + +Dimension = 3; +x[1] = eta; x[2] = q; x[3] = phi +qf[eta_,q_,phi_] := amp (Exp[-(eta-eta0)^2/sigma^2]+Exp[-(eta+eta0)^2/sigma^2]) Sin[q]^n (1+c Cos[phi]^2) + +md = { +{Exp[2 qf[eta,q,phi]],0,0}, +{0,Exp[2 qf[eta,q,phi]],0}, +{0,0,Sin[q]^2}} psisph[eta,q,phi]^4; +InitializeMetric[md]; + +Clear[exc]; +DefineTensor[exc]; +SetTensor[exc[la,lb],{{0,0,0},{0,0,0},{0,0,0}}]; + +tmp = RicciR[la,lb] Metricg[ua,ub]+exc[la,lb] Metricg[ua,ub] exc[lc,ld] Metricg[uc,ud]- + exc[la,lb] exc[lc,ld] Metricg[ua,uc] Metricg[ub,ud]; +tmp = RicciToAffine[tmp]; +tmp = EvalMT[tmp]; +tmp = ExpandAll[-Exp[2 qf[eta,q,phi]] psisph[eta,q,phi]^5/8 tmp] +sav=tmp +tmp = SubFun[sav,psisph[eta,q,phi],2 Cosh[eta/2]+psisph[eta,q,phi]] + +(* Make the stencil... *) + +stencil = ExpandAll[tmp /. { + D[psisph[eta,q,phi],eta]->(psisph[i+1,j,k]-psisph[i-1,j,k])/(2 deta), + D[psisph[eta,q,phi],eta,eta]->(psisph[i+1,j,k]+psisph[i-1,j,k]-2 psisph[i,j,k])/(deta deta), + D[psisph[eta,q,phi],q]->(psisph[i,j+1,k]-psisph[i,j-1,k])/(2 dq), + D[psisph[eta,q,phi],q,q]->(psisph[i,j+1,k]+psisph[i,j-1,k]-2 psisph[i,j,k])/(dq dq), + D[psisph[eta,q,phi],phi]->(psisph[i,j,k+1]-psisph[i,j,k-1])/(2 dphi), + D[psisph[eta,q,phi],phi,phi]->(psisph[i,j,k+1]+psisph[i,j,k-1]-2 psisph[i,j,k])/(dphi dphi), + psisph[eta,q,phi]->psisph[i,j,k] + }]; + +ac = Coefficient[stencil,psisph[i,j,k]] +an = Coefficient[stencil,psisph[i+1,j,k]] +as = Coefficient[stencil,psisph[i-1,j,k]] +ae = Coefficient[stencil,psisph[i,j,k+1]] +aw = Coefficient[stencil,psisph[i,j,k-1]] +aq = Coefficient[stencil,psisph[i,j+1,k]] +ab = Coefficient[stencil,psisph[i,j-1,k]] +rhs = -SubFun[tmp,psisph[eta,q,phi],0] + +FortranOutputOfDepList = "(i,j,k)"; +$FortranReplace = Union[{ + "UND"->"_", + "(eta,q,phi)"->"(i,j,k)" +}]; +fd = FortranOpen["bhbrill3d.x"]; +FortranWrite[fd,An[i,j,k],an ]; +FortranWrite[fd,As[i,j,k],as ]; +FortranWrite[fd,Ae[i,j,k],ae ]; +FortranWrite[fd,Aw[i,j,k],aw ]; +FortranWrite[fd,Aq[i,j,k],aq ]; +FortranWrite[fd,Ab[i,j,k],ab ]; +FortranWrite[fd,Ac[i,j,k],ac ]; +FortranWrite[fd,Rhs[i,j,k],rhs ]; +FortranClose[fd]; + +(* Next part, write out conformal g's and d's *) + + +xv = Exp[eta] Sin[q] Cos[phi]; +yv = Exp[eta] Sin[q] Sin[phi]; +zv = Exp[eta] Cos[q]; + +mc = Table[ D[ {xv,yv,zv}[[i]], {eta,q,phi}[[j]] ],{i,1,3},{j,1,3}]; +mci = Simplify[Inverse[mc]]; + +Clear[mct]; +DefineTensor[mct,{{1,2},1}]; +Iter[mct[ua,lb], + mct[ua,lb]=mc[[ua,-lb]]; +]; + +Clear[mcti]; +DefineTensor[mcti,{{1,2},1}]; +Iter[mcti[ua,lb], + mcti[ua,lb]=mci[[ua,-lb]]; +]; + +gijtmp = Exp[2 eta]/psisph[eta,q,phi]^4 Metricg[lc,ld] mcti[uc,la] mcti[ud,lb] + +Clear[i2]; +DefineTensor[i2,{{2,1},1}]; + +fd = FortranOpen["gij.x"]; +Iter[i2[ua,ub], + v1 = {x,y,z}[[ua]]; + v2 = {x,y,z}[[ub]]; + metv = ToExpression["g"<>ToString[v1]<>ToString[v2]<>"[i,j,k]"]; + gg[v1,v2]=Simplify[EvalMT[gijtmp,{la->-ua,lb->-ub}]]; + FortranWrite[fd,metv,gg[v1,v2]]; + For[ii=1,ii<=3,ii++, + v3 = {x,y,z}[[ii]]; + dmetv = ToExpression["d"<>ToString[v3]<>ToString[metv]]; + res = OD[gg[v1,v2],lc] mcti[uc,ld]/2; + res = EvalMT[res,ld-> -ii]; + res = Simplify[res]; + FortranWrite[fd,dmetv,res]; + ]; +]; +FortranClose[fd]; + +$FortranReplace = { + "UND"->"_", + "(eta,q,phi)"->"" +}; + +fd = FortranOpen["psi_1st_deriv.x"]; +For[ii=1,ii<=3,ii++, + v1 = {x,y,z}[[ii]]; + psv =ToExpression["psi"<>ToString[v1]<>"[i,j,k]"]; + rhs = CD[Exp[-eta/2] psi3d[eta,q,phi],lc] mcti[uc,la]; + rhs = EvalMT[rhs,{la->-ii}]/(Exp[-eta/2] psi3d[eta,q,phi]); + FortranWrite[fd,psv,rhs]; +]; +FortranClose[fd]; + +fd = FortranOpen["psi_2nd_deriv.x"]; +Iter[i2[ua,ub], + v1 = {x,y,z}[[ua]]; + v2 = {x,y,z}[[ub]]; + psv = ToExpression["psi"<>ToString[v1]<>ToString[v2]<>"[i,j,k]"]; + rhs = OD[OD[Exp[-eta/2] psi3d[eta,q,phi],lc] mcti[uc,la],ld] mcti[ud,lb]; + rhs = EvalMT[rhs,{la->-ua,lb->-ub}]/(Exp[-eta/2] psi3d[eta,q,phi]); + FortranWrite[fd,psv,rhs]; +]; +FortranClose[fd]; |