sequential_distribute_dump fails

[ninnghazad]

sequential_distribute_dump fails for large domains in cPickle with

File "/usr/local/lib/python2.7/dist-packages/anuga/parallel/sequential_distribute.py", line 246, in sequential_distribute_dump
    cPickle.dump( tostore, f, protocol=cPickle.HIGHEST_PROTOCOL)
SystemError: error return without exception set

Probably related:
numpy/numpy#2396
https://bugs.python.org/issue11564
https://stackoverflow.com/questions/34091717/avoid-cpickle-error-alternatives-to-cpickle-on-python-2-7

As using python3 isn't easily done it seems, i am unable to save large domains. (roughly >100,000,000 tris) - which is rather annoying when trying multiple variants on same large mesh.

My only advice to other users at this moment would be to not use sequential_distribute_dump to save the combined domain and the parts, but only the partial (per-proc) domains. these might just be small enough to work.

As to how to circumvent the error without upgrading anuga to python3, using numpy.save()/numpy.savez() is the most often suggested solution. (not sure if anuga uses numpy arrays for fields internally)

[stoiver]

@ninnghazad it wouldn't be too hard to store the information in a netcdf file. Most of the data is either numpy arrays or strings. Maybe boundary_map is a dictionary.

See line 200 of sequential_distribute.py in anuga/parallel, where the data to be stored is accumulated into an object called tostore, which we then pickle it in the code around line 242 and load around line 272.

Could make tostore into a class that collects the info and have alternative ways to store and load, ie pickle or netcdf or a combination of numpy save and pickle for the strings.

[ninnghazad]

I am currently trying a combination of numpy.save to seperate files and cPickle for the rest. Have to do some bigger tests to see if that is enough for me. However, other parts of the subdivide/dump code seem to not like big meshes either, like pmesh_divide_metis_with_map:
sequential_distribute: Subdivide mesh Error! ***Memory allocation failed for AllocateWorkSpace: edegrees. Requested size: -1415511280 bytes[HOSTNAME:69273] *** Process received signal ***
That is with 512G Ram and a lot of disk to spare, so i guess it's some 32bit-limit.

Has anybody successfully done simulations using large (as in > 100Mio tris) meshes before?

EDIT: Just noticed there also is #15

[stoiver]

@ninnghazad what system are you working on? I am sure we have had a user run with 50M triangles, but I thought he had also gone higher.

[ninnghazad]

@stoiver it's some custom numbercruncher, uses two Xeon E5 18c/36t, 512G mem, infiniband and whatnot. it's running debian on a 4.9 kernel and it's a headless system.

my goal is to get to above 500Mio tris to simulate strong rainfall over a medium sized city.

[stoiver]

@ninnghazad sounds like an interesting project. It would be great to insure anuga can run on such large problems. You have a large machine, but my guess is that you will need to run such large problems on the nci machine.

But I am happy to help. An obvious first step would be to change over to metis 5 as suggested in issue #15

[ninnghazad]

@stoiver it surely is. ANUGA could be a rather big thing around here, the flexibility as well as being public open source.

I'll see how goes, but for a start i want to test the whole thing using openmpi without the added hassle of networking and go from there. I assume a simulation working in local parallel is easier to migrate to other nodes than one where i have no idea if the problem is the simulation or the setup.

Anyway, i appreciate you taking the time to talk about this. The original reason for this issue i seem to have circumvented by using numpy's save(), but i'll test more before considering that solved and possibly uploading the rather naive changes i made.

About #15, i actually am trying to shoehorn latest pymetis into distribute_mesh.py. I am still testing, and the way i did it might be too simple, but first results look like it might work. Had to convert domain.neighbours from ndarray to list-copy though, and that might be a rather large amount of not really needed memory-usage.

Back to testing.

[stoiver]

@ninnghazad, I had a look at pymetis, it is probably a good way to go, but an alternative would be to use ctypes with a standard installation of metis 5.1.

ctypes was not available back when we started the parallel version, back in 2005.

[ninnghazad]

some diff

diff --git a/anuga/parallel/sequential_distribute.py b/anuga/parallel/sequential_distribute.py
index 66939130..35184224 100644
--- a/anuga/parallel/sequential_distribute.py
+++ b/anuga/parallel/sequential_distribute.py
@@ -5,7 +5,6 @@
 
 import numpy as num
 
-
 from anuga import Domain
 
 from anuga.parallel.distribute_mesh  import send_submesh
@@ -217,7 +216,7 @@ def sequential_distribute_dump(domain, numprocs=1, verbose=False, partition_dir=
     """
 
     from os.path import join
-    
+
     partition = Sequential_distribute(domain, verbose, debug, parameters)
 
     partition.distribute(numprocs)
@@ -234,17 +233,25 @@ def sequential_distribute_dump(domain, numprocs=1, verbose=False, partition_dir=
             if exception.errno != errno.EEXIST:
                 raise
 
-    
+    import cPickle
     for p in range(0, numprocs):
 
         tostore = partition.extract_submesh(p) 
 
-        import cPickle
         pickle_name = partition.domain_name + '_P%g_%g.pickle'% (numprocs,p)
         pickle_name = join(partition_dir,pickle_name)
         f = file(pickle_name, 'wb')
-        cPickle.dump( tostore, f, protocol=cPickle.HIGHEST_PROTOCOL)
 
+       lst = list(tostore)
+
+       num.save(pickle_name+".np1",tostore[1],allow_pickle=False) # this append .npy to filename
+       lst[1] = pickle_name+".np1.npy"
+       num.save(pickle_name+".np2",tostore[2],allow_pickle=False)
+       lst[2] = pickle_name+".np2.npy"
+
+       tostore = tuple(lst)
+
+        cPickle.dump( tostore, f, protocol=cPickle.HIGHEST_PROTOCOL)
     return
 
 
@@ -257,7 +264,7 @@ def sequential_distribute_load(filename = 'domain', partition_dir = '.', verbose
 
     pickle_name = filename+'_P%g_%g.pickle'% (numprocs,myid)
     pickle_name = join(partition_dir,pickle_name) 
-    
+
     return sequential_distribute_load_pickle_file(pickle_name, numprocs, verbose = verbose)
 
 
@@ -265,9 +272,9 @@ def sequential_distribute_load_pickle_file(pickle_name, np=1, verbose = False):
     """
     Open pickle files
     """
-    
-    import cPickle    
+
     f = file(pickle_name, 'rb')
+    import cPickle
 
     kwargs, points, vertices, boundary, quantities, boundary_map, \
                    domain_name, domain_dir, domain_store, domain_store_centroids, \
@@ -276,6 +283,9 @@ def sequential_distribute_load_pickle_file(pickle_name, np=1, verbose = False):
                    domain_quantities_to_be_stored, domain_smooth = cPickle.load(f)
     f.close()
 
+    points = num.load(points)
+    vertices = num.load(vertices)
+
     #---------------------------------------------------------------------------
     # Create domain (parallel if np>1)
     #---------------------------------------------------------------------------

@stoiver this works for me, but it's kinda meh. also didn't look into domain.quantities to see if there are more ndarrays. i would assume that the above would still be problematic if one would not as i do set quantities after the distributed mesh is loaded from cache (so quantities aren't all saved in the cache files, making them smaller). have not yet looked into the save-as-ncdf way you mentioned.

however, a 120mio triangle test worked (see #15). took ~210sec walltime per 1sec simtime on a wet start and with a bit of rain. merging swws at the end of the simulation isn't quite working yet, but sww_merge.py did the job, so it's probably just some usage-error.
as i somewhat expected the merged sww is kind of hard to handle. crayfish won't display it (but gives a proper "std::bad_alloc"), so i will use the to_raster stuff from one of the examples.

and some stats (1sec yields):

Time = 60.0000, delta t in [0.04204673, 0.04808688], steps=22 (207s)
Communication time 1330.96 seconds
Reduction Communication time 1986.94 seconds
Broadcast time 0.00 seconds
Total time 12425.78 seconds
Finished simulation on 70 processors ...

that's a lot of comm-overhead right there - suspect my openmpi might be faulty, but that's another issue.

[stoiver]

@ninnghazad great work. I'll checkout your #15 work. But a comment on your timings. I have found that there is not much advantage to using more than the number of physical cores when running in parallel. So I would suggest trying something like 36 processors and see what timings you get.

[stoiver]

@ninnghazad, by the way can you provide some more info on your project. It would be good to point to your projects using these large grids.

[ninnghazad]

@stoiver The project just started and will probably take at least a bunch of months to finish. There will be some kind of publication, but i can't yet say if all of the resulting data will be public. It entails multiple simulations for a whole city. I will however make sure ANUGA is prominently named in whatever publications or reports there will be, and get back to you with any material i am allowed to share once the project is further along the timeline. Apart from actual results i am of course happy to share the hows and tips i can give. But there still are ways to go and stuff to figure out before i would advice anybody to use my run.py and setup. Just ran into some problem somewhere above 120 mio triangles when doing:

domain = anuga.shallow_water.shallow_water_domain.Domain(p,v,b)

On the other hand converting the 120tri-results to raster (GeoTIFF) for viewing was slow but worked.

[stoiver]

@ninnghazad, I was just curious about your project. There is no expectation for you to make the data public. But updates to the code would be much appreciated!

[ninnghazad]

@stoiver the goal of the project is to create risk-assessment maps for Dortmund, a city of about 280km2 in western Germany, incorporating local pumping-stations, surface friction as well as different time/rain-sets using a resolution of at most 1m2.

[ninnghazad]

@stoiver by the way, updates to the code, do you prefer PRs or diffs?

[stoiver]

@ninnghazad our preferred method is PRs. Thanks for the info about your project.

[Girishchandra-Yendargaye]

@ninnghazad Hello Sir...Could you please tell How do you overcome the issues of Large mesh generation and Domain Creation which is sequential step ...We are also trying the same....

[ninnghazad]

@Girishchandra-Yendargaye
I am currently using rectangular meshes, to generate these:

• use a custom tool to generate .json input from large TIF files. creates triangle,point and boundary files.
• use a python script to read these into numpy and pass the arrays to anuga.shallow_water.shallow_water_domain.Domain(p,v,b)
• use machines with at least 512GB RAM

This method is not optimal and could easily optimized by for example using numpy's memmap format instead of json as intermediary.

Also remember that depending on what you are simulating, calculating on domains that large will take a lot of time. A lot of fast CPUs (with AVX512 preferably) and using properly optimized python/numpy helps. Be prepared to wait weeks for simulations to finish.

[Girishchandra-Yendargaye]

@Girishchandra-Yendargaye
I am currently using rectangular meshes, to generate these:

• use a custom tool to generate .json input from large TIF files. creates triangle,point and boundary files.
• use a python script to read these into numpy and pass the arrays to anuga.shallow_water.shallow_water_domain.Domain(p,v,b)
• use machines with at least 512GB RAM

This method is not optimal and could easily optimized by for example using numpy's memmap format instead of json as intermediary.

Also remember that depending on what you are simulating, calculating on domains that large will take a lot of time. A lot of fast CPUs (with AVX512 preferably) and using properly optimized python/numpy helps. Be prepared to wait weeks for simulations to finish.

Can you please provide some sample code... As I dint find such method in shallow water which accepts only triangle,points and boundary Also What about elevation setting in next step which is sequential again...

[ninnghazad]

To build the domain i generate (custom) json and use that like this:

                if verbose:
			print 'Parse boundaries ...'
		tmp = json.load(open('input/mesh_'+str(project.gebiet)+'.boundaries.json'))
		b = {}
		for i in tmp:
			b[(i[0],i[1])] = 'exterior'
		tmp = None
		gc.collect();
	
		if verbose:
			print 'Parse vertices ...'
		p = json.load(open('input/mesh_'+str(project.gebiet)+'.points.json'))
		p = numpy.array(p, numpy.float)
		print "point: ",p[0]
	
		if verbose:
			print 'Parse elements ...'
		v = json.load(open('input/mesh_'+str(project.gebiet)+'.triangles.json'))
		v = numpy.array(v, numpy.int)
	
		gc.collect();
		if verbose:
			print 'Creating domain ...'
		domain = anuga.shallow_water.shallow_water_domain.Domain(p,v,b)
	
		if verbose:
			print "Free some memory ..."
		p = None
		v = None
		b = None
		gc.collect();

I do the above in sequential, and have one big "master-domain", which i then split into however many i need in parallel. When parallel simulation is done, i merge the pieces again, in sequential cause of RAM.

To apply heights from GeoTIFF to domain i use something like this:

def demValue(x,y):
	global minx,miny,maxx,maxy,w,h,rangex,rangey
	r = []
	print "DEM local extent: ",domain.get_extent()

	for xx,yy in zip(x,y):
		px = int(((xx-minx)/rangex)*(w))
		px = min(px,w-1)
		py = int(h-((yy-miny)/rangey)*(h))
		py = min(py,h-1)
		if px < 0 or py < 0 or px >= w or py >= h:
			print "DEM OOA: ",xx,yy,px,py
			r.append(project.noDataValue)
			continue
		pixel = dem.ReadRaster(px , py, 1, 1, 1, 1)
		value = struct.unpack(demDataType, pixel)[0]
		r.append(value)
	return r


import osgeo.gdal
import osgeo.gdalconst
indataset = osgeo.gdal.Open( project.dem_filename, osgeo.gdal.GA_ReadOnly )
dem = indataset.GetRasterBand(1)
gt = indataset.GetGeoTransform()
if gt is None or gt == (0.0, 1.0, 0.0, 0.0, 0.0, 1.0):
	sys.stdout.write("rasterfile has no srs set");
	sys.exit(1)
w = dem.XSize
h = dem.YSize
demDataType = pt2fmt(dem.DataType)
minx = gt[0]
maxx = gt[0]+(gt[1] * w)
rangex = maxx - minx
miny = gt[3]+(gt[5] * h)
maxy = gt[3]
rangey = maxy - miny

if myid == 0:
	print "DEM Geotransform: ",minx,maxx,rangex,miny,maxy,rangey
domain.set_quantity('elevation',demValue);

So the elevation setting can be sequential - but i doesn't have to be, i have use it like that on about 70 cores in parallel. It's not optimized, but as the time it takes compared to the simulation itself is tiny, i didn't mind.

Please remember that these are just pieces of code i grabbed out of some of my projects - this is not official ANUGA code and i am not responsible if it burns down your datacenter.

Oh, and the JSON format is as simple as can be:

$cat input/mesh_0.boundaries.json | xxd | head
00000000: 5b5b 302c 325d 2c5b 302c 305d 2c5b 322c  [[0,2],[0,0],[2,
00000010: 325d 2c5b 342c 325d 2c5b 362c 325d 2c5b  2],[4,2],[6,2],[
00000020: 382c 325d 2c5b 392c 305d 2c5b 3130 2c32  8,2],[9,0],[10,2
00000030: 5d2c 5b31 302c 305d 2c5b 3132 2c32 5d2c  ],[10,0],[12,2],
00000040: 5b31 342c 325d 2c5b 3136 2c32 5d2c 5b33  [14,2],[16,2],[3
00000050: 302c 325d 2c5b 3332 2c32 5d2c 5b33 342c  0,2],[32,2],[34,
00000060: 325d 2c5b 3336 2c32 5d2c 5b33 382c 325d  2],[36,2],[38,2]
00000070: 2c5b 3430 2c32 5d2c 5b34 322c 325d 2c5b  ,[40,2],[42,2],[
00000080: 3434 2c32 5d2c 5b34 352c 305d 2c5b 3436  44,2],[45,0],[46
00000090: 2c32 5d2c 5b34 362c 305d 2c5b 3438 2c32  ,2],[46,0],[48,2
$ cat input/mesh_0.points.json | xxd | head
00000000: 5b5b 3337 3338 3938 2e35 3030 2c35 3637  [[373898.500,567
00000010: 3739 3434 2e35 3030 5d2c 5b33 3733 3839  7944.500],[37389
00000020: 392e 3530 302c 3536 3737 3934 342e 3530  9.500,5677944.50
00000030: 305d 2c5b 3337 3339 3030 2e35 3030 2c35  0],[373900.500,5
00000040: 3637 3739 3434 2e35 3030 5d2c 5b33 3733  677944.500],[373
00000050: 3930 312e 3530 302c 3536 3737 3934 342e  901.500,5677944.
00000060: 3530 305d 2c5b 3337 3339 3032 2e35 3030  500],[373902.500
00000070: 2c35 3637 3739 3434 2e35 3030 5d2c 5b33  ,5677944.500],[3
00000080: 3733 3930 332e 3530 302c 3536 3737 3934  73903.500,567794
00000090: 342e 3530 305d 2c5b 3337 3338 3934 2e35  4.500],[373894.5
$ cat input/mesh_0.triangles.json | xxd | head
00000000: 5b5b 312c 302c 3130 5d2c 5b31 302c 3131  [[1,0,10],[10,11
00000010: 2c31 5d2c 5b32 2c31 2c31 315d 2c5b 3131  ,1],[2,1,11],[11
00000020: 2c31 322c 325d 2c5b 332c 322c 3132 5d2c  ,12,2],[3,2,12],
00000030: 5b31 322c 3133 2c33 5d2c 5b34 2c33 2c31  [12,13,3],[4,3,1
00000040: 335d 2c5b 3133 2c31 342c 345d 2c5b 352c  3],[13,14,4],[5,
00000050: 342c 3134 5d2c 5b31 342c 3135 2c35 5d2c  4,14],[14,15,5],
00000060: 5b37 2c36 2c32 395d 2c5b 3239 2c33 302c  [7,6,29],[29,30,
00000070: 375d 2c5b 382c 372c 3330 5d2c 5b33 302c  7],[8,7,30],[30,
00000080: 3331 2c38 5d2c 5b39 2c38 2c33 315d 2c5b  31,8],[9,8,31],[
00000090: 3331 2c33 322c 395d 2c5b 3130 2c39 2c33  31,32,9],[10,9,3

But you could use anything that you can parse and convert to numpy arrays really.

[Girishchandra-Yendargaye]

@ninnghazad Thank you for this information....I assume lines in /mesh_0.triangles.json and mesh_0.boundaries.json are index numbers of mesh_0.points.json starting with 0 correct?

[ninnghazad]

@Girishchandra-Yendargaye
You are correct. Remember that i chose this way because i needed large rectangular meshes - ANUGA's builtin capabilities might provide easier methods to achieve what you may be looking for.

[Girishchandra-Yendargaye]

@ninnghazad I am getting below error...What does it mean..

RuntimeError: Error at anuga\abstract_2d_finite_volumes\neighbour_mesh_ext.c:55:
Segment (336, 337) does not exist

Am I missing something....

[ninnghazad]

i would have to see some data and code to tell you more than the error does. glancing at the source tells me you are trying to construct a boundary with bad parameters.

[Girishchandra-Yendargaye]

@ninnghazad Have you tried 12 crores triangles?
Does distribute has limitations?