Guest User

Log of FreeFEM++ for cylindrical waveguide

a guest
Jul 12th, 2023
163
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 4.64 KB | None | 0 0
  1. do script " cd '.'; /usr/local/ff++/mpich3/bin/FreeFem++ 'test.edp' -glut '/usr/local/ff++/mpich3/bin/ffglut' ;"
  2. -- FreeFem++ v4.9 (Mer 21 avr 2021 16:09:39 CEST - git v4.9)
  3. Load: lg_fem lg_mesh lg_mesh3 eigenvalue
  4. 1 : load "MUMPS"init MUMPS_SEQ: MPI_Init
  5.  
  6. 2 : load "iovtk" load: iovtk
  7.  
  8. 3 : load "medit"
  9. 4 : load "msh3"
  10. 5 : load "Element_Mixte3d"
  11. 6 :
  12. 7 : bool debug = true;
  13. 8 :
  14. 9 : real r = 0.00127;
  15. 10 : real d = r;
  16. 11 : int nx = 20;
  17. 12 : int nz = 20;
  18. 13 : real zmin = 0.0;
  19. 14 : real zmax = 0.01;
  20. 15 :
  21. 16 : // build the mesh
  22. 17 : border C(t=0.0, 2.0*pi) {
  23. 18 : x=r*cos(t);
  24. 19 : y=r*sin(t);
  25. 20 : label=1;
  26. 21 : };
  27. 22 : mesh Thx = buildmesh(C(nx));
  28. 23 :
  29. 24 : int [int] rup = [0, 1], rdown = [0, 2], rmid = [1, 3];
  30. 25 : int guide = 3, in = 1, out = 2;
  31. 26 : mesh3 Th = buildlayers(Thx, nz, zbound=[zmin, zmax]
  32. 27 : , coef = 1.
  33. 28 : , labelup = rup
  34. 29 : , labelmid = rmid
  35. 30 : , labeldown = rdown
  36. 31 : );
  37. 32 : //medit("Cyl", Th, wait=debug);
  38. 33 : //plot(Th, wait=debug);
  39. 34 :
  40. 35 : cout << "labels(Th) = " << labels(Th) << endl;
  41. 36 :
  42. 37 : real f = 94*10^9;
  43. 38 : real er = 1;
  44. 39 : real c0 = 299792458;
  45. 40 : real k = 2*pi*f*sqrt(er) / c0;
  46. 41 : int m = 1, n = 1;
  47. 42 : real beta = sqrt(k^2 - (m*pi/d)^2 - (n*pi/d)^2);
  48. 43 : real Z = sqrt(er) * 120 * pi;
  49. 44 :
  50. 45 : real ukb = 1 / (k^2 - beta^2);
  51. 46 : func expbz = exp(-1i*beta*z);
  52. 47 : func ExTE = (1i*k*Z) * ukb * (n*pi) / d * cos(m*pi*x/d) * sin(n*pi*y/d) * expbz;
  53. 48 : func EyTE = -(1i*k*Z) * ukb * (n*pi) / d * sin(m*pi*x/d) * cos(n*pi*y/d) * expbz;
  54. 49 :
  55. 50 : func Gix = +2*1i*beta*ExTE;
  56. 51 : func Giy = +2*1i*beta*EyTE;
  57. 52 : func Giz = 0;
  58. 53 :
  59. 54 :
  60. 55 : fespace Nh(Th, Edge13d);
  61. 56 : Nh<complex> [Ex, Ey, Ez], [vx, vy, vz];
  62. 57 :
  63. 58 : // Macros
  64. 59 : macro Curl(ux,uy,uz) [dy(uz)-dz(uy),dz(ux)-dx(uz),dx(uy)-dy(ux)] ) // EOM
  65. 60 : macro Nvec(ux,uy,uz) [uy*N.z-uz*N.y,uz*N.x-ux*N.z,ux*N.y-uy*N.x] ) // EOM //uxN
  66. 61 : macro Curlabs(ux,uy,uz) [abs(dy(uz)-dz(uy)),abs(dz(ux)-dx(uz)),abs(dx(uy)-dy(ux))] ) //EOM
  67. 62 :
  68. 63 : problem waveguide([Ex, Ey, Ez], [vx, vy, vz], solver=sparsesolver) =
  69. 64 : int3d(Th)(Curl(Ex,Ey,Ez) [dy(Ez)-dz(Ey),dz(Ex)-dx(Ez),dx(Ey)-dy(Ex)] ' * Curl(vx,vy,vz) [dy(vz)-dz(vy),dz(vx)-dx(vz),dx(vy)-dy(vx)] )
  70. 65 : - int3d(Th)(k^2*[Ex,Ey,Ez]' * [vx, vy, vz])
  71. 66 : + int2d(Th,in,out)(1i * beta * Nvec(Ex,Ey,Ez) [Ey*N.z-Ez*N.y,Ez*N.x-Ex*N.z,Ex*N.y-Ey*N.x] ' * Nvec(vx,vy,vz) [vy*N.z-vz*N.y,vz*N.x-vx*N.z,vx*N.y-vy*N.x] )
  72. 67 : - int2d(Th,in)([vx,vy,vz]' * [Gix,Giy,Giz])
  73. 68 : + on(guide,Ex=0,Ey=0,Ez=0);
  74. 69 :
  75. 70 : waveguide;
  76. 71 :
  77. 72 : cout << "r = " << r << " Integrated Ex = " << int3d(Th)(abs(Ex)^2) << endl;
  78. 73 : medit("real", Th, [real(Ex),real(Ey),real(Ez)]);
  79. 74 : //cout << (abs(Ex)^2) << end;
  80. 75 : //cout << (abs(Ey)^2) << end;
  81. 76 : //cout << (abs(Ez)^2) << end;
  82. 77 :
  83. 78 : //int[int] Order = [1];
  84. 79 : //savevtk("result_x.vtu", Th, abs(real(Ex)), dataname="RealEx", order=Order);
  85. 80 : //savevtk("result_y.vtu", Th, abs(real(Ey)), dataname="RealEy", order=Order);
  86. 81 : //savevtk("result_z.vtu", Th, abs(real(Ez)), dataname="RealEz", order=Order);
  87. 82 : //medit("real", Th, [real(Ex),real(Ey),real(Ez)]);
  88. 83 : sizestack + 1024 =9384 ( 8360 )
  89.  
  90. -- mesh: Nb of Triangles = 68, Nb of Vertices 45
  91. labels(Th) = 3
  92. 1 2 3
  93. -- Build Nodes/DF on mesh : n.v. 945, n. elmt. 4080, n b. elmt. 936
  94. nb of Nodes 14120 nb of DoF 28240 DFon=0220
  95. -- FESpace: Nb of Nodes 14120 Nb of DoF 28240
  96. -- Solve :
  97. min (0,0) max (0,0)
  98. -- Medit, Release 3.0a (Nov. 30, 2007)
  99. Copyright (c) LJLL, 1999-2007.
  100. compiled: Mer 20 mai 2020 16:08:51 CEST (with ff++ 4.6).
  101.  
  102. medit with binary version of popen : mesh(es) and solution(s)
  103. mesh_name= real
  104. Loading data file(s)
  105. End of mesh
  106. .sol: Dimension 3 (mesh)3 (lecture)1
  107. Input seconds: 0.00
  108.  
  109. medit1()
  110.  
  111. Building scene(s)
  112. Creating scene 1
  113. Loading default options
  114. Scene seconds: 0.03
  115.  
  116. Rendering scene(s)
  117.  
  118. Total running seconds: 0.40
  119. Thank you for using Medit.
  120. r = 0.00127 Integrated Ex = nan
  121. version de medit ffmedit -popen -filebin -addsol 1 real
  122. times: compile 0.022873s, execution 8.47454s, mpirank:0
  123. ######## We forget of deleting 6 Nb pointer, 0Bytes , mpirank 0, memory leak =44677024
  124. CodeAlloc : nb ptr 4277, size :496712 mpirank: 0
  125. Ok: Normal End
  126. close MUMPS_SEQ: MPI_Finalize
  127.  
Tags: FreeFEM++
Advertisement
Add Comment
Please, Sign In to add comment