Skip to content
Snippets Groups Projects
Commit 023c9150 authored by Nicolas Marsic's avatar Nicolas Marsic
Browse files

remove old photon trap stuff

parent a24f61c9
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
tests/haroche/cavity_haroche_getdp.py deleted 100644 → 0
+ 0
− 101
View file @ a24f61c9
import os
import scipy as sc
import scipy.io as sio
import numpy as np
import pylab as pl
import sys
pi = np.pi
# os.system('rm '+str_path+'/*.pos')
# os.system('rm '+str_path+'/*.dat')
# os.system('rm '+str_path+'/*.out')
# os.system('rm '+str_path+'/*.txt')
# os.system('rm '+str_path+'/*.pre')
# os.system('rm '+str_path+'/*.res')
#### your configuration (nproc is your number of available processors)
str_gmsh_path = '~/programs/gmsh/Gmsh-2.8.4-stable.app/Contents/MacOS/'
str_getdp_path = '~/programs/getdp/getdp_build_svn/'
str_mesh_filename = 'haroche.msh'
str_pro_filename = 'cavity_haroche.pro'
#### scaling the problem (*10^12)
nm = 1.e-9
epsilon0 = 8.854187817e-3*nm
mu0 = 400.*pi*nm
cel = 1.0/(np.sqrt(epsilon0 * mu0))
#### mesh parameters : MeshElementSize = lambda/(paramaille*n)
lambda_m = 5.8e-3
paramaille_air = lambda_m/3.
paramaille_pml = lambda_m/3.
paramaille_mir = lambda_m/3.
#### opto-geometric parameters
#lambda_vp = 5.8e-3
#eig_target= (2.*np.pi*cel/lambda_vp)**2
freq_target = 51.099e9
lambda_vp = cel/freq_target
eig_target= (2.*np.pi*cel/lambda_vp)**2
neig = 70
#### write all above variables in a text file, used as an input to both gmsh and getdp
par_gmsh_getdp =open('parameters_gmsh_getdp.dat', 'w')
par_gmsh_getdp.write('nm = %3.15e;\n'%(nm))
par_gmsh_getdp.write('epsilon0 = %3.15e;\n'%(epsilon0))
par_gmsh_getdp.write('mu0 = %3.15e;\n'%(mu0))
par_gmsh_getdp.write('cel = %3.15e;\n'%(cel))
par_gmsh_getdp.write('lambda_vp = %3.15e;\n'%(lambda_vp))
par_gmsh_getdp.write('lambda_m = %3.15e;\n'%(lambda_m))
par_gmsh_getdp.write('paramaille_air = %3.15e;\n'%(paramaille_air))
par_gmsh_getdp.write('paramaille_pml = %3.15e;\n'%(paramaille_pml))
par_gmsh_getdp.write('paramaille_mir = %3.15e;\n'%(paramaille_mir))
par_gmsh_getdp.write('eig_target = %3.15e;\n'%(eig_target))
par_gmsh_getdp.write('neig = %d ;\n'%(neig))
par_gmsh_getdp.close()
### define strings to be "oscommanded"
str_getdp1 = str_getdp_path+'getdp '+str_pro_filename+' -pre all -msh '+str_mesh_filename+' -cal -pos postop_eigenvectors_full -v 50 -bin'
str_slepc_options = ' -eps_monitor -eps_view \
-eps_type krylovschur \
-st_ksp_type preonly \
-st_pc_type lu \
-st_pc_factor_mat_solver_package mumps \
-eps_max_it 100 \
-eps_mpd 400 \
-eps_target '+ "%3.5e" % eig_target +' \
-eps_target_real \
-eps_nev '+"%d" % (neig)
# \
#-v 0'
# ### calling gmsh : mesh & save mesh!
os.system(str_gmsh_path+'gmsh geometry_haroche_realistic.geo -3 -o '+str_mesh_filename)
# ### calling getdp : preprocess and solve, using mpi if required
os.system(str_getdp1+str_slepc_options)
vp_real = np.loadtxt('./EigenValuesReal.txt',usecols=[5])
vp_imag = np.loadtxt('./EigenValuesImag.txt',usecols=[5])
print 'EigenFreq (GHz)', 2e9*np.pi*vp_real.transpose()
print 'vp imag', vp_imag.transpose()
os.system(str_gmsh_path+'gmsh haroche_postplot.geo &')
# os.system(str_gmsh_path+'gmsh eigenVectors_CompX_faces.pos geometry_haroche_realistic.geo &')
# omega2=(vp_real+1j*vp_imag)**2
# pl.plot(np.sqrt(-omega2.real).sort())
# print np.sort(2*pi*cel/vp_real*1e9)
# pl.figure()
# pl.plot(vp_real,vp_imag,'o')
# pl.figure()
# pl.plot(-vp_real,-vp_imag,'o')
# pl.show()
# omega_real_sort = np.sort(-omega2.real)
# pl.figure()
# pl.plot(omega_real_sort,'o')
# pl.figure()
# pl.plot(omega2.imag,'-')
# pl.show()
# vp_real_adj = np.loadtxt('./EigenValuesReal_adj.txt',usecols=[5])
# vp_imag_adj = np.loadtxt('./EigenValuesImag_adj.txt',usecols=[5])
tests/haroche/cavity_haroche_sf.py deleted 100644 → 0
+ 0
− 75
View file @ a24f61c9
import os
import sys
import time
import numpy
## Physical / Math constants ##
pi = numpy.pi
nm = 1.e-9
epsilon0 = 8.854187817e-3 * nm
mu0 = 400. * pi * nm
cel = 1.0 / (numpy.sqrt(epsilon0 * mu0))
## User data ##
if(len(sys.argv) != 6):
raise ValueError('Bad argument: '
'cavity_haroche_sf meshName symmetry femOrder nEig nProc')
meshName = str(sys.argv[1])
symmetry = int(sys.argv[2])
femOrder = int(sys.argv[3])
nEig = int(sys.argv[4])
nProc = int(sys.argv[5])
## My Data ##
tol = 1e-6
maxIt = 100
## Eigen Problem Shift ##
freqTarget = 51.099e9
lambdaVp = cel / freqTarget
eigTarget = (2. * pi * cel / lambdaVp)**2
## Start ##
start = time.time()
print '## Haroche'
print ' -- Process : %d' %nProc
print ' -- Mesh : %s' %meshName
print ' -- FEM order : %d' %femOrder
print ' -- N. Eigen : %d' %nEig
print ' -- Shift : %e' %eigTarget
print ' -- Symmetry : %e' %symmetry
print ' -- Tolerance : %e' %tol
print ' -- Max Iteration : %d' %maxIt
print ' '
## Calling small_fem ##
print '## Simulating'
os.system('mpirun -np %d' %nProc + \
' har' + \
' -msh %s' %meshName + \
' -o %d' %femOrder + \
' -n %d' %nEig + \
' -shift %e' %eigTarget + \
' -sym %d' %symmetry + \
' -tol %e' %tol + \
' -maxit %d' %maxIt + \
' -solver -eps_monitor -eps_view')
## Renaming Output ##
for i in range(nProc):
os.rename('harocheModes_proc%d.msh' %i + \
'post_proc%d_' %i + \
'femorder_%d_' %femOrder + \
'sym_%d_' %symmetry + \
meshName)
os.rename('harocheValues.txt', \
'eig_femorder_%d_' %femOrder + \
'sym_%d_' %symmetry + \
meshName + '.txt')
## Done ##
stop = time.time()
print '## Done in %e' %(stop - start) + ' [s]'
tests/haroche/geometry_haroche.py deleted 100644 → 0
+ 0
− 223
View file @ a24f61c9
#!/usr/bin/env python
from gmshpy import *
import numpy
import sys
import os
## Pyhsics ##
#############
## Scaling
mm = 1.e-3
nm = 1.e-9
## Speed of light
epsilon0 = 8.854187817e-3 * nm
mu0 = 400. * numpy.pi * nm
c = 1.0 / (numpy.sqrt(epsilon0 * mu0))
## Haroche Frequency, Wavelength & Wavenumber
freq_haroche = 51.099e9
lambda_haroche = c / freq_haroche
k_haroche = freq_haroche * 2 * numpy.pi / c
## Geomtrical Parameters ##
###########################
## Mirror
r = 39.4 * mm
R = 40.6 * mm
L_cav = 27.57 * mm
deltaz = 5. * mm
radius_mirror = 25. * mm
thick_mirror_atcenter = 1.415 * mm
dist = 1.00
dist2PML_xy = dist * lambda_haroche
dist2PML_z = dist * lambda_haroche
## Box
box_x = radius_mirror + dist2PML_xy
box_y = radius_mirror + dist2PML_xy
box_z = L_cav / 2. + thick_mirror_atcenter + dist2PML_z
## PML
thick = 1.00
pml_x = thick * lambda_haroche
pml_y = thick * lambda_haroche
pml_z = thick * lambda_haroche
## Dump PML data on disk
pml = open('pml.dat' ,'w')
pml.write(str( pml_x) + '\n')
pml.write(str( pml_y) + '\n')
pml.write(str( pml_z) + '\n')
pml.write(str( box_x) + '\n')
pml.write(str( box_y) + '\n')
pml.write(str( box_z) + '\n')
pml.write(str(k_haroche) + '\n')
os.fsync(pml)
pml.close()
## Mesh Parameters ##
#####################
if(len(sys.argv) != 5):
raise ValueError('Bad argument: '
'geometry_haroche air pml mir order')
refinement_air = float(sys.argv[1])
refinement_pml = float(sys.argv[2])
refinement_mir = float(sys.argv[3])
order = int(sys.argv[4])
paramaille_air = lambda_haroche / refinement_air
paramaille_pml = lambda_haroche / refinement_pml
paramaille_mir = lambda_haroche / refinement_mir
## Options ##
#############
# 2D mesh algorithm (1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
# 3D mesh algorithm (1=Delaunay, 4=Frontal, 5=Frontal Delaunay, 6=Frontal Hex, 7=MMG3D)
GmshSetOption('Mesh' , 'Algorithm' , 1.0)
GmshSetOption('Mesh' , 'Algorithm3D' , 1.0)
GmshSetOption('Mesh' , 'HighOrderOptimize', 1.0)
GmshSetOption('Mesh' , 'ElementOrder' , float(order))
GmshSetOption('General', 'Verbosity' , 4.0)
## Geomtrical Model ##
######################
myModel1 = GModel()
myModel2 = GModel()
myModel3 = GModel()
myModel4 = GModel()
myModel1.addCylinder([0, 0, -R + L_cav / 2.], \
[0, 0, L_cav / 2. + thick_mirror_atcenter], \
radius_mirror)
## Torus - almost spherical - centered in 0,0,0
myModel2.addTorus([0, 0, -R + L_cav / 2.], [0, 1, -R + L_cav / 2.], R-r, r)
myModel3.addSphere(0, 0, -R + L_cav / 2., r+r / 1000)
myModel4.addBlock([0, 0, 0], [box_x, box_y, box_z])
myModel1.computeBooleanDifference(myModel2);
myModel1.computeBooleanDifference(myModel3);
myModel4.computeBooleanDifference(myModel1);
## Extrusion of boundary surface to get PML region
myModel4.extrude(myModel4.getFaceByTag(9) , [0, 0, 0], [pml_x, 0, 0])
myModel4.extrude(myModel4.getFaceByTag(12), [0, 0, 0], [0, 0, pml_z])
myModel4.extrude(myModel4.getFaceByTag(8) , [0, 0, 0], [0, 0, pml_z])
myModel4.extrude(myModel4.getFaceByTag(23), [0, 0, 0], [0, pml_y, 0])
myModel4.extrude(myModel4.getFaceByTag(14), [0, 0, 0], [0, pml_y, 0])
myModel4.extrude(myModel4.getFaceByTag(25), [0, 0, 0], [0, pml_y, 0])
myModel4.extrude(myModel4.getFaceByTag(19), [0, 0, 0], [0, pml_y, 0])
myModel4.occconnect()
## Mesh size at Vertices
for i in range(1, myModel4.getNumVertices()):
myModel4.getVertexByTag(i).setPrescribedMeshSizeAtVertex(paramaille_pml)
myModel4.getVertexByTag(28).setPrescribedMeshSizeAtVertex(paramaille_mir)
myModel4.getVertexByTag(29).setPrescribedMeshSizeAtVertex(paramaille_mir)
myModel4.getVertexByTag(33).setPrescribedMeshSizeAtVertex(paramaille_mir)
myModel4.getVertexByTag(1).setPrescribedMeshSizeAtVertex(paramaille_air)
myModel4.getVertexByTag(3).setPrescribedMeshSizeAtVertex(paramaille_air)
myModel4.getVertexByTag(5).setPrescribedMeshSizeAtVertex(paramaille_air)
myModel4.getVertexByTag(7).setPrescribedMeshSizeAtVertex(paramaille_air)
myModel4.getVertexByTag(9).setPrescribedMeshSizeAtVertex(paramaille_air)
myModel4.getVertexByTag(11).setPrescribedMeshSizeAtVertex(paramaille_air)
myModel4.getVertexByTag(18).setPrescribedMeshSizeAtVertex(paramaille_air)
myModel4.getVertexByTag(27).setPrescribedMeshSizeAtVertex(paramaille_air)
myModel4.getVertexByTag(30).setPrescribedMeshSizeAtVertex(paramaille_air)
myModel4.getVertexByTag(31).setPrescribedMeshSizeAtVertex(paramaille_air)
myModel4.getVertexByTag(32).setPrescribedMeshSizeAtVertex(paramaille_air)
## Physical Entites ##
######################
myModel4.getRegionByTag(8).addPhysicalEntity(138) # Air
myModel4.getRegionByTag(1).addPhysicalEntity(139) # PML X
myModel4.getRegionByTag(4).addPhysicalEntity(141) # PML Y
myModel4.getRegionByTag(2).addPhysicalEntity(142) # PML Z
myModel4.getRegionByTag(6).addPhysicalEntity(140) # PML XY
myModel4.getRegionByTag(3).addPhysicalEntity(144) # PML XZ
myModel4.getRegionByTag(5).addPhysicalEntity(145) # PML YZ
myModel4.getRegionByTag(7).addPhysicalEntity(143) # PML XYZ
myModel4.getFaceByTag( 1).addPhysicalEntity(146) # XOZ
myModel4.getFaceByTag(10).addPhysicalEntity(146) # XOZ
myModel4.getFaceByTag(14).addPhysicalEntity(146) # XOZ
myModel4.getFaceByTag(35).addPhysicalEntity(146) # XOZ
myModel4.getFaceByTag( 7).addPhysicalEntity(147) # YOZ
myModel4.getFaceByTag(20).addPhysicalEntity(147) # YOZ
myModel4.getFaceByTag(24).addPhysicalEntity(147) # YOZ
myModel4.getFaceByTag(34).addPhysicalEntity(147) # YOZ
myModel4.getFaceByTag( 4).addPhysicalEntity(149) # XOY
myModel4.getFaceByTag(17).addPhysicalEntity(149) # XOY
myModel4.getFaceByTag(28).addPhysicalEntity(149) # XOY
myModel4.getFaceByTag(39).addPhysicalEntity(149) # XOY
myModel4.getFaceByTag( 6).addPhysicalEntity(150) # Outer PML
myModel4.getFaceByTag(12).addPhysicalEntity(150) # Outer PML
myModel4.getFaceByTag(15).addPhysicalEntity(150) # Outer PML
myModel4.getFaceByTag(16).addPhysicalEntity(150) # Outer PML
myModel4.getFaceByTag(22).addPhysicalEntity(150) # Outer PML
myModel4.getFaceByTag(25).addPhysicalEntity(150) # Outer PML
myModel4.getFaceByTag(26).addPhysicalEntity(150) # Outer PML
myModel4.getFaceByTag(29).addPhysicalEntity(150) # Outer PML
myModel4.getFaceByTag(30).addPhysicalEntity(150) # Outer PML
myModel4.getFaceByTag(31).addPhysicalEntity(150) # Outer PML
myModel4.getFaceByTag(32).addPhysicalEntity(150) # Outer PML
myModel4.getFaceByTag(33).addPhysicalEntity(150) # Outer PML
myModel4.getFaceByTag(36).addPhysicalEntity(148) # Mirror
myModel4.getFaceByTag(37).addPhysicalEntity(151) # Frame
myModel4.getFaceByTag(38).addPhysicalEntity(151) # Frame
myModel4.getVertexByTag(1).addPhysicalEntity(1000000) # Dummy point for GetDP
## Mesh & Save ##
#################
brepName = 'haroche.brep'
meshName = 'haroche' + \
'_geoorder_%d' %order + \
'_air_%d' %refinement_air + \
'_pml_%d' %refinement_pml + \
'_mir_%d' %refinement_mir + '.msh'
myModel4.mesh(3)
myModel4.save(meshName)
myModel4.save(brepName)
## Display ##
#############
print('Data used: ')
print(' ** Air : ' + str(refinement_air))
print(' ** PML : ' + str(refinement_pml))
print(' ** Mirror : ' + str(refinement_mir))
print(' ** Order : ' + str(order))
## Dump on disk
log = open('mesh.log' ,'w')
log.write('Data used: \n')
log.write(' ** Air : ' + str(refinement_air) + '\n')
log.write(' ** PML : ' + str(refinement_pml) + '\n')
log.write(' ** Mirror : ' + str(refinement_mir) + '\n')
log.write(' ** Order : ' + str(order) + '\n')
os.fsync(log)
log.close()
tests/haroche/geometry_haroche_full.py deleted 100644 → 0
+ 0
− 69
View file @ a24f61c9
#!/usr/bin/env python
from gmshpy import *
import numpy
import sys
import os
## Pyhsics ##
#############
## Scaling
mm = 1.e-3
nm = 1.e-9
## Speed of light
epsilon0 = 8.854187817e-3 * nm
mu0 = 400. * numpy.pi * nm
c = 1.0 / (numpy.sqrt(epsilon0 * mu0))
## Haroche Frequency, Wavelength & Wavenumber
freq_haroche = 51.099e9
lambda_haroche = c / freq_haroche
k_haroche = freq_haroche * 2 * numpy.pi / c
## Geomtrical Parameters ##
###########################
r = 39.4 * mm
R = 40.6 * mm
L_cav = 27.57 * mm
deltaz = 5. * mm
radius_mirror = 25. * mm
thick_mirror_atcenter = 1.415 * mm
## Geomtrical Model ##
######################
myModel11 = GModel()
myModel21 = GModel()
myModel31 = GModel()
myModel12 = GModel()
myModel22 = GModel()
myModel32 = GModel()
myModel11.addCylinder([0, 0, -R + L_cav / 2.], \
[0, 0, +L_cav / 2. + thick_mirror_atcenter], \
radius_mirror)
myModel12.addCylinder([0, 0, +R - L_cav / 2.], \
[0, 0, -L_cav / 2. - thick_mirror_atcenter], \
radius_mirror)
## Torus - almost spherical - centered in 0,0,0
myModel21.addTorus([0, 0, -R + L_cav / 2.], [0, 1, -R + L_cav / 2.], R-r, r)
myModel22.addTorus([0, 0, +R - L_cav / 2.], [0, 1, +R - L_cav / 2.], R-r, r)
myModel31.addSphere(0, 0, -R + L_cav / 2., r+r / 1000)
myModel32.addSphere(0, 0, +R - L_cav / 2., r+r / 1000)
myModel11.computeBooleanDifference(myModel21);
myModel11.computeBooleanDifference(myModel31);
myModel12.computeBooleanDifference(myModel22);
myModel12.computeBooleanDifference(myModel32);
myModel11.computeBooleanUnion(myModel12);
## Save ##
##########
brepName = 'haroche.brep'
myModel11.save(brepName)
tests/haroche/haroche_postplot.geo deleted 100644 → 0
+ 0
− 28
View file @ a24f61c9
Merge "eigenVectors.pos" ;
View[0].Visible = 0;
Plugin(MathEval).View = 0;
Plugin(MathEval).Expression0 = 'v0' ;
Plugin(MathEval).Run ;
View[1].Visible = 0;
Plugin(CutPlane).A = 1;
Plugin(CutPlane).B = 0;
Plugin(CutPlane).C = 0;
Plugin(CutPlane).D = 0;
Plugin(CutPlane).View = 1 ;
Plugin(CutPlane).Run ;
Plugin(CutPlane).A = 0;
Plugin(CutPlane).B = 1;
Plugin(CutPlane).C = 0;
Plugin(CutPlane).D = 0;
Plugin(CutPlane).View = 1 ;
Plugin(CutPlane).Run ;
Plugin(CutPlane).A = 0;
Plugin(CutPlane).B = 0;
Plugin(CutPlane).C = 1;
Plugin(CutPlane).D = 0;
Plugin(CutPlane).View = 1 ;
Plugin(CutPlane).Run ;
tests/haroche/partition.py deleted 100644 → 0
+ 0
− 28
View file @ a24f61c9
#!/usr/bin/env python
from gmshpy import *
import sys
import os
## Get Data ##
##############
if(len(sys.argv) != 3):
raise ValueError('Bad argument: '
'partition filename number_of_partitions')
name = str(sys.argv[1])
part = int(sys.argv[2])
## Read Mesh ##
###############
model = GModel()
model.readMSH(name)
## Partition ##
###############
partitionOpt = meshPartitionOptions()
partitionOpt.setNumOfPartitions(part)
PartitionMesh(model, partitionOpt)
## Save ##
##########
model.save(os.path.splitext(name)[0] + '_part_' + str(part) + '.msh')
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment