Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

added print_chunk_costs #983

Merged
merged 6 commits into from
Aug 20, 2019
Merged

added print_chunk_costs #983

merged 6 commits into from
Aug 20, 2019

Conversation

stevengj
Copy link
Collaborator

@stevengj stevengj commented Aug 7, 2019

Add a routine to print the estimated chunk costs. Run it with sim.structure.print_chunk_costs() sim.structure.print_estimated_costs().

Update: now prints as part of choose_chunkdivision if !quiet

@stevengj
Copy link
Collaborator Author

stevengj commented Aug 7, 2019

Might not be quite what you want — we probably want to sum the chunks on each process.

Updated: PR now prints cost per process.

for (int i = 0; i < count_processors(); i++)
costs[i] = 0;
for (int i = 0; i < num_chunks; i++)
costs[chunks[i]->n_proc()] += chunk_cost(i);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For some reason, this line seems to be producing a segmentation fault at runtime even though make for this branch/PR finishes without errors.

[simpetus:16699] [ 0] /lib/x86_64-linux-gnu/libpthread.so.0(+0x12890)[0x7f85427c3890]
[simpetus:16699] [ 1] /usr/local/lib/libctlgeom.so.7(+0x3bc4)[0x7f8540759bc4]
[simpetus:16699] [ 2] /usr/local/lib/libctlgeom.so.7(geom_fix_object_ptr+0x479)[0x7f8540761549]
[simpetus:16699] [ 3] /usr/local/lib/libctlgeom.so.7(geom_get_bounding_box+0x34)[0x7f8540761da4]
[simpetus:16699] [ 4] /usr/local/lib/libctlgeom.so.7(overlap_with_object+0x123)[0x7f8540763103]
[simpetus:16699] [ 5] /home/oskooi/install/meep6/src/.libs/libmeep.so.16(_ZN9meep_geom14fragment_stats13compute_statsEv+0x9a)[0x7f8540c156aa]
[simpetus:16699] [ 6] /home/oskooi/install/meep6/src/.libs/libmeep.so.16(_ZN9meep_geom14fragment_stats7computeEv+0x9)[0x7f8540c16269]
[simpetus:16699] [ 7] /home/oskooi/install/meep6/src/.libs/libmeep.so.16(_ZNK4meep11grid_volume8get_costEv+0x55)[0x7f8540be5555]
[simpetus:16699] [ 8] /home/oskooi/install/meep6/src/.libs/libmeep.so.16(_ZNK4meep9structure21print_estimated_costsEv+0xc3)[0x7f8540bab773]
[simpetus:16699] [ 9] /home/oskooi/install/meep6/python/meep/_meep.so(+0x62533)[0x7f8540ebc533]
[simpetus:16699] [10] corrupted double-linked list

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe add a master_printf("n_proc[%d] = %d\n", i, chunks[i]->n_proc()); line above this line to see what the value of n_proc is?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

n_proc() is not the issue (i.e., output is as expected), rather it's the function call chunks[i]->gv.get_cost() inside chunk_cost(int i) which is causing the error. Could this have something to do with get_cost() being a const function?

@oskooi
Copy link
Collaborator

oskooi commented Aug 16, 2019

empty cell

import meep as mp

cell_size = mp.Vector3(10,10,10)

sim = mp.Simulation(resolution=10,
                    cell_size=(10,10,10),
                    split_chunks_evenly=False.
                    collect_stats=True)

sim.init_sim()
sim.structure.print_estimated_costs()

output

$ mpirun -n 5 python empty.py
Using MPI version 3.1, 5 processes
-----------
Initializing structure...
Splitting into 5 chunks by cost
time for choose_chunkdivision = 0.000195965 s
Working in 3D dimensions.
Computational cell is 10 x 10 x 10 with resolution 10
time for set_epsilon = 0.506503 s
-----------
estimated costs per process: 0, 0, 0, 0, 0
estimated cost mean = 0, stddev = 0

Elapsed run time = 0.5075 s

@stevengj
Copy link
Collaborator Author

If I call compute_fragment_stats explicitly then it works (it gives a nonzero cost):

import meep as mp

cell_size = mp.Vector3(10,10,10)

sim = mp.Simulation(resolution=10,
                    cell_size=(10,10,10),
                    split_chunks_evenly=False)

sim.init_sim()
sim._compute_fragment_stats(sim.structure.user_volume)
sim.structure.print_estimated_costs()

@ChristopherHogan
Copy link
Contributor

ChristopherHogan commented Aug 17, 2019

When python passes a Medium to a C++ function, a typemap converts it to a material_type via pymaterial_to_material, and when that function exits, the C++ material_type is destroyed. The issue in this PR is that computing fragment statistics relies on the global libctl variable default_material, but by the time print_estimated_costs is called, the material_type it points to has been deleted.

Two solutions come to mind:

  1. Print the costs within choose_chunk_division if a flag is present.
  2. Write a python wrapper around structure::print_estimated_costs that passes sim.default_material and sim.geometry to a function like fragment_stats::init_libctl.

@stevengj stevengj merged commit 6cb30a0 into master Aug 20, 2019
@stevengj stevengj deleted the print_chunk_costs branch August 20, 2019 01:28
@oskooi
Copy link
Collaborator

oskooi commented Aug 21, 2019

In the following example involving a grating outcoupler, the standard deviation value of the chunk costs is exactly 0 over the entire range of chunk counts when the individual chunk costs are not equal. The standard deviation should be a non-zero value. (Note that the sum of the costs is a constant, ~47186.6, which is correct.)

$ for i in `seq 11`; do echo "n = "$i; mpirun -n $i python grating.py -res 30 -N 6 -nfreq 30 |grep estimated; done
n = 1
estimated costs per process: 47186.6
estimated cost mean = 47186.6, stddev = 0
n = 2
estimated costs per process: 23554, 23632.7
estimated cost mean = 23593.3, stddev = 0
n = 3
estimated costs per process: 15727.7, 15697.3, 15761.7
estimated cost mean = 15728.9, stddev = 0
n = 4
estimated costs per process: 11753.1, 11800.9, 11800.8, 11831.7
estimated cost mean = 11796.6, stddev = 0
n = 5
estimated costs per process: 9378.84, 9432.75, 9445.12, 9438.18, 9491.9
estimated cost mean = 9437.36, stddev = 0
n = 6
estimated costs per process: 7847.11, 7880.62, 7809.78, 7887.5, 7842.1, 7919.65
estimated cost mean = 7864.46, stddev = 0
n = 7
estimated costs per process: 7847.11, 7880.62, 7809.78, 7887.5, 7842.1, 7919.65, 0
estimated cost mean = 6740.97, stddev = 0
n = 8
estimated costs per process: 5874.55, 5878.61, 5888.39, 5912.54, 5893.24, 5909.2, 5903.4, 5928.3
estimated cost mean = 5898.53, stddev = 0
n = 9
estimated costs per process: 5227.39, 5248.88, 5251.49, 5203.48, 5239.18, 5254.66, 5225.56, 5258.03, 5278.2
estimated cost mean = 5242.99, stddev = 0
n = 10
estimated costs per process: 4678.67, 4700.19, 4707.26, 4725.51, 4696.16, 4748.99, 4716.49, 4721.66, 4741, 4750.93
estimated cost mean = 4718.69, stddev = 0
n = 11
estimated costs per process: 4678.67, 4700.19, 4707.26, 4725.51, 4696.16, 4748.99, 4716.49, 4721.66, 4741, 4750.93, 0
estimated cost mean = 4289.71, stddev = 0

grating.py

import meep as mp
from meep.materials import Si, SiO2_aniso
import math
import argparse

def main(args):
    h = args.hh
    w = args.w
    a = args.a
    d = args.d
    N = args.N
    N = N+1

    wvl_cen = 1.55
    dair_sub = 0.5

    dpml = wvl_cen
    boundary_layers = [mp.PML(dpml)]

    sxy = 2*(N+dpml)
    sz = dpml+dair_sub+h+dair_sub+dpml
    cell_size = mp.Vector3(sxy,sxy,sz)

    geometry = []

    # rings of Bragg grating                                                                                                                                                                 
    for n in range(N,0,-1):
        geometry.append(mp.Cylinder(material=Si,
                                    center=mp.Vector3(0,0,0),
                                    radius=n*a,
                                    height=h))
        geometry.append(mp.Cylinder(material=mp.air,
                                    center=mp.Vector3(0,0,0),
                                    radius=n*a-d,
                                    height=h))

    # remove left half of Bragg grating rings to form semi circle                                                                                                                            
    geometry.append(mp.Block(material=mp.air,
                             center=mp.Vector3(-0.5*N*a,0,0),
                             size=mp.Vector3(N*a,2*N*a,h)))
    geometry.append(mp.Cylinder(material=Si,
                                center=mp.Vector3(0,0,0),
                                radius=a-d,
                                height=h))

    # angle sides of Bragg grating                                                                                                                                                           

    # rotation angle of sides relative to Y axis (degrees)                                                                                                                                   
    rot_theta = -math.radians(args.rot_theta)

    pvec = mp.Vector3(0,0.5*w,0)
    cvec = mp.Vector3(-0.5*N*a,0.5*N*a+0.5*w,0)
    rvec = cvec-pvec
    rrvec = rvec.rotate(mp.Vector3(0,0,1), rot_theta)

    geometry.append(mp.Block(material=mp.air,
                             center=pvec+rrvec,
                             size=mp.Vector3(N*a,N*a,h),
                             e1=mp.Vector3(1,0,0).rotate(mp.Vector3(0,0,1),rot_theta),
                             e2=mp.Vector3(0,1,0).rotate(mp.Vector3(0,0,1),rot_theta),
                             e3=mp.Vector3(0,0,1)))

    pvec = mp.Vector3(0,-0.5*w,0)
    cvec = mp.Vector3(-0.5*N*a,-(0.5*N*a+0.5*w),0)
    rvec = cvec-pvec
    rrvec = rvec.rotate(mp.Vector3(0,0,1),-rot_theta)

    geometry.append(mp.Block(material=mp.air,
                             center=pvec+rrvec,
                             size=mp.Vector3(N*a,N*a,h),
                             e1=mp.Vector3(1,0,0).rotate(mp.Vector3(0,0,1),-rot_theta),
                             e2=mp.Vector3(0,1,0).rotate(mp.Vector3(0,0,1),-rot_theta),
                             e3=mp.Vector3(0,0,1)))

    # input waveguide                                                                                                                                                                        
    geometry.append(mp.Block(material=mp.air,
                             center=mp.Vector3(-0.25*sxy,0.5*w+0.5*a,0),
                             size=mp.Vector3(0.5*sxy,a,h)))
    geometry.append(mp.Block(material=mp.air,
                             center=mp.Vector3(-0.25*sxy,-(0.5*w+0.5*a),0),
                             size=mp.Vector3(0.5*sxy,a,h)))
    geometry.append(mp.Block(material=Si,
                             center=mp.Vector3(-0.25*sxy,0,0),
                             size=mp.Vector3(0.5*sxy,w,h)))

    # substrate                                                                                                                                                                              
    geometry.append(mp.Block(material=SiO2_aniso,
                             center=mp.Vector3(0,0,-0.5*sz+0.25*(sz-h)),
                             size=mp.Vector3(mp.inf,mp.inf,0.5*(sz-h))))

    # mode frequency                                                                                                                                                                         
    fcen = 1/wvl_cen

    sim = mp.Simulation(resolution=args.res,
                        cell_size=cell_size,
                        boundary_layers=boundary_layers,
                        geometry=geometry,
                        dimensions=3,
                        split_chunks_evenly=False,
                        eps_averaging=False)

    nearfield = sim.add_near2far(fcen, 0.2*fcen, args.nfreq, mp.Near2FarRegion(mp.Vector3(0,0,0.5*sz-dpml), size=mp.Vector3(sxy-2*dpml,sxy-2*dpml,0)))

    sim.init_sim()

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('-hh', type=float, default=0.22, help='wavelength height (default: 0.22 um)')
    parser.add_argument('-w', type=float, default=0.50, help='wavelength width (default: 0.50 um)')
    parser.add_argument('-a', type=float, default=1.0, help='Bragg grating periodicity/lattice parameter (default: 1.0 um)')
    parser.add_argument('-d', type=float, default=0.5, help='Bragg grating thickness (default: 0.5 um)')
    parser.add_argument('-N', type=int, default=5, help='number of grating periods')
    parser.add_argument('-rot_theta', type=float, default=20, help='rotation angle of sides relative to Y axis (default: 20 degrees)')
    parser.add_argument('-res', type=int, default=30, help='resolution (default: 30 pixels/um)')
    parser.add_argument('-nfreq', type=int, default=50, help='number of frequency bins (default: 50)')
    args = parser.parse_args()
    main(args)

stevengj added a commit that referenced this pull request Aug 23, 2019
@stevengj
Copy link
Collaborator Author

Should be fixed by 7f7e85d

@stevengj
Copy link
Collaborator Author

Now that #994 is merged you will need to call mp.verbosity(2) in order to see this output.

bencbartlett pushed a commit to bencbartlett/meep that referenced this pull request Sep 9, 2021
* added print_chunk_costs

* print estimated cost per process, not per chunk

* initialize to zero

* tweak

* print in choose_chunkdivision

* only print cost estimates for split by cost
bencbartlett pushed a commit to bencbartlett/meep that referenced this pull request Sep 9, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants