Julia, Threads Working, Climate Skeleton Run 6 cores

Just a bit of a personal Woo-Hoo! moment. I just now got threads to stop sending me stupid error messages about ‘world times’ being different (though I’m not sure exactly how ;-) and I’ve run the skeleton Cell Climate Model top level flow of control on 6 processes (threads sort of) on the Odroid N2.

One unexpected problem was the difficulty in slowing the processing down enough to see the cores load up…

So, OK, this code is UGLY and still has a lot of commented out bits I was playing with to see what caused which broken thing, so not laughing, OK? (Well, at least not too much…)

I had to resort to “stubs” to test things. Seems attempting to call a function inside the place you wrote it can cause differences…

So first up, the stub to run the spitting out of tasks bit:

ems@OdroidN2:~/Cmodel/ZZTest$ cat stubRCells.jl             




Note the “sleep(30)” is commented out. Briefly used when I was playing around with the @async macro to delay finish of the stub so the child tasks could complete (not yet up on locking and such…)

This basically says “bring in the RCells.jl program, then run the “RCells()” function from it.

Then, the thing that does spit out the tasks.

ems@OdroidN2:~/Cmodel/ZZTest$ cat RCells.jl 
# Climate Model via Cells
# 16 Dec 2020 - E.M.SMith

using Distributed;


# @everywhere function RCells()

function RCells()
	println(" ")
	println("This is Cells - It launches per cell model runs ")
	println(" ")
	global celllist = (1:3)	
	for year in 1:2
		for month in 3:4
			for day in 5:6
#		@sync     begin
#		@parallel	for cell in celllist
		@sync	for cell in celllist
	                         	println("Spawn Cell ",cell," in year ", year," in month ", month," on day ", day)
					@time Dcell(cell)
#                          end

“using Distributed” pulls in the Julia libraries with all the parallel code stuff so you can spit out threads on multiple processes. Rcell.jl (note lower case C and singular) is the main worker process that gets spit out for each cell. We pull that code in and say to make it available.

“addprocs(4)” says to use 4 new processes in which to run things (default is 1). I don’t know if it wins, or the “julia -p 6” wins. Testing / tuning “soon” ;-)

This first ran as an @everywhere function, but when I backed that out, it still ran, so isn’t needed. @everwhere says make this function available to all threads, but since it is just the top level “out spitter”, isn’t necessary.

The function has the for loops cut back from HUGE to dinky so the test output is manageable. I used unique integers for each loop just to make it easier to see what was working right.

I was playing with @parallel, @sync, and @async. For @async, if the top level ends before the workers are done, they get scrapped (thus the timer). @sync worked fine. I’ve not done anything with @parallel yet and need to re-read exactly what it does.

Then, you can see I spit out a message saying what it’s doing in that loop cycle, and then with the @time measure how long a cell takes as I invoke the cell itself. I still need to get all the parameter passing stuff in there, but at least now things cycle.

You invoke it with a statement as to number of processes available to run at one time:

time julia -p 7 stubRCells.jl 

Here’s the end of the run statistics on time used:

    After all of that, I'm done with Tropospheric for a while.
  0.163094 seconds (94 allocations: 61.038 MiB, 17.92% gc time)

real	0m21.130s
user	0m57.412s
sys	0m6.604s

The 0.163094 etc. is a @time macro inside one of the Julia programs for just it. The “time” in front of the julia -p 6 at the Linux level gives us the summary for the whole run at 21 seconds. Note we used 57 CPU seconds in that 21 seconds. Parallel is good ;-)

The Rcell.jl does the work, one per cell. What I think got rid of the “world time” conflict was putting the include() statements just after using Distributed and with @everywhere in front of them. I’m going to try commenting out “using Distributed” as it ought to work based on the calling program (I think, maybe). Then also the @everywhere (possibly)…

I briefly thought pf playing with this as a “module” instead of a function, so there is a commented out “module” definition line. (Playing for later ;-) Then there’s some commented out stuff from earlier trials.

ems@OdroidN2:~/Cmodel/ZZTest$ cat Rcell.jl 
# Rcell.jl
# Climate Model via Cells
# 16 Dec 2020 - E.M.SMith

using Distributed

@everywhere include("Surface.jl")
@everywhere include("SubSurface.jl")
@everywhere include("AirLayers.jl")

function Dcell(cellno)
#module Dcell
	println("You have a cell! ",cellno)

# Commenting out Surface() fixes world age for that one fn.  Why happening?
#	Surface()
# This didn't change anything #
#	Base.invokelatest(Surface())

So now that works, and it properly calls surface, subsurface, and airlayers that call all their functions as well.

I briefly tried telling things to sleep to get different completion times:

ems@OdroidN2:~/Cmodel/ZZTest$ grep sleep *
Stratos.jl:#	sleep(rand(0:2))
SubSurface.jl:#	sleep(rand(0:4))
Surface.jl:#	sleep(rand(0:4))
Tropos.jl:#	sleep(rand(0:8))

Which it did, but without any clue as to CPUs in use. So I added some workload. I also removed printed blank lines in the output as I was getting a LOT in testing ;-) You can also see the prior run “sleep(foo)” commented out as load was being added.

ems@OdroidN2:~/Cmodel/ZZTest$ cat Surface.jl 
# Climate Model via Cells
# 16 Dec 2020 - E.M.SMith

function Surface()
#	println(" ")
	println("You have reached Surface Physics. ")
#	sleep(rand(0:4))
#	println(" ")
	A * 2.4

# Debugging self call. Works.

The odd thing is the “self call” works in single threads, but seems to cause some kind of world time issue in multi-threads. Concurrency and multiple instances can be weird.

Adding math like that in 4 places only barely slowed it down. Create the array A and stuff it with random numbers between 1:100, in a 1000 x 1000 array. So 1 million elements. Then multiply each element by 2.4 (and take an int to float hit too).

I tried bumping that up to 10000 x 10000 but it turns out that 100,000,000 x 8 bytes for 64 bit is just a shade under a GB, then having a half dozen of those… well, I hit a memory / swap wall ;-}

So added the math work.

In Conclusion

So now, over THAT particular hump, I’m at the point where I need to start filling in this “toy” framework with actual physics. One Lump At A Time.

I’m going to start with adding parameter passing (a bit trickier than usual as Julia really loves to make everything local to functions and inside loops and… So a bit of work just to figure out where your variable FOO turned into a 2nd FOO that doesn’t hold what you put into FOO as it is local (again) where you wanted to work on FOO #1 data…

Then it will be Sun Time. So I’ll work on “build world” to make an actual world, and then put sun into each panel (world not turning yet, just doing trig work). Then set it in motion in an airless dry world of more or less hexagonal facets jumping one hour at a time…

Then add radiative physics in the airless world.

After that, it will be just one step at a time adding substance to each of the other chunks called (AND working out all the interprocess communications…)

FWIW, run time for this was still somewhat dominated by Julia set up / compile time. Had to add those array math bits to get it slow enough. You get a big load of all cores for about 2 seconds as it has a bit of a think, then it starts cranking on work product. It took 4 of them Array / Math blocks to slow it down enough to get decent cpu use observable, and even then it looks like it finishes the math and hangs around a while as the print to screen happens.

ems@OdroidN2:~/Cmodel/ZZTest$ grep "A=rand" *
Stratos.jl:	A=rand(1:100,1000,1000)
SubSurface.jl:	A=rand(1:100,1000,1000)
Surface.jl:        A=rand(1:100,1000,1000)
Tropos.jl:        A=rand(1:100,1000,1000)

So, with that, I’m off to bed. (Always head to bed after a big success at coding. If you decide to “try one more thing”, before you know it it will be 6 am and the spouse is asking “Have you been up all night again?”. Why? Because you are SURE it is just one more little “one line fix” and you will have it!!! And it isn’t.

Subscribe to feed

About E.M.Smith

A technical managerial sort interested in things from Stonehenge to computer science. My present "hot buttons' are the mythology of Climate Change and ancient metrology; but things change...
This entry was posted in AGW Science and Background, Tech Bits. Bookmark the permalink.

38 Responses to Julia, Threads Working, Climate Skeleton Run 6 cores

  1. V.P. Elect Smith says:

    I changed “addprocs” to 6 and re-ran it:

      0.162071 seconds (92 allocations: 61.037 MiB, 17.52% gc time)
    real	0m21.646s
    user	1m1.224s
    sys	0m7.672s

    System time up a second, user time up about 4 seconds. Real just a tad more. Setting it down to 2:

        After all of that, I'm done with Tropospheric for a while.
      0.162100 seconds (94 allocations: 61.038 MiB, 17.46% gc time)
    real	0m19.514s
    user	0m46.548s
    sys	0m5.932s

    Sys time drops about .7 seconds (from the 4 case) and user time reduces as well. Very small change in elapsed (‘real’) time, but still a bit (~1.5 sec.) faster.

    So the “addprocs” is controlling inside the setting of “julia -p 6” telling it to use all 6 cores on this board.

    Setting both to “9” (so over the number of total cores by 50%) the htop listing seems to show peak memory used higher and all cores fully saturated during the compute phase. But it goes by fast and is hard to say for sure:

      0.158502 seconds (94 allocations: 61.039 MiB, 17.28% gc time)
    real	0m27.523s
    user	0m14.016s
    sys	0m1.496s

    Real time gets an added 6 to 7 seconds, user time (user CPU time) drops and System time goes way down. Odd that. I’d expect more system time dealing with more jobs. Then again, it’s spitting them out in bigger batches? Less time spent waiting on @sync to say “OK, spit out another one”?

    So a couple of conclusions:

    1) This toy system is still way too lite on absolute computes per cell to really amortize the full set-up and deploy cost of parallel execution effectively. That will change as the cells grow in complexity.

    2) Parallel execution can memory limit much faster as N cells at once uses N x memory at once. Watching memory use in parallel operation is important. Boards with extra memory as important as boards with added cores. Memory per Core is the metric to look at when buying compute boards. At least 1 GB.

    2.b) Large arrays eat a LOT of memory FAST. Avoid a lot of big arrays as shared data structures. Arrays / memory can be traded for disk I/Os (and vice versa).

    3) Tuning matters a LOT. Just that one parameter, processes, gave a bunch of different results. When it gets to having separate threads inside some of the sub-processes and / or parceling them out to local vs remote compute boards, then the communications load will need to be considered as well. Julia has facilities for choosing threads vs processes (shared memory / resources vs not) and various types of communication between processes. All that will need / benefit from tuning selections.

    Time for more morning coffee ;-)

  2. V.P. Elect Smith says:

    Just for grins, I set processes to 12 for both the command line and addprocs, this resulted in minor memory saturation and swapping of about 400 MB.

      0.142997 seconds (92 allocations: 61.037 MiB, 5.85% gc time)
    real	0m51.882s
    user	1m33.636s
    sys	0m14.784s

    The “last cell run” didn’t take much different time at 14.29, but overall, the “real” is almost double, pushing a minute. User time blew out to 1.5 minutes (so about 50% more) and Systems time roughly doubled from most of the prior runs (but still is a small part of the total).

    Clearly running too many processes at once for a given bit of hardware is a Very Bad Thing and to be avoided. Running too few (at least with these light weight processes) is not as much of an issue.


    Oh, and it is great fun to watch HTOP show the machine 100% to the wall on computes, in any of the process number settings ;-) Julia (the compiler) seems well aware of multiple cores and “loads them up” during the “Just In Time” compile step at launch.

    OK, I’m going to set it back to 5 or 6 (for this 6 core machines) and leave it there while I work on the next step / iteration of the code. In theory, 5 is the right number. That leaves one core for the Operating System, your terminal session, and other stuff. On a dedicated compute engine board with most likely 4 cores, likely best is 4 (or equal to number of cores) as it isn’t doing much else.

  3. V.P. Elect Smith says:

    I’ve added a passed parameter to turn on / off the printing of stuff. Run the code each way. Not much changes in the statistics (even with the ‘if then’ being done):

    With printing on:

       You are in the Troposphere Convective Zone. 
            You have reached Troposheric Self Absorption.
            So just deal with it.  All that IR and stuff
            You are a rising ball of hot air. 
            Send in the Clouds, there must be Clouds. 
            You are all wet, or getting there. 
            You have reached Air RH & T Adjustment. 
            After all we've been through, we need to adjust.
            You have ridden the Winds... 
        After all of that, I'm done with Tropospheric for a while.
      0.163244 seconds (94 allocations: 61.038 MiB, 17.61% gc time)
    real	0m20.806s
    user	0m50.772s
    sys	0m6.980s


    Spawn Cell 2 in year 2 in month 4 on day 5
      0.143977 seconds (16 allocations: 61.036 MiB, 4.91% gc time)
    Spawn Cell 3 in year 2 in month 4 on day 5
      0.155338 seconds (16 allocations: 61.036 MiB, 16.74% gc time)
    Spawn Cell 1 in year 2 in month 4 on day 6
      0.369785 seconds (16 allocations: 61.036 MiB, 65.11% gc time)
    Spawn Cell 2 in year 2 in month 4 on day 6
      0.143047 seconds (16 allocations: 61.036 MiB, 4.96% gc time)
    Spawn Cell 3 in year 2 in month 4 on day 6
      0.155367 seconds (16 allocations: 61.036 MiB, 16.73% gc time)
    real	0m20.817s
    user	0m50.824s
    sys	0m7.060s

    I left on the printing that a cell had been spawned. Need some kind of info that it really is doing something ;-)

    The stub for testing now passes in a parameter p for printing:

    ems@OdroidN2:~/Cmodel/ZZTest$ cat stubRCells.jl 

    This, then gets passed on down the stack of things. I can, if desired, force on any one function printing via resetting p internal to it. This means I can focus “diagnostic prints” on just one function under work. Or run it all wide open and chatty.

    Here’s what “Clouds” looks like:

    ems@OdroidN2:~/Cmodel/ZZTest$ cat Clouds.jl 
    # Climate Model via Cells
    # Cloud Function
    # 16 Dec 2020 - E.M.SMith
    function Clouds(p::Char)
    	if p == 'y'
    #	println(" ")
    	println("        Send in the Clouds, there must be Clouds. ")
    #	println(" ")

    No, I didn’t change the indent level of the actual print statements. A lot of pointless effort doing nothing that matters so didn’t bother.

    OK, parameter passing (both cell number as a non-type stated number and print or not as a single typed Char) being done. Guess it’s time to take on “BuildWorld” and start some real physics…

  4. V.P. Elect Smith says:

    I’m not sure, but I may have the hard bit of BuildWorld done.

    Here’s the stub for testing:

    ems@OdroidN2:~/Cmodel/ZZTest$ cat stubBuild.jl 

    Called with p as y to print diagnostics, and with a world size of 100 points, it gave me 100 points (both in polar and x,y,z form – still need to convert to Lat and Long). And I need to do a QA pass. And all that “who’s your neighbor” stuff, but the basic “dot placement” seems right, and looks like a decent interpretation of:

    Click to access sphere_equi.pdf

    Here’s the top and bottom of a run:

    ems@OdroidN2:~/Cmodel/ZZTest$ julia stubBuild.jl 
    You have reached BuildWorld. 
    For each Cell, place it on the Globe 
    Calculate the neighbors, angles, inclination,etc 
    Stuff the database with calculated values 
    just before for loop
    type IncHalf Float64
    type Theta Float64
    type Temp Float64
    type Mphi Float64
    Just before inner loop:
    X= 0.17364817766693033   Y= 0.0   Z= 0.984807753012208
    X= -0.08682408883346512   Y= 0.1503837331804353   Z= 0.984807753012208
    X= -0.08682408883346525   Y= -0.15038373318043524   Z= 0.984807753012208
    type IncHalf Float64
    type Theta Float64
    type Temp Float64
    type Mphi Float64
    Just before inner loop:
    X= 0.49999999999999994   Y= 0.0   Z= 0.8660254037844387
    X= 0.38302222155948895   Y= 0.32139380484326957   Z= 0.8660254037844387
    X= 0.0868240888334652   Y= 0.49240387650610395   Z= 0.8660254037844387
    X= -0.24999999999999986   Y= 0.4330127018922193   Z= 0.8660254037844387
    X= -0.4698463103929541   Y= 0.1710100716628344   Z= 0.8660254037844387
    X= -0.46984631039295416   Y= -0.1710100716628343   Z= 0.8660254037844387
    X= -0.25000000000000017   Y= -0.4330127018922192   Z= 0.8660254037844387
    X= 0.08682408883346497   Y= -0.492403876506104   Z= 0.8660254037844387
    X= 0.38302222155948884   Y= -0.32139380484326974   Z= 0.8660254037844387
    type IncHalf Float64
    type Theta Float64
    type Mphi Float64
    Just before inner loop:
    X= 0.5000000000000003   Y= 0.0   Z= -0.8660254037844385
    X= 0.3830222215594893   Y= 0.32139380484326985   Z= -0.8660254037844385
    X= 0.08682408883346526   Y= 0.49240387650610434   Z= -0.8660254037844385
    X= -0.25000000000000006   Y= 0.43301270189221963   Z= -0.8660254037844385
    X= -0.4698463103929545   Y= 0.17101007166283455   Z= -0.8660254037844385
    X= -0.46984631039295455   Y= -0.17101007166283444   Z= -0.8660254037844385
    X= -0.2500000000000004   Y= -0.4330127018922195   Z= -0.8660254037844385
    X= 0.08682408883346504   Y= -0.4924038765061044   Z= -0.8660254037844385
    X= 0.3830222215594892   Y= -0.32139380484327   Z= -0.8660254037844385
    type IncHalf Float64
    type Theta Float64
    type Temp Float64
    type Mphi Float64
    Just before inner loop:
    X= 0.17364817766693028   Y= 0.0   Z= -0.984807753012208
    X= -0.0868240888334651   Y= 0.15038373318043524   Z= -0.984807753012208
    X= -0.08682408883346522   Y= -0.15038373318043521   Z= -0.984807753012208

    I still have a type checking block printing (had a ‘type’ issue earlier that was in fact a syntax error).

    The commented out stuff in the middle is where I had NumCells hard set at 100 prior to putting the parameter in place, and set “pie” to a fixed value of Pi. Julia was giving me an irrational error message ;-) Something like “Irrational can not be called as a function”. Julia also lets you leave out the “*” for multiply operator, and I’d decided to try that. Seems you can NOT do that before an expression in parenthesis as then it interprets the prior variable (or constant) name() as a function call… So I had:

    pi(stuff operator stuff) and it took that as a function instead of a multiply. OK, “don’t do that” if it hurts… So pie(stuff operator stuff) changed to “float64 can not be called…” and I got a whack with the clue stick… Thus some “left overs” from my sorting that out (along with pulling (stuff op stuff) out to a separate operation). So some ‘re-condensing’ to be done now too.

    just before for loop
    type IncHalf Float64
    ERROR: LoadError: MethodError: objects of type Irrational{:��} are not callable
     [1] BuildWorld(::Char) at /home/ems/Cmodel/ZZTest/BuildWorld.jl:69
     [2] top-level scope at none:0

    I decorated it with ‘typeof’ statements looking for what was irrational, forgetting about pi… a built in of type irrational. Oh, and I originally had “Mtheta” as what is now marginaltheta. So the comments a bit curious talking about Mtheta. Did I mention this is not cleaned up at all?

    What I’ve go so far:

    ems@OdroidN2:~/Cmodel/ZZTest$ cat BuildWorld.jl 
    # Climate Model via Cells
    # BuildWorld 
    # 16 Dec 2020 - E.M.SMith
    function BuildWorld(p::Char,NumCells::Int)
    	if p == 'y'
    	println(" ")
    	println("You have reached BuildWorld. ")
    	println("For each Cell, place it on the Globe ")
    	println("Calculate the neighbors, angles, inclination,etc ")
    	println("Stuff the database with calculated values ")
    	println(" ")
    #global pie=3.1415926
    # Begin actual devo code:
    # Counter is number of dots placed.
    # Take the area of your sphere (via 4Pi*R^2) and divide it by 
    # the number of cells to make to get area / cell.
    # Distance between cells is the square root of the Area per cell.
    # (Really? Why again?)
    	Counter= 0
    	AreaCell= 4pi*R^2/NumCells
    	DistCell= sqrt(AreaCell)
    # Physics uses Theta as Z axis or Up/down and 
    # Phi as x / y axis or round the middle
    # Math Folks use Theta are round the middle and Theta as up / dn
    # This probaly uses the Physics / iso std whth Theta Polar angle 
    # or angle from axis of rotation and Phi the longitude rotation
    # about the axis.
    # What the hell is Mtheta?  The increment in Theta per cell?
    # Marginal Theta.  Pi is from pole to pole.  Pi / Disance is 
    # unit of increment poleward.
    # Then dTheta is pi / Mtheta?  Why again?  and if dPhi is 
    # roughly equal to dTheta, why find both?
    # Is Mtheta the number to plant in any one latitude band?
    # Or the number or latitude bands?  Or both?
    	marginaltheta = round(pi/DistCell)
    	dTheta = pi/marginaltheta
    	dPhi = AreaCell/dTheta
    #	d in Theta direction ~= d in Phi direction
    # so we fina a new Theta for a given Mtheta (count of theta increments?)
    # by taking how many we did, plus a half, times pi.  Then dividing by
    # the increments.
    # then Mphi looks like a simlar increment in phi.
    println("just before for loop")
    	global Theta
    	global phi
    #	Theta=.4
    	global marginaltheta
    	for m in 0:limit
    		println("type IncHalf ",typeof(IncHalf))
    		println("type Theta ",typeof(Theta))
    		println("type Temp ",typeof(Temp))
    		global Mphi= round(Temp)
    		println("type Mphi ",typeof(Mphi))
    		println("Just before inner loop:")
    			for n in 0:limit2 
    				phi = 2pi*n/Mphi
    				x= sin(Theta)*cos(phi)
    				y= sin(Theta)*sin(phi)
    				z= cos(Theta)
    				println("X= ",x,"   Y= ",y,"   Z= ",z)
    				Counter += 1
    #	Cracked Egg / Backwards G is Theta
    #	Lollypop / Circular T is Phi

    Oh, and a reminder to myself what the “cursive” Greek letters in the cited paper translate to in more familiar Greek. ;-)

    Basically I’m getting 10 per line of Lat instead of variable / line of Lat.

    I wonder if it is time for the run to the store for a celebration bottle of vino?

    Hmmm… Maybe… (votes welcome ;-)

  5. V.P. Elect Smith says:

    Gee… Building a “World” is remarkably fast. Just to see how much all that Trig was costing, I decided to build a world of 10,000 cells. That wasn’t enough to notice, so went to 100,000:

    ems@OdroidN2:~/Cmodel/ZZTest$ ./dobuild 
    World of 10000 count cells
    real	0m8.172s
    user	0m21.048s
    sys	0m4.020s
    ems@OdroidN2:~/Cmodel/ZZTest$ vi stubBuild.jl 
    ems@OdroidN2:~/Cmodel/ZZTest$ ./dobuild 
    World of 100000 count cells
    real	0m8.344s
    user	0m21.360s
    sys	0m4.100s

    Looks like building the Cell World is a cheap operation even for many cells.

    I’m still going to limit early runs to smaller numbers. Both so that I can get many runs in even as code complexity builds, but also so some cases are of a scale where “doing it by hand to check” still can work ;-)

    Well, I’m going to just putter with the code, shrinking and polishing it for a while, then trying to figure out a good way to add “what cells are adjacent” (and track them) along with that whole LAT LONG thing…

  6. Compu Gator says:

    Wey-ell, I can at least congratulate you.

    But your reports start me thinking analytically at times when what I really need to do is let my brain slow down for drifting off to sleep. So my congratulations include mixing & hoisting an adult beverage in your honor.

    Like how to subdivide the globe for calculating those nonEuclidean angles of direction from soap-bubble to soap-bubble. Or not.[🌐]

    So indeed, beware the inclination “hey!  I’ve still got enough time to try one more interesting little thing”. So you start it before you realize that today, you haven’t actually fixed yourself dinner. Or lunch nor breakfast.

    Note 🌐 :  https://chiefio.wordpress.com/2020/11/08/there-is-no-good-coordinate-system-for-climate-models/.

  7. V.P. Elect Smith says:

    The author of the paper noted that ~”it looks like another dot could be put at the pole”. I think this is a fault in his algorithm not handling the single point case properly. I’ve put in “edge case handling” where he just finagled the math some to avoid it. I now get ONE dot at each pole, six around it (if you use a modulo 6 number of cells, at least so far that has worked) and then various numbers in between.

    Here’s the (somewhat) cleaned up and extended BuildWorld code and result.

    First, the top and end of the run for 215 or 216 target cells (Hey, I was incrementing, but it didn’t seem to change much). I also added a printout of the values of Theta, Phi, and the cell count. Notice the first cell is at x=y=0 or North Pole, and Z=1 or full radius straight up. Then you get 6 cells in a nice hexagonal grid around it, then 12. (Then it goes to 17, 21, 24, 26, 26, 24, 21, 17, 12, 6, 1 and you are done. Only the temperate to tropical bands have a “hexagonal discontinuity” to deal with, but they also have linear bands for linear wind flow.


    ems@OdroidN2:~/Cmodel/ZZTest$ ./dobuild 
    World of 215 count cells
    Theta 0.0 Mtheta 13.0 Mphi 1
    Just before inner loop:
    Theta 0.0 Phi 0.0 Count 1
    X= 0.0   Y= 0.0   Z= 1.0
    Theta 0.241660973353061 Mtheta 13.0 Mphi 6.0
    Just before inner loop:
    Theta 0.241660973353061 Phi 0.0 Count 2
    X= 0.23931566428755774   Y= 0.0   Z= 0.970941817426052
    Theta 0.241660973353061 Phi 1.0471975511965976 Count 3
    X= 0.1196578321437789   Y= 0.20725344479657334   Z= 0.970941817426052
    Theta 0.241660973353061 Phi 2.0943951023931953 Count 4
    X= -0.11965783214377881   Y= 0.20725344479657337   Z= 0.970941817426052
    Theta 0.241660973353061 Phi 3.141592653589793 Count 5
    X= -0.23931566428755774   Y= 2.9307716225558014e-17   Z= 0.970941817426052
    Theta 0.241660973353061 Phi 4.1887902047863905 Count 6
    X= -0.11965783214377898   Y= -0.20725344479657332   Z= 0.970941817426052
    Theta 0.241660973353061 Phi 5.235987755982989 Count 7
    X= 0.1196578321437789   Y= -0.20725344479657334   Z= 0.970941817426052
    Theta 0.483321946706122 Mtheta 13.0 Mphi 12.0
    Just before inner loop:
    Theta 0.483321946706122 Phi 0.0 Count 8
    X= 0.4647231720437685   Y= 0.0   Z= 0.8854560256532099
    Theta 0.483321946706122 Phi 0.5235987755982988 Count 9
    X= 0.40246207271719014   Y= -0.23236158602188473   Z= -0.8854560256532096
    Theta 2.8999316802367323 Mtheta 13.0 Mphi 6.0
    Just before inner loop:
    Theta 2.8999316802367323 Phi 0.0 Count 208
    X= 0.23931566428755768   Y= 0.0   Z= -0.970941817426052
    Theta 2.8999316802367323 Phi 1.0471975511965976 Count 209
    X= 0.11965783214377887   Y= 0.2072534447965733   Z= -0.970941817426052
    Theta 2.8999316802367323 Phi 2.0943951023931953 Count 210
    X= -0.11965783214377879   Y= 0.20725344479657332   Z= -0.970941817426052
    Theta 2.8999316802367323 Phi 3.141592653589793 Count 211
    X= -0.23931566428755768   Y= 2.930771622555801e-17   Z= -0.970941817426052
    Theta 2.8999316802367323 Phi 4.1887902047863905 Count 212
    X= -0.11965783214377895   Y= -0.20725344479657326   Z= -0.970941817426052
    Theta 2.8999316802367323 Phi 5.235987755982989 Count 213
    X= 0.11965783214377887   Y= -0.2072534447965733   Z= -0.970941817426052
    Theta 3.1415926535897936 Mtheta 13.0 Mphi -0.0
    Just before inner loop:
    Theta 3.1415926535897936 Phi 0 Count 214
    X= -3.216245299353273e-16   Y= -0.0   Z= -1.0
    real	0m8.585s
    user	0m21.304s
    sys	0m4.088s

    Note that last X at e^-16. That’s essentially zero. And Z=-1.0 is straight south.

    Here’s the code as it stands now:

    ems@OdroidN2:~/Cmodel/ZZTest$ cat BuildWorld.jl 
    # Climate Model via Cells
    # BuildWorld 
    # 22 Dec 2020 - E.M.SMith
    function BuildWorld(p::Char,NumCells::Int)
    	if p == 'y'
    	println(" ")
    	println("You have reached BuildWorld. ")
    	println("For each Cell, place it on the Globe ")
    	println("Calculate the neighbors, angles, inclination,etc ")
    	println("Stuff the database with calculated values ")
    	println(" ")
    # Begin actual devo code:
    # Counter is number of dots placed.
    # Take the area of your sphere (via 4Pi*R^2) and divide it by 
    # the number of cells to make to get area / cell.
    # Distance between cells is the square root of the Area per cell.
    # (Really? Why again?)
    	Counter= 0
    	AreaCell= 4pi*R^2/NumCells
    	DistCell= sqrt(AreaCell)
    # Physics uses Theta as Z axis or Up/down and 
    # Phi as x / y axis or round the middle
    # Math Folks use Theta are round the middle and Theta as up / dn
    # This probaly uses the Physics / iso std whth Theta Polar angle 
    # or angle from axis of rotation and Phi the longitude rotation
    # about the axis.
    # Mtheta is the increment in Theta per row of cells
    # Marginal Theta.  Pi is from pole to pole.  Pi / Distance is 
    # unit of increment poleward.
    # Then dTheta is pi / Mtheta. Pi is one arc, pi/Mtheta is 
    # the delta of theta in the polward direction in one band of cells
    # and if dPhi is roughly equal to dTheta, why find both?
    	Mtheta = round(pi/DistCell)
    	dTheta = pi/Mtheta
    	dPhi = AreaCell/dTheta
    #	d in Theta direction ~= d in Phi direction
    # so we fina a new Theta for a given Mtheta (count of theta increments)
    # by taking how many we did, plus a half, times pi.  Then dividing by
    # the increments.
    # then Mphi looks like a simlar increment in phi.
    if p == 'y'
    println("just before for loop")
    	global Theta
    	global phi
    	global Mtheta
    	for m in 0:(Mtheta)
    #		Theta=(pi*(m+0.5))/Mtheta
    		if m == 0
    			global Mphi= round(2pi*sin(Theta)/dPhi)
    		println(" ")
    		println("Theta ",Theta," Mtheta ",Mtheta," Mphi ",Mphi)
    #		if p == 'y'
    		println("Just before inner loop:")
    #		end
    		if Mphi == 0
    			x= sin(Theta)*cos(phi)
    			y= sin(Theta)*sin(phi)
    			z= cos(Theta)
    			Counter += 1
    			println("Theta ",Theta," Phi ",phi," Count ",Counter)
    			println("X= ",x,"   Y= ",y,"   Z= ",z)
    		for n in 0:Mphi-1
    			phi = 2pi*n/Mphi
    			x= sin(Theta)*cos(phi)
    			y= sin(Theta)*sin(phi)
    			z= cos(Theta)
    			Counter += 1
    			println("Theta ",Theta," Phi ",phi," Count ",Counter)
    			println("X= ",x,"   Y= ",y,"   Z= ",z)
    #	Cracked Egg / Backwards G is Theta
    #	Lollypop / Circular T is Phi

    I need to think on this a bit to see if I can find a way to “clean it up” in terms of incorporating the exception cases into the regular loops, but “we’ll see”.

    Also, all this polar coordinate math is nice and all, but now I’ve got to (maybe) take a conversion to LAT and LON for compatibility with all the site data and general understanding of the world by mere mortals… So I’m thinking maybe a “do over” with the same idea but using LAT and LON as the fundamental units… and R as actual Earth R and lines of Latitude as actual Earth size… Program gets a bit messier with that approach, but the result is LAT and LONG and areas you can use on a map…

  8. V.P. Elect Smith says:


    As you can see, I just did my (Best Colombo voice:) “One More Thing”…

    I’ve got about 3/4 of a bottle of COSTCO Pinot Grigio down the gullet to wash down the microwaved Pot Pie I had for dinner (Yeah, coding frenzy, toss some crap from the freezer in the microwave while you noodle out the next change and eat it while typing new code with one hand and… Sigh. I am a Geek…)

    So unlikely to try one more one more thing tonight… other than one more glass…

  9. V.P. Elect Smith says:

    I just realized I’ve not updated the comments in the code from when I was puzzling out what the article Author’s Greek might mean. In the printout of the results, you can see that Mtheta is the number of bands of latitude over which dots are to be distributed, and Mphi is the number of dots worth of space in each band of latitude. The rest is just bin filling really. Why did he use M for number? Who knows.

    But by handling the ‘edge cases’ where Mphi is zero, at the poles, we get a very nice single dot at each pole, then a nice 6 hexagons close packing around it, and then 12 as hexagons double at each row out. Until the curvature of a sphere requires the tiling to curve from a nearly flat plain into a downward direction. Then the count breaks from a doubling to less than double. As we approach the equator, the number of dots / latitude band changes much less, as we are approaching a cylinder aspect (which would have equal tiles per row).

    This is essentially what I was doing on paper when pondering No Good Grid in that prior posting. I’d not recognized it as the article Author didn’t start with a dot at the pole, so didn’t get a hex tiling near the pole.

    Which brings up a minor pondering…

    With “one at the pole”, you can have air descending and being sent in all directions, all points of the compass. That’s a little hard to handle in code looking mostly at lateral motion. “Negative Convection” will matter a great deal and generally it is not handled as a primary thing in winds code.

    With “several near the pole”, you have more opportunities for wind to “cross the polar area” in all directions between those cells. It doesn’t need to be handled inside one cell, but can be emergent based on many cells.

    Which is better? Which is easier to code?

    At this point, I don’t know. I’m just noting that there are two different ways to layout the dots on the globe and the results might change as well as the difficulty of “getting it right” in the code to model air flow.

    Basically, does introducing an “edge case” of one dot AT the pole, cause edge cases to crop up in all the future code to handle air flow at/near the pole?

  10. p.g.sharrow says:

    @EMSmith, thinking about this thing, it appears to me that your fluid flow at “0” latitude would be biased down from vertical toward Equator then vector modified by drag from Earth’s rotation. Inertia/gravity from space always pulls against the fluid spinup so there is always is a constant there. Those flow vectors are then modified by terrain followed by energy in/out that cause local changes, more vertical then horizontal. The sun – input, and space – output, of energy are fairly constant but locally modified.
    So there are a lot of fixed physics that are locally modified, all vector modified by the Earth’s rotation and convection.
    So starting with a smooth bare rock, rotating, under vacuum, near a star/furnace, in dead cold space is the place to start.

  11. p.g.sharrow says:

    Maybe even a rotational locked planet in circle orbit around it’s star would be the starting point as rotation, distance and inclination all must at some point be included.
    This sounds like calculating electronic flows in circuits, Inductance, resistance, capacitance, magnetics all have their place and formula,
    So each must have it’s place, even if at the start it is a “constant” of “1” to be formulated later.

  12. V.P. Elect Smith says:

    Hey, I got Lat and Lon into it now:

    Lat=  -62.3  Lon=  210.0  Count=      203 
    Lat=  -62.3  Lon=  240.0  Count=      204 
    Lat=  -62.3  Lon=  270.0  Count=      205 
    Lat=  -62.3  Lon=  300.0  Count=      206 
    Lat=  -62.3  Lon=  330.0  Count=      207 
    Theta 2.900 Mtheta     13 Mphi      6
    Lat=  -76.2  Lon=  0.0  Count=      208 
    Lat=  -76.2  Lon=  60.0  Count=      209 
    Lat=  -76.2  Lon=  120.0  Count=      210 
    Lat=  -76.2  Lon=  180.0  Count=      211 
    Lat=  -76.2  Lon=  240.0  Count=      212 
    Lat=  -76.2  Lon=  300.0  Count=      213 
    Theta 3.142 Mtheta     13 Mphi      0
    Lat= -90.00000000000003 Lon= 0
    Theta 3.142 Phi 0.000        Count     214 
    X= -0.000    Y=  -0.000    Z=  -1.000 
    real	0m10.064s
    user	0m23.088s
    sys	0m3.812s


    Yeah, a lot of stuff “TBD” later…

    Initially, I’m just working on making the facets (cells) on the globe. (Still to do is figure their “tilt” vs the polar axis, then place a Sun in one direction and figure their angle to the solar illumination).

    At first there will be no orbit (no orbital mechanics). So things like nutation and change of orbital eccentricity will be zero. Also, first cut will have zero rotation of the Earth during development, so essentially your rotational locked planet (with a near infinite orbital period ;-)

    You can think of it as a ball made of flat nearly hexagonal tiles sitting in space.

    I’ll then try to make it rotate properly so tiles angle to sun changes in steps (likely 24 of them to start, so “jumps” 24 times to 24 different rotational positions, like an infinitely fast step motor) “Someday” I might try to work out some continuous function math instead of steps… maybe… someday… 8-{

    Once I have a rotating ball of facets and a Sun, then I’ll work on adding rock heat sink, subsurface heat flow outward, and radiation to space. Basically a Moon Analog. Once that’s working, comes the “fun” part. Adding all the working fluids (and occasional solids – snow & ice) on our planet… I could be iterating on that bit of programming for a decade…

    At some point the runs will slow down so much that development on one SBC isn’t feasible, and then I’ll actually start distributing cells to other boards. When that happens, I’ll also play with a planet of much larger cell count (more facets at way more angles…). Then I’ll be in a race condition of precision vs my income to buy more boards ;-) (But I’ve got about a dozen to get busy before that’s an issue …)

    One bit I ought to look at, but haven’t, is the degree to which Julia cares about bit length in your computer. I think I can distribute work to a heterogeneous environment (32 bit, 64 bit) and maybe even a mix of PC and ARM. But that’s TBD much much later… about the time my 8? 64 Bit SBCs are all loaded up and cranking.

    When I think about the scope of all of it, it can be a bit intimidating. Like all big software projects, you just pick a little corner, work on it, and start building out… only on rare occasions looking wayyyyy over there at the far corner…

  13. V.P. Elect Smith says:

    Running BuildWorld with N=6000 still gives “one at the pole and 6 around it” so that can hold well.

    Theta 3.096 Mtheta     69 Mphi      6
    Lat=  -87.4  Lon=  0.0  Count=     5996 
    Lat=  -87.4  Lon=  60.0  Count=     5997 
    Lat=  -87.4  Lon=  120.0  Count=     5998 
    Lat=  -87.4  Lon=  180.0  Count=     5999 
    Lat=  -87.4  Lon=  240.0  Count=     6000 
    Lat=  -87.4  Lon=  300.0  Count=     6001 
    Theta 3.142 Mtheta     69 Mphi      0
    Lat= -90.0 Lon= 0
    Theta 3.142 Phi 0.000        Count    6002 
    X= 0.000    Y=  0.000    Z=  -1.000 
    real	0m11.249s
    user	0m23.884s
    sys	0m4.220s

    It looks like about 2.6 degrees on a side.

    Lat=  -30.0  Lon=  344.7  Count=     4556 
    Lat=  -30.0  Lon=  347.8  Count=     4557 
    Lat=  -30.0  Lon=  350.8  Count=     4558 
    Lat=  -30.0  Lon=  353.9  Count=     4559 
    Lat=  -30.0  Lon=  356.9  Count=     4560 
    Theta 2.140 Mtheta     69 Mphi    115
    Lat=  -32.6  Lon=  0.0  Count=     4561 
    Lat=  -32.6  Lon=  3.1  Count=     4562 
    Lat=  -32.6  Lon=  6.3  Count=     4563 

    Here’s the distribution by latitude band:

    ems@OdroidN2:~/Cmodel/ZZTest$ ./dobuild 
    World of 6000 count cells
    Theta 0.000 Mtheta     69 Mphi      1
    Theta 0.046 Mtheta     69 Mphi      6
    Theta 0.091 Mtheta     69 Mphi     12
    Theta 0.137 Mtheta     69 Mphi     19
    Theta 0.182 Mtheta     69 Mphi     25
    Theta 0.228 Mtheta     69 Mphi     31
    Theta 0.273 Mtheta     69 Mphi     37
    Theta 0.319 Mtheta     69 Mphi     43
    Theta 0.364 Mtheta     69 Mphi     49
    Theta 0.410 Mtheta     69 Mphi     54
    Theta 0.455 Mtheta     69 Mphi     60
    Theta 0.501 Mtheta     69 Mphi     66
    Theta 0.546 Mtheta     69 Mphi     71
    Theta 0.592 Mtheta     69 Mphi     76
    Theta 0.637 Mtheta     69 Mphi     81
    Theta 0.683 Mtheta     69 Mphi     86
    Theta 0.728 Mtheta     69 Mphi     91
    Theta 0.774 Mtheta     69 Mphi     95
    Theta 0.820 Mtheta     69 Mphi    100
    Theta 0.865 Mtheta     69 Mphi    104
    Theta 0.911 Mtheta     69 Mphi    108
    Theta 0.956 Mtheta     69 Mphi    112
    Theta 1.002 Mtheta     69 Mphi    115
    Theta 1.047 Mtheta     69 Mphi    118
    Theta 1.093 Mtheta     69 Mphi    121
    Theta 1.138 Mtheta     69 Mphi    124
    Theta 1.184 Mtheta     69 Mphi    126
    Theta 1.229 Mtheta     69 Mphi    129
    Theta 1.275 Mtheta     69 Mphi    131
    Theta 1.320 Mtheta     69 Mphi    132
    Theta 1.366 Mtheta     69 Mphi    134
    Theta 1.411 Mtheta     69 Mphi    135
    Theta 1.457 Mtheta     69 Mphi    136
    Theta 1.503 Mtheta     69 Mphi    136
    Theta 1.548 Mtheta     69 Mphi    137
    Theta 1.594 Mtheta     69 Mphi    137
    Theta 1.639 Mtheta     69 Mphi    136
    Theta 1.685 Mtheta     69 Mphi    136
    Theta 1.730 Mtheta     69 Mphi    135
    Theta 1.776 Mtheta     69 Mphi    134
    Theta 1.821 Mtheta     69 Mphi    132
    Theta 1.867 Mtheta     69 Mphi    131
    Theta 1.912 Mtheta     69 Mphi    129
    Theta 1.958 Mtheta     69 Mphi    126
    Theta 2.003 Mtheta     69 Mphi    124
    Theta 2.049 Mtheta     69 Mphi    121
    Theta 2.094 Mtheta     69 Mphi    118
    Theta 2.140 Mtheta     69 Mphi    115
    Theta 2.185 Mtheta     69 Mphi    112
    Theta 2.231 Mtheta     69 Mphi    108
    Theta 2.277 Mtheta     69 Mphi    104
    Theta 2.322 Mtheta     69 Mphi    100
    Theta 2.368 Mtheta     69 Mphi     95
    Theta 2.413 Mtheta     69 Mphi     91
    Theta 2.459 Mtheta     69 Mphi     86
    Theta 2.504 Mtheta     69 Mphi     81
    Theta 2.550 Mtheta     69 Mphi     76
    Theta 2.595 Mtheta     69 Mphi     71
    Theta 2.641 Mtheta     69 Mphi     66
    Theta 2.686 Mtheta     69 Mphi     60
    Theta 2.732 Mtheta     69 Mphi     54
    Theta 2.777 Mtheta     69 Mphi     49
    Theta 2.823 Mtheta     69 Mphi     43
    Theta 2.868 Mtheta     69 Mphi     37
    Theta 2.914 Mtheta     69 Mphi     31
    Theta 2.959 Mtheta     69 Mphi     25
    Theta 3.005 Mtheta     69 Mphi     19
    Theta 3.051 Mtheta     69 Mphi     12
    Theta 3.096 Mtheta     69 Mphi      6
    Theta 3.142 Mtheta     69 Mphi      0
    Lat= -90.0 Lon= 0
    Theta 3.142 Phi 0.000        Count    6002 
    X= 0.000    Y=  0.000    Z=  -1.000 
    real	0m10.179s
    user	0m23.740s
    sys	0m3.512s

    So is 137 around the equator enough to capture everything? 292 km long? 182 miles?

    That’s about 2 x as long as the entire San Francisco Bay Area, where San Jose can be 105 F and San Francisco 50 F on the same day ( I know, I drove it on that day and arrived without a coat… just about 60 miles apart).

    Then, 182 miles inland puts you well into the Sierra Nevada Mountains and Snow country…

    I think we’re gonna need more cells…

  14. V.P. Elect Smith says:

    To get a 29 mile cell size, long axis of the hexagon, it takes 240,000 cells…

    Theta 3.120 Mtheta    434 Mphi     19
    Theta 3.127 Mtheta    434 Mphi     13
    Theta 3.134 Mtheta    434 Mphi      6
    Theta 3.142 Mtheta    434 Mphi      0
    Lat= -90.0 Lon= 0
    Theta 3.142 Phi 0.000        Count  240009 
    X= 0.000    Y=  0.000    Z=  -1.000 
    real	0m10.201s
    user	0m23.360s
    sys	0m3.912s

    Each cutting in half of the length, takes a 4 x of number of cells. So a 15 mile cell will need about a million cells.

    Well, this program is already good for something. I can play with potential costing for each degree of granularity if one cell runs on one core… Or the inverse: How much granularity can I afford to buy, or how much time will expand if running multiple cells per core…

    As a 1/4 million was “only” 10 seconds, I think I try it with a million… Ought to be about 40 seconds or a bit less as “compile and go” only happen once per run and it’s a few seconds.

  15. V.P. Elect Smith says:

    Interesting… A lot of the elapsed time seems to be the compile and then print out steps. Likely waiting for the terminal window to scroll…

    Theta 3.120 Mtheta    886 Mphi     38
    Theta 3.124 Mtheta    886 Mphi     31
    Theta 3.127 Mtheta    886 Mphi     25
    Theta 3.131 Mtheta    886 Mphi     19
    Theta 3.135 Mtheta    886 Mphi     13
    Theta 3.138 Mtheta    886 Mphi      6
    Theta 3.142 Mtheta    886 Mphi      0
    Lat= -90.0 Lon= 0
    Theta 3.142 Phi 0.000        Count  999981 
    X= 0.000    Y=  0.000    Z=  -1.000 
    real	0m11.567s
    user	0m24.796s
    sys	0m3.956s

    So set at 1 Million, it doesn’t quite use them all. 999982 of them (one cell just after the last count increment). Mtheta is the number of latitude bands, and there are 2 x as many longitude steps at the equator ( Pi vs 2 Pi for pole to pole half circle vs around the equator full circle). So 886 + 886 = 1772 to circle the globe or about:

    40,000 / 1772 = 22.6 km long axis of hexagon
    24,900 / 1772 = 13.5 miles long axis of hexagon

    I think that’s likely fine enough to tell San Francisco from San Jose from Sacramento from Reno Nevada…

    So somewhere between 1/4 million cells and 1 million is probably about right.

    IF I ran 2 cells / core on 4 core compute boards, that would be 8 / board. So 1/32 million to 1/8 million boards. $40/board all up and running, $1.25 to $5 Million.

    I think I’ll start off with fewer cells and a lot more of them per core ;-)

  16. V.P. Elect Smith says:

    Let’s say I sink, oh, $1000 into boards. That’s 25 boards. 100 cores.

    Using the lower end 1/4 million cells, or 125,000, that’s 1,250 cells / core. Run sequentially, that’s 1250/24*60 = 1250/ 1440 or .86 minutes / cell. 52 seconds.

    So any run time over 52 seconds per time cycle-cell, I can’t run every cell even once per day. As it is currently looking like 1 model hour per cell-cycle, that’s expanding time not compressing it. I’ll need to either have run times in the single seconds (not going to happen) or a whole lot more money (also unlikely to happen) or use a lot fewer cells (certainly the case until all the code is written).

    Well, nice to know, anyway.

  17. V.P. Elect Smith says:

    Looks like a count of 10369 gives a 2 degree cell size (360 / (2 x 90)) at the equator or 180/90 for Mtheta latitude bands. That’s 120 Nautical Miles on a side as there are 60 N.Miles in a degree (one per minute of arc).

    Theta 3.072 Mtheta 90 Mphi 13
    Theta 3.107 Mtheta 90 Mphi 6
    Theta 3.142 Mtheta 90 Mphi 0
    Lat= -90.0 Lon= 0
    Theta 3.142 Phi 0.000 Count 10369
    X= 0.000 Y= 0.000 Z= -1.000

    real 0m12.583s
    user 0m22.524s
    sys 0m4.092s

    What prompted this? I was pondering what could I run for testing in a cluster. Figured I’ve got about 36 cores usable in the cluster, so take 24×60 = 1440 minutes in a day, timex 36 is 51840 core-minutes. Run each cell for 5 minutes, 51840 / 5 = 10368.

    So once I set up all my ARM boards in one big cluster (or at least 9 of them), I can run each cell for 5 minutes if I let the cluster run over a full day. That’s inside the realm of “doable” for testing / evaluating things and likely debugging / enhancing / refining the basics.

    Not going to knock anyone’s socks off with a 5 minute cell run time per physical day. OTOH, IFF I get one hour model time per time tick, and can get, oh, 1 full time tick computed per minute (per cell per core), that 5 minute run is 5 hours of “model time”. A full day in 5 days of total model run across the cluster. (5 model-hours x 5 run days = 25 model hours completed) Enough to see if things look sane or not as the day progresses.

    IFF I can make it very efficient, and get a full 1 hour model time tick computed in 12 seconds, then that one 24 hour cluster run time will make one full model day too. I hit time-sync with reality.

    So, OK, I’ve got a “target” to aim at while working on things. See if I can compute one time tick of a cell in 12 seconds elapsed time or less. Might be doable…

    I think that will be my target for testing purposes. 10368 (that uses 10369 or 10370 actual dots to finish the globe), a 1 hour model time increment (so 24 ticks in a model day), and a 12 second run time per tick/cycle. That ought to be quite good enough for testing and likely also for actually observing things of interest. (Heck, I’m already noodling out interesting scale / cost issues from this little scrap of code…)

    So that’s what I’m going to aim for first. A 2 degree cell size and 120 nautical miles across it. Ought to be some interesting stuff show up at that size Fractal Ruler ;-)

  18. V.P. Elect Smith says:

    I’ve upgraded the “print” flag to an integer, so I can turn on / off various blocks of diagnostic or informative printing with steps. (Source code later after I’m done fiddling with it today or tomorrow…)

    I really like this grid size for testing, with the exact unit degrees increments of Latitude.

    Lat=  -84.0  Lon=  284.2  Count=    10346 
    Lat=  -84.0  Lon=  303.2  Count=    10347 
    Lat=  -84.0  Lon=  322.1  Count=    10348 
    Lat=  -84.0  Lon=  341.1  Count=    10349 
    Lat=  -86.0  Lon=  0.0  Count=    10350 
    Lat=  -86.0  Lon=  27.7  Count=    10351 
    Lat=  -86.0  Lon=  55.4  Count=    10352 
    Lat=  -86.0  Lon=  83.1  Count=    10353 
    Lat=  -86.0  Lon=  110.8  Count=    10354 
    Lat=  -86.0  Lon=  138.5  Count=    10355 
    Lat=  -86.0  Lon=  166.2  Count=    10356 
    Lat=  -86.0  Lon=  193.8  Count=    10357 
    Lat=  -86.0  Lon=  221.5  Count=    10358 
    Lat=  -86.0  Lon=  249.2  Count=    10359 
    Lat=  -86.0  Lon=  276.9  Count=    10360 
    Lat=  -86.0  Lon=  304.6  Count=    10361 
    Lat=  -86.0  Lon=  332.3  Count=    10362 
    Lat=  -88.0  Lon=  0.0  Count=    10363 
    Lat=  -88.0  Lon=  60.0  Count=    10364 
    Lat=  -88.0  Lon=  120.0  Count=    10365 
    Lat=  -88.0  Lon=  180.0  Count=    10366 
    Lat=  -88.0  Lon=  240.0  Count=    10367 
    Lat=  -88.0  Lon=  300.0  Count=    10368 
    Lat= -90.0 Lon= 0
    Theta 3.142 Phi 0.000        Count   10369 
    X= 0.000    Y=  0.000    Z=  -1.000 The special case 
    real	0m11.501s
    user	0m23.880s
    sys	0m4.740s

    Longitude nicely a 60 degree hex packing just one away from the pole. 0.0 repeating often.

    So that’s settled (for testing / devo purposes…). Though I’ll likely cut it back to a much smaller number for a test run on the “tilt” calculation so I can run a couple by hand from various aspects of the tesselation…

    Well, back to the coding & polishing…

  19. Compu Gator says:

    (I drafted this text on 1 December, as a comment to the blog entry cited below [**]. I have no idea why it remained in my log, un posted. A casual scrolling shows no obvious evidence of me tearing apart what I had already written, e.g., to prepare me to rewrite it. Then today, trying to catch up, I experienced more than 1/2 dozen of the exasperating “Sorry, this comment could not be posted” messages from WordPress in 1/2 hour. So this time, I’m trying to post it as a comment here.)

    E.M. Smith posted on 8 November 2020 at (no time shown) GMT [**]:
    Per the Earth Wiki, the surface area of the Earth is: “510.072.000 km2 (196,940,000 sq mi)”.  So using 10,000 “grid cells” each one would be 19,694 square miles in area or 51,007 sq.km. That’s 140 miles on a side or 225 km. Any object smaller than 225 x 225 km can not even show up. Your coastline will have 225 km or 140 mile excursions in it. The English Channel can not exi
    st. Sicily, if you are lucky, will get one cell. To get any useful precision and accuracy, you will require hundreds of thousands of cells.

    I’ve been fascinated by reading this particular blog entry, because I’ve been nibbling at devising a scheme for linking global historical events to geography. My ideas had progressed only as far as expecting to cover the globe with hexagons, which wasn’t much of a stretch from playing the Avalon-Hill board games in which movement was based on simplified terrain overlaid with unit-hexagons. So thanks for providing the example < https://i.stack.imgur.com/L6kmP.jpg > (original posting, above). I wonder whether the 6-sided polygons would create a preöccupation with their shared edges, and discourage, or at least complicate, the use of the 8 cardinal directions. The tighly packed “soap” bubbles seem to be an improvement, in part because their variably shaped sides could encourage thinking in terms of directions from a center (e.g., of volume?), albeit of a distorted sphere.

    Alas, I too would have comparable issues with selecting a tractable geographic scale. Evaluating the hypothetical scale using “10,000 ‘grid cells’“, thus “140 miles on a side“,  I entercountered numerous issues, e.g.:
    • It’s 21 mi. across the English Channel at the Strait of Dover, and it’s approx. 60 mi. across the Channel from both Saint-Valery-sur-Somme to Hastings, and St. Alban’s Head to Cap de la Hague. So forget William the Conqueror’s Norman fleet being stuck in that French port for more than a month by unfavorable winds; his forces would just ride or march across dry land to King Harold’s Anglodanish England. And forget the anguish by Supreme Allied Commander Ike over waiting for a forecast of a necessary break in stormy weather, which would allow him a “go” for the D-Day invasion.
    • Ireland might be approximated as 300×170 mi., so it probably deserves 2 bubbles. Its loosely hourglass-shaped neighbor Wales might be boiled down to an approx. 100×150 mi. rectangle, and 1 bubble. The Irish Sea, which separates Ireland from Britain, ranges 47–120 mi. wide. Granting the Irish Sea a bubble to prevent merging Ireland’s land with Britain would require catastrophic subsidence flooding of N.W. Wales, but it’s gotta be done to keep the islands separate.
    • San Francisco, both the Peninsula and its Bay, are too small to be represented, so they’d both be combined into an uninterrupted expanse of dry land. So what the overland expedition of Juan Bautista de Anza II discovered in 1776 couldn’t be depicted. Ahhh, what the hey-ell, the U.S. National Park Sevice categorizes it as “American Latino Heritage”,  and that year’s Declaration of Independence is mucho more importante, si?.
    • Florida Current begins its flow more-or-less between the Florida Keys and Cuba, a channel known as the Straits of Florida; initially approx. 80 mi. wide and eastward, it narrows and curves into a northeastward coast-hugging flow to pass the Bahamas.
    • Cuba is approx. 140 mi. wide near its wide E. end, but narrows overall toward its W. end; e.g., it’s only approx. 30 mi. straight across the island at Havana. It would be an offense to geographic accuracy to grant both undersized geographic features: The Strait of Florida, famous as the path of Spanish Galleons to the Atlantic; and Cuba, notorious as the home to Kennedy’s defanged Bay of Pigs invasion and to Khruschëv’s risky nuke-missile installations.
    • Yucatan Channel is a 120-mi.-wide gap between the Yucatan Peninsula and the W. end of Cuba, so it’s nearly large enough to be assigned its own bubble; but locating it precisely might be relatively unimportant?
    • Isthmus of Panama is approx. 40 mi. of pestiferous jungle wide at Vasco Núñez de Balboa’s famous crossing, and 51 mi. along the Canal, so at the hypothetical scale, it would disappear, restoring the prehistoric Central American Seaway, simplifying “discovery of the Pacific” for Balboa, who could have just sailed thro’ the Central American Seaway. Pres. Teddy Roosevelt’s Gunboat Diplomacy would’ve had nothing to steal from Colombia in 1903, and Pres. Jimmy Carter would’ve been left without a subject for his foolish treaty with the Panamanian de facto dictator in 1977, and Hutchison Whampoa would have no land for operating (wink, wink, nudge, nudge) “container shipping ports” at the ends of the Canal. If as currently supposed, the prehistoric closing of the isthmus created the Gulf Stream, what happens to climate forecasting when scaling eliminates the isthmus, hmmm?
    • Crete is barely long enough: approx. 165 mi., but way too narrow: approx. 35 mi., so we’d lose the Mycenaeans and their Linear-B Greek, plus the Minoans and their Linear-A whatever. So it’s obvious that the Aegean Islands would all disappear or be artificially run aground on the W. coast of Anatolia.
    • Mainland cities that were historically significant as ports in ancient times are now sometimes inland, sometimes abandoned, notably Miletus, silted in by the Maeander River on the W. coast of Anatolia. This is an issue even for ports that don’t qualify as “ancient”, notably Pevensey (Sussex), which was the Channel port where William the Conquerer landed his invasion forces, is now 1 mi. inland. Among the issues is the alteration(s) it can cause in trade routes.
    • Bering Strait is 51 mi. at its narrowest, which normally would disqualify it from a bubble; but I dare not recreate the Bering Land Bridge thro’ stubborn inflexibility in granting bubbles.
    • Alaskan Peninsula (a.k.a. Aleutian Penin.) requires further pondering; it extends nearly 500 mi., but narrows to approx. only 50 mi. wide somewhere W. of Kodiak I.; in combination with its pendant Aleutian and Commander Is., it’s considered the N. limit to the Pacific. All those islands would disappear, so the Capt.-Cdr. Vitus Bering would die at sea instead of on his eponymous island, stubborn to his last breath. Or he and early Rossiyskaya-Amerikanskaya history would practically disappear. It seems to me that merging the newly in-name-only Bering “Sea” with the mighty Pacific and its water circulation would matter quite a lot to Arctic weather & climate. At least the more southerly combination of Admiralty, Baranof, and Chichagof Is. (Alexander Archipelago, S.E. Alaska), which includes the site of Novo-Arkhangelsk: hq. of the Russian-American Company, a place now known as Sitka: the original capital of U.S. Territorial Alaska, might be large enough to be granted a bubble or two.
    • Corsica, at approx. 50×100 mi., would disappear, which would eliminate the deaths & destruction in the wars led by its native son Napoleone di Buonaparte.
    • Biblical scholars will be pleased that the Red Sea is both long enough: 1200 mi., and wide enough at its widest: approx. 200 mi., but how accurately can its shape & direction be assigned?

    I’m either blessed or fortunate that my goal for a global model is a representation that’s more-or-less static, and isn’t the basis for iterative computation, unlike the subject hypothetical climate model. So I ”only” need to store the characteristics for each of the hypothetical “10,000 ‘grid cells’“,  and assign the historical events as they’re accumulated, not iterate over them.

    Halving the “10,000 ‘grid cells’” at “140 miles on a side” to 70 miles on a side (i.e., each 4,900 sq. mi.) would require 40,200 cells; shrinking them farther to 50 miles on a side (i.e., each 2,500 sq. mi.) would require 78,800 of them. Halving the latter cells begins to make the numbers scary: 25 miles on a side (i.e., each 625 sq. mi.) would require 315,000 of them.[#]  With their columns for lat., long., level, encoded physical characterics, and connections or transfers to adjacent cells, both in & out. Then to use them for history, during which armies or earthquakes can level cities, and emperors might rebuild them, I’d need to identify the applicable time.

    Note ** : https://chiefio.wordpress.com/2020/11/08/there-is-no-good-coordinate-system-for-climate-models/ (in original posting for that blog entry; the last substantive comment was posted on 10 November).

    Note ⛵ : Arrrgh!  Nearly every search I try for “monsoon navigation ‘south asia’” snags a page lamenting deaths of children in Monsoon floods. Wikipedia mentions only agriculture among economic impacts. The only useful & straightforwardly worded page I’ve found tells readers: “The wind blows from the northeast (towards the sea) in winter (the dry-monsoon) and from the southwest (towards the land) in summer (the wet-monsoon).”: http://en.banglapedia.org/index.php?title=Monsoon. Consider that having the wind straight behind the stern along a vessel’s centerline and straight toward its destination is not necessarily the best wind direction for a given sailing vessel. In fact, that situation creates a continual pain in the aß on a lateen-rigged vessel.

    Note # : Cell counts derived simply from E.M.’s already posted area of Earth’s surface (“196,940,000 sq mi”) as divided by the area of hypothetical square cells, approximated as planar shapes.

  20. V.P. Elect Smith says:


    I’m glad someone gets something out of my more, er, “technical” code postings.

    FWIW, near as I can tell, the “unable to post comment” thing comes when other comments have accumulated, since you last opened the page, in sufficient numbers and have been posted, such that your open page is now “stale”.

    Many times I’ve “select all” then “copy” then “reload” the page, paste back in the comment block, and away it goes!

    I’ve got no metrics for how stale is stale enough… It is more than one comment though.

    Yes, the scale issue is a giant one. As pointed out, you need on the order of 1 MILLION cells scale for anything close to reality (and likely a billion would be better for smaller scale weather issues), yet our data for temperature is on the order of 1000 to at most 6000 cells, scattered with high bias toward the USA, then Europe, then everywhere else for a few…

    Simply put: We don’t have enough historical data to START a valid model, nor enough to VALIDATE a model even if it had a million cell granularity, and the climate models we do run are in the 10K+ scale, so worthless.

    But hey, I like a challenge ;-)

    So I’m laying out a design approach that CAN scale to a million cells via one cell / core on a Beowulf like cluster.

    I have 2 “hopes” for this.

    1) Someone with 1/10000 the “Bribo / Covid Stimulus” package can use the ideas and code to do something useful.

    2) At least I can illustrate the stupidity being done now in models.

  21. V.P. Elect Smith says:

    FWIW, I’m pretty much done with the “generate cells” step now. I can generate them in any number in a quasi-close packed nearly hexagonal array.

    Next “hard bit” is figuring out how to identify all the “neighbors” of each cell. The easy ones are the E W (latitude band) neighbors as each one follows the other. The hard bit are the NE NW SE SW neighbors as since this is NOT a flat plane hex tile, but spherical, some will not always be there… So I get to figure out how to account for that and for the variable angle to neighbor…

    Once that is done, I can proceed to the rest of the stuff. (Oh, and the angle to sun, too… basically the spacial orientation of each tile.

  22. Simon Derricutt says:

    EM – from 23 December 2020 at 6:09 pm:
    The “winds across the single cell at the poles” shouldn’t really be any different than anywhere else. For any cell at the “ground” level, there may be an updraught or downdraught, and so if that cell is the centre of the low or high barometric pressure it can be having air coming into or going out to all 6 of its neighbours at the same time. In fact, if you don’t make that possible, then the air circulation won’t work in the same way as the physical world. The directions of the input or output air masses will be modified by “the Coriolis force” which is actually just conservation of linear momentum.

    It might be worth mentioning that the net mass of air entering and leaving a cell is unlikely to be zero anyway, since though I figure you’ve recognised that, it may not be as much to the front as it needs to be. The evaporated water will at the same time add to the mass of air whilst reducing the density of it. Since you’re handling the water content as separate from the non-condensing atmosphere, it really ought to drop out of the calculations without too much extra work, but the density change will be important. Also affects the local pressure inside the cell, so how much of that mass goes into other cells will depend of the pressure differences between them (and the viscosity) as well as the amount of updraught or downdraught in that cell. Updraught or downdraught will depend on the density relative to those around it. Precisely how that works when the pressures all around it are not equal (and of course they won’t be equal) doesn’t look that simple….

    The sheer size of the problem once you get down to a cell small enough to get close-enough to reality (individual clouds down to 100m or so) really shows that the larger scales actually used in the “official” GCMs haven’t a hope of getting the physics accurate and will need to have a lot of fudge built-in to what happens in a cell. Of course, you can use pattern-matching (AI) to look at the history and say, effectively, “last time the pressure map was like this, it progressed to this, so we predict that will happen this time, too”. Still, that’s the reason weather forecasting is still only accurate (enough) for 3 days or so. Plus maybe only that 1000 or so temperature measurement points. Just too sparse a grid.

  23. E.M.Smith says:


    Yeah, in general I’m hoping emergent behavior from basic physics as you described. My special concern about the pole cell is due to the night jet / polar vortex. There’s a whole lot more descending air there than anywhere else. It will be more complicated…

  24. p.g.sharrow says:

    Simon is right, establishing the actual conditions in any one cell based on observation will be daunting, A very few reporting stations in random places to “guess” the weather in the center of specific cells. You might find this software useful in the free student version.
    this is an Acad type program for land/surface research and planning by government and private agencies. And loke Acad is over $3,000 a year for the full blown version, but might give you Ideas on surface break downs.( or useful code)….pg

  25. p.g.sharrow says:

    Once you establish each of your cells center, by longitude and latitude, you will need an army of help establishing the physical conditions of that point. The above program might help with that interface. Once standard conditions for that point are determined you can then modify them with observed weather conditions from known points.
    The Bay area looks to me to be a good place to start as it is varied, well known with several reporting stations. If your model works there and then California, Mapping the planet should be simple…pg

  26. V.P. Elect Smith says:

    real topology will come at the end. Initially I’ll just load cells with a more generic physicality. Basically, rough altitude and a guess at dominant surface type from looking at a map. (Not much else you can do with a 120 nMile hexagon)

    “Someday” after fully written and running OK in a 10k test, someone with $Bucks can pay a bunch of interns to load 1000000 data points of land data… if anyone wants to run the thing…

    Or, it might be possible to scrape a maps system (like Google Maps) with an automated tool. Maybe…

    FWIW, I was figuring mostly that I’d just look at the land information in GHCN or similar and extrapolate for testing purposes (once past the “make a guess” phase…)

    I’ve also got some worries about altitude changes in landforms. Layers and all that. How to get the Rocky Mnts. to show up as something other than a 10k foot tall solid wall. Granularity just not there to do it (laterally or vertically). So if the surface of a cell is marked as “troposphere” height, it looks like a wall to the adjacent cell. But not having it there makes N. America one giant plain coast to coast…

    But that’s all programming problems for later steps. Right now I’m just trying to get facets (cells) and their orientations done ;-)

    “Writing Code, one bug at a time” 8-)

  27. p.g.sharrow says:

    might be a useful tool to gather resources, I’ve tried Gofundme. I going to examine what is available. If guys can fund themselves constructing big wooden sailboats there must be something that will work…pg

  28. V.P. Elect Smith says:


    I’ve wondered about doing something like that. Figured maybe I ought to work out a bit more than “Send me money so I can play with computers” first. So I put off the exploration of such sites until I’ve got a general framework and can show some actual work. That it isn’t a pipe dream, but a work in progress with progress…

    I could see doing that about the time I’ve got the tessellated tile world glistening in the sun, airless and dry. Then it’s “just atmosphere and oceans to go, then running it at high precision”… sort of glossing over that being the other 90% of the job ;-)

    I also didn’t want to be saying something like “I’m comfortable with FORTRAN but think I might learn some new fangleg language while I’m at it”… where saying “First 3 components written in Julia, the language being used for new Climate Models by..” then cite that lab using it (who’s name escapes me at the moment…).

    Basically I’m trying to knock off the worst “red flags” that I would call, prior to asking others for money. It’s just, IMHO, a whole lot better to say “Have been writing this in the target language for a year, and have the first 3 steps done. 3 more to go. Funding?” Or something like that.

    I also wanted to scope out my ability and skill at doing this in the target language so I could more properly “size the job”. I’d not want to be promising “Done in 6 months” when it was then discovered to be 2 years…

    So far, my skill and speed at writing Juila are rising rapidly. (Unlike Python where it was painful to get graphing right) One recent example: Adding formatted printing. Turns out it is just the Printf library you suck in then it’s the same format syntax as C printf (that I’ve already used for years). Took about 15 minutes to find that out and start doing formatted prints.

    Next stop is using arrays (that I’ve seen in a lot of examples, but not done yet). I’m going to try adding a BIG array to this code to hold: Cell#, Lat, Lon, Theta, Phi, Distance (though it may be stand alone one value not in the array), and 6 neighbor cells Cell# and vector in degrees toward it. Basically that “identify neighbors” part. Likely also the ‘inclination to pole axis” since it ought to be easy to compute at the same time.

    Once I have all that in an Array, I can then do other fun things with it in terms of building the planet. (Like save it in a file or database…)

    At that point, I’ve got a real globe of tiles centered on points, and the vector to all neighbor points, along with their relative distances and inclination angles. At that point, it’s “Here Comes The Sun!…” and I can start playing with dirt and temperatures…

    That’s my immediate big target / checkpoint. When I get there, it will be a bit of a party time. That’s likely when I’ll feel like I’m clearly “good enough” in the language and problem space (and demonstrated it…) and can properly budget / staffing plan the rest of it enough, that I could “make a decent RFC” (Request For Cash ;-)

  29. V.P. Elect Smith says:

    Finding that the algorithm from the paper is not that good at the edge cases… First it didn’t handle “one at the pole”, so I added that. Now I find out it isn’t “staggering” the rows of latitude as one would expect for hexagonal packing…

    I started to do the “find your neighbors” array packing. This is a print out of mostly that array and a few other useful items. I also cut the world down to 103 count (that gives 106 final cells). Notice that each latitude band starts at longitude 0.

    The layout of the long lines starting with “Any” (the “type” of the array) is a set of numbers for:
    CellNumber, Lat, Long, Theta, Phi, prior CellNumber (ought to always be a neighbor, either just above you when the other band ended, or just before you in a band), the cell just after you (same reasoning), then the 4 other “neighbor buckets” that I’ve stuffed a decremented number into just as filler until I finish finding them. So:

    ACells[Counter , :]=[Counter Lat Lon Theta phi Counter-1 Counter+1 Counter-3 Counter-4 Counter-5 Counter-6 ]

    I also print out distance in Theta (dTheta) and distance in Phi (dPhi) to the next cell. This stays constant through any given run of a given number of cells.

    In Theory, all I need do is “Find the dots inside about 1.5 dTheta of you, they are your neighbors”.

    The two polar dots have the following / prior 6 as their neighbors, so easy to handle that edge case.

    But notice that the Latitude Bands continue to start at 0.0 longitude. That’s an issue. So to get pseudo hex packing, I need to address that before I go on to working out the rest of the “neighbors” per cell. I can’t just take a dTheta (or dPhi that’s almost identical), bump it up to cover diagonal distances, and look for who’s that close, if the cells are really in a pseudo square packing.

    I think I now know why the author just used sqrt(area) for distance. It’s a square packing not a hex, that he got. It just didn’t stand out when he packed a small set around an empty pole at the start, as then the rest sort of skewed around that (rhomboid packing?)

    Things you learn when you get into the weeds of the code…

    ems@OdroidN2:~/Cmodel/ZZTest$ ./dobuild 
    World of 103 count cells
    dTheta 0.349 dPji 0.350 
    Any[1.0, 90.0, 0.0, 0.0, 0.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]
    dTheta 0.349 dPji 0.350 
    Any[2.0, 70.0, 0.0, 0.349066, 0.0, 1.0, 3.0, -1.0, -2.0, -3.0, -4.0]
    Any[3.0, 70.0, 60.0, 0.349066, 1.0472, 2.0, 4.0, 0.0, -1.0, -2.0, -3.0]
    Any[4.0, 70.0, 120.0, 0.349066, 2.0944, 3.0, 5.0, 1.0, 0.0, -1.0, -2.0]
    Any[5.0, 70.0, 180.0, 0.349066, 3.14159, 4.0, 6.0, 2.0, 1.0, 0.0, -1.0]
    Any[6.0, 70.0, 240.0, 0.349066, 4.18879, 5.0, 7.0, 3.0, 2.0, 1.0, 0.0]
    Any[7.0, 70.0, 300.0, 0.349066, 5.23599, 6.0, 8.0, 4.0, 3.0, 2.0, 1.0]
    dTheta 0.349 dPji 0.350 
    Any[8.0, 50.0, 0.0, 0.698132, 0.0, 7.0, 9.0, 5.0, 4.0, 3.0, 2.0]
    Any[9.0, 50.0, 30.0, 0.698132, 0.523599, 8.0, 10.0, 6.0, 5.0, 4.0, 3.0]
    Any[10.0, 50.0, 60.0, 0.698132, 1.0472, 9.0, 11.0, 7.0, 6.
    dTheta 0.349 dPji 0.350 
    Any[88.0, -50.0, 0.0, 2.44346, 0.0, 87.0, 89.0, 85.0, 84.0, 83.0, 82.0]
    Any[89.0, -50.0, 30.0, 2.44346, 0.523599, 88.0, 90.0, 86.0, 85.0, 84.0, 83.0]
    Any[90.0, -50.0, 60.0, 2.44346, 1.0472, 89.0, 91.0, 87.0, 86.0, 85.0, 84.0]
    Any[91.0, -50.0, 90.0, 2.44346, 1.5708, 90.0, 92.0, 88.0, 87.0, 86.0, 85.0]
    Any[92.0, -50.0, 120.0, 2.44346, 2.0944, 91.0, 93.0, 89.0, 88.0, 87.0, 86.0]
    Any[93.0, -50.0, 150.0, 2.44346, 2.61799, 92.0, 94.0, 90.0, 89.0, 88.0, 87.0]
    Any[94.0, -50.0, 180.0, 2.44346, 3.14159, 93.0, 95.0, 91.0, 90.0, 89.0, 88.0]
    Any[95.0, -50.0, 210.0, 2.44346, 3.66519, 94.0, 96.0, 92.0, 91.0, 90.0, 89.0]
    Any[96.0, -50.0, 240.0, 2.44346, 4.18879, 95.0, 97.0, 93.0, 92.0, 91.0, 90.0]
    Any[97.0, -50.0, 270.0, 2.44346, 4.71239, 96.0, 98.0, 94.0, 93.0, 92.0, 91.0]
    Any[98.0, -50.0, 300.0, 2.44346, 5.23599, 97.0, 99.0, 95.0, 94.0, 93.0, 92.0]
    Any[99.0, -50.0, 330.0, 2.44346, 5.75959, 98.0, 100.0, 96.0, 95.0, 94.0, 93.0]
    dTheta 0.349 dPji 0.350 
    Any[100.0, -70.0, 0.0, 2.79253, 0.0, 99.0, 101.0, 97.0, 96.0, 95.0, 94.0]
    Any[101.0, -70.0, 60.0, 2.79253, 1.0472, 100.0, 102.0, 98.0, 97.0, 96.0, 95.0]
    Any[102.0, -70.0, 120.0, 2.79253, 2.0944, 101.0, 103.0, 99.0, 98.0, 97.0, 96.0]
    Any[103.0, -70.0, 180.0, 2.79253, 3.14159, 102.0, 104.0, 100.0, 99.0, 98.0, 97.0]
    Any[104.0, -70.0, 240.0, 2.79253, 4.18879, 103.0, 105.0, 101.0, 100.0, 99.0, 98.0]
    Any[105.0, -70.0, 300.0, 2.79253, 5.23599, 104.0, 106.0, 102.0, 101.0, 100.0, 99.0]
    dTheta 0.349 dPji 0.350 
    Any[106.0, -90.0, 0.0, 3.14159, 0.0, 105.0, 104.0, 103.0, 102.0, 101.0, 100.0]

    So that’s where I’ve got to tonight…

    I think I can just offset the start point of even latitude bands by 1/2*dPhi and be good… but might need to shift the bounds on the dot placing loop so it doesn’t stop too soon…

  30. V.P. Elect Smith says:

    Well, other folks have the kitchen busy for the next few minutes, so while I wait to get started as Chef, might as well update as Chief Programmer.

    I think I’ve got the offset worked out. (Needs noodling through a QA pass…).

    On a flat tile, the next row out from the pole WILL have the same alignment (whatever alignment you use) as it is close packed in all directions. Only as you get further out does an offset show up. (The flat tiling must be looked at as rows of rings as it gets draped over a globe). It will not be ‘exactly 1/2’ offset due to that curvature, getting more ‘bent’ as you get away from the poles. So here’s what I’ve go so far:

    Interesting that the “start odd at 1/2 increment” now puts the Pole Longitude at 180. Doesn’t matter as “at the pole” all longitudes are the same. Also, row 3 starts at a 15 degree offset where row 2 at 60 degrees would imply a 30 degree offset. I need to run a very large cell run to see if that “goes away” at scale, or if I need to start the offset based on the LAST latitude band Mphi instead of the current lattitude band Mphi up to the equator, then (somehow…) use the NEXT latitude Mphi instead of the current.

    So a bit of polish needed on figuring that offset. Perhaps just use 1/2 the distance between the prior latitude band two ‘first dots’…

    ems@OdroidN2:~/Cmodel/ZZTest$ ./dobuild
    World of 180 count cells

    dTheta 0.262 dPji 0.267
    Any[1.0, 90.0, 180.0, 0.0, 3.14159, 0.0, 2.0, -2.0, -3.0, -4.0, -5.0]

    dTheta 0.262 dPji 0.267
    Any[2.0, 75.0, 0.0, 0.261799, 0.0, 1.0, 3.0, -1.0, -2.0, -3.0, -4.0]
    Any[3.0, 75.0, 60.0, 0.261799, 1.0472, 2.0, 4.0, 0.0, -1.0, -2.0, -3.0]
    Any[4.0, 75.0, 120.0, 0.261799, 2.0944, 3.0, 5.0, 1.0, 0.0, -1.0, -2.0]
    Any[5.0, 75.0, 180.0, 0.261799, 3.14159, 4.0, 6.0, 2.0, 1.0, 0.0, -1.0]
    Any[6.0, 75.0, 240.0, 0.261799, 4.18879, 5.0, 7.0, 3.0, 2.0, 1.0, 0.0]
    Any[7.0, 75.0, 300.0, 0.261799, 5.23599, 6.0, 8.0, 4.0, 3.0, 2.0, 1.0]

    dTheta 0.262 dPji 0.267
    Any[8.0, 60.0, 15.0, 0.523599, 0.261799, 7.0, 9.0, 5.0, 4.0, 3.0, 2.0]
    Any[9.0, 60.0, 45.0, 0.523599, 0.785398, 8.0, 10.0, 6.0, 5.0, 4.0, 3.0]
    Any[10.0, 60.0, 75.0, 0.523599, 1.309, 9.0, 11.0, 7.0, 6.0, 5.0, 4.0]
    Any[11.0, 60.0, 105.0, 0.523599, 1.8326, 10.0, 12.0, 8.0, 7.0, 6.0, 5.0]
    Any[12.0, 60.0, 135.0, 0.523599, 2.35619, 11.0, 13.0, 9.0, 8.0, 7.0, 6.0]
    Any[13.0, 60.0, 165.0, 0.523599, 2.87979, 12.0, 14.0, 10.0, 9.0, 8.0, 7.0]
    Any[14.0, 60.0, 195.0, 0.523599, 3.40339, 13.0, 15.0, 11.0, 10.0, 9.0, 8.0]
    Any[15.0, 60.0, 225.0, 0.523599, 3.92699, 14.0, 16.0, 12.0, 11.0, 10.0, 9.0]
    Any[16.0, 60.0, 255.0, 0.523599, 4.45059, 15.0, 17.0, 13.0, 12.0, 11.0, 10.0]
    Any[17.0, 60.0, 285.0, 0.523599, 4.97419, 16.0, 18.0, 14.0, 13.0, 12.0, 11.0]
    Any[18.0, 60.0, 315.0, 0.523599, 5.49779, 17.0, 19.0, 15.0, 14.0, 13.0, 12.0]
    Any[19.0, 60.0, 345.0, 0.523599, 6.02139, 18.0, 20.0, 16.0, 15.0, 14.0, 13.0]

    dTheta 0.262 dPji 0.267
    Any[20.0, 45.0, 0.0, 0.785398, 0.0, 19.0, 21.0, 17.0, 16.0, 15.0, 14.0]
    Any[21.0, 45.0, 21.1765, 0.785398, 0.369599, 20.0, 22.0, 18.0, 17.0, 16.0, 15.0]
    Any[22.0, 45.0, 42.3529, 0.785398, 0.739198, 21.0, 23.0, 19.0, 18.0, 17.0, 16.0]
    Any[23.0, 45.0, 63.5294, 0.785398, 1.1088, 22.0, 24.0, 20.0, 19.0, 18.0, 17.0]
    Any[144.0, -30.0, 315.0, 2.0944, 5.49779, 143.0, 145.0, 141.0, 140.0, 139.0, 138.0]
    Any[145.0, -30.0, 333.0, 2.0944, 5.81195, 144.0, 146.0, 142.0, 141.0, 140.0, 139.0]
    Any[146.0, -30.0, 351.0, 2.0944, 6.12611, 145.0, 147.0, 143.0, 142.0, 141.0, 140.0]

    dTheta 0.262 dPji 0.267
    Any[147.0, -45.0, 0.0, 2.35619, 0.0, 146.0, 148.0, 144.0, 143.0, 142.0, 141.0]
    Any[148.0, -45.0, 21.1765, 2.35619, 0.369599, 147.0, 149.0, 145.0, 144.0, 143.0, 142.0]
    Any[149.0, -45.0, 42.3529, 2.35619, 0.739198, 148.0, 150.0, 146.0, 145.0, 144.0, 143.0]
    Any[150.0, -45.0, 63.5294, 2.35619, 1.1088, 149.0, 151.0, 147.0, 146.0, 145.0, 144.0]
    Any[151.0, -45.0, 84.7059, 2.35619, 1.4784, 150.0, 152.0, 148.0, 147.0, 146.0, 145.0]
    Any[152.0, -45.0, 105.882, 2.35619, 1.848, 151.0, 153.0, 149.0, 148.0, 147.0, 146.0]
    Any[153.0, -45.0, 127.059, 2.35619, 2.21759, 152.0, 154.0, 150.0, 149.0, 148.0, 147.0]
    Any[154.0, -45.0, 148.235, 2.35619, 2.58719, 153.0, 155.0, 151.0, 150.0, 149.0, 148.0]
    Any[155.0, -45.0, 169.412, 2.35619, 2.95679, 154.0, 156.0, 152.0, 151.0, 150.0, 149.0]
    Any[156.0, -45.0, 190.588, 2.35619, 3.32639, 155.0, 157.0, 153.0, 152.0, 151.0, 150.0]
    Any[157.0, -45.0, 211.765, 2.35619, 3.69599, 156.0, 158.0, 154.0, 153.0, 152.0, 151.0]
    Any[158.0, -45.0, 232.941, 2.35619, 4.06559, 157.0, 159.0, 155.0, 154.0, 153.0, 152.0]
    Any[159.0, -45.0, 254.118, 2.35619, 4.43519, 158.0, 160.0, 156.0, 155.0, 154.0, 153.0]
    Any[160.0, -45.0, 275.294, 2.35619, 4.80479, 159.0, 161.0, 157.0, 156.0, 155.0, 154.0]
    Any[161.0, -45.0, 296.471, 2.35619, 5.17439, 160.0, 162.0, 158.0, 157.0, 156.0, 155.0]
    Any[162.0, -45.0, 317.647, 2.35619, 5.54399, 161.0, 163.0, 159.0, 158.0, 157.0, 156.0]
    Any[163.0, -45.0, 338.824, 2.35619, 5.91359, 162.0, 164.0, 160.0, 159.0, 158.0, 157.0]

    dTheta 0.262 dPji 0.267
    Any[164.0, -60.0, 15.0, 2.61799, 0.261799, 163.0, 165.0, 161.0, 160.0, 159.0, 158.0]
    Any[165.0, -60.0, 45.0, 2.61799, 0.785398, 164.0, 166.0, 162.0, 161.0, 160.0, 159.0]
    Any[166.0, -60.0, 75.0, 2.61799, 1.309, 165.0, 167.0, 163.0, 162.0, 161.0, 160.0]
    Any[167.0, -60.0, 105.0, 2.61799, 1.8326, 166.0, 168.0, 164.0, 163.0, 162.0, 161.0]
    Any[168.0, -60.0, 135.0, 2.61799, 2.35619, 167.0, 169.0, 165.0, 164.0, 163.0, 162.0]
    Any[169.0, -60.0, 165.0, 2.61799, 2.87979, 168.0, 170.0, 166.0, 165.0, 164.0, 163.0]
    Any[170.0, -60.0, 195.0, 2.61799, 3.40339, 169.0, 171.0, 167.0, 166.0, 165.0, 164.0]
    Any[171.0, -60.0, 225.0, 2.61799, 3.92699, 170.0, 172.0, 168.0, 167.0, 166.0, 165.0]
    Any[172.0, -60.0, 255.0, 2.61799, 4.45059, 171.0, 173.0, 169.0, 168.0, 167.0, 166.0]
    Any[173.0, -60.0, 285.0, 2.61799, 4.97419, 172.0, 174.0, 170.0, 169.0, 168.0, 167.0]
    Any[174.0, -60.0, 315.0, 2.61799, 5.49779, 173.0, 175.0, 171.0, 170.0, 169.0, 168.0]
    Any[175.0, -60.0, 345.0, 2.61799, 6.02139, 174.0, 176.0, 172.0, 171.0, 170.0, 169.0]

    dTheta 0.262 dPji 0.267
    Any[176.0, -75.0, 0.0, 2.87979, 0.0, 175.0, 177.0, 173.0, 172.0, 171.0, 170.0]
    Any[177.0, -75.0, 60.0, 2.87979, 1.0472, 176.0, 178.0, 174.0, 173.0, 172.0, 171.0]
    Any[178.0, -75.0, 120.0, 2.87979, 2.0944, 177.0, 179.0, 175.0, 174.0, 173.0, 172.0]
    Any[179.0, -75.0, 180.0, 2.87979, 3.14159, 178.0, 180.0, 176.0, 175.0, 174.0, 173.0]
    Any[180.0, -75.0, 240.0, 2.87979, 4.18879, 179.0, 181.0, 177.0, 176.0, 175.0, 174.0]
    Any[181.0, -75.0, 300.0, 2.87979, 5.23599, 180.0, 182.0, 178.0, 177.0, 176.0, 175.0]

    dTheta 0.262 dPji 0.267
    Any[182.0, -90.0, 0.0, 3.14159, 0.0, 181.0, 180.0, 179.0, 178.0, 177.0, 176.0]

    So some progress, but still not a really nice hex packing.

  31. p.g.sharrow says:

    I don’t know. Try as I might I can’t visualize flat representations on a round globe. I think you need to work from DATA points rather then AREAs of data. Then your areas become circles that overlap as need be and you collect data within the area to calculate the values at that areas central point…pg

  32. V.P. Elect Smith says:

    That’s basically what I’m doing. Area is only calculated as the “average size of a cell” that is then used to find “how far apart to place the dots”.

    The dots are basically the center of “soap bubbles” of air with the edges determined by the minimum energy interface that a soap bubble would make.

    What the program does is try to place those dots “equally distant”. The problem I’ve found with the Author’s work (of the original document) is that while they are equally spaced in one latitude band, and each latitude band is equally distant from the others, his method does not try to “equally space” at all between dots in two different bands. So instead of dots in a triangle with one over the space between the other two, you end up with dots directly above each other (originally seen as ‘all zero’ start cells).

    I’ve partially improved that with offsetting every even row of dots “start position”, but it has a discontinuity with the row that is on the other side, as the two are different “counts” and I’m only using the current row count in my “adjustment”. So “better” but not quite “right” (yet).

    Then, when the counts are different, there will be some dots that end up “right above each other…

    Partly this is an artifact of my changing the method to have ONE dot at the pole, that then forces the next ring to be ‘exactly 6′, which is great that close to the pole, and at some numbers of dots, the next row even becomes 12. BUT, eventually you must curve to the globe and the soap bubbles just can’t be squashed that much without removing a few… (can’t increase by 2 x at each ’round’ outward). So the dots no longer have an exact alignment with the gaps in the next row up.

    The original method has 3 dots around the pole, then arcs of dots spreading down from there. Better in the mid latitudes but at the cost of a slightly dodgy pole with no dot at all. (OTOH, wind between 3 cells around the exact center of the pole would be easier to get right in a polar vortex…)

    NO Method will be exactly right, as you can not tile a sphere with regular polygons without a discontinuity. The goal is just to get as close as possible.

    So I need to try a few more things with “one polar cell” and then with the original method of ‘no polar cell’ and pick whatever looks best.

  33. V.P. Elect Smith says:

    For anyone wondering what I’m going on about with soap bubbles:

    They form a flat line between them where they meet, and spherical otherwise, as that is the least energy film. More here:


    Double Bubble in the Plane: In 1993, Alfaro, Brock, Foisy, Hodges, and Zimba showed the configuration of shapes with the minimum perimeter enclosing two equal areas is achieved for two intersecting circles separated by a line such that the arcs all meet at 120 degree angles:

    [picture of two bubbles with flat interface between them]

    The red dots indicate the angles meeting at 120 degrees.
    The problem of finding the optimal configurations with minimum perimeter for four or more bubbles in the plane is currently still open. The conjectured configuration of shapes follows the patterns above with all intersecting circles meeting at 120 degree angles:
    [another picture…]
    In general, soap bubbles always meet in groups of threes at equal angles of 120 degrees. Showing configurations of four or bubbles are optimal is still an open problem — perhaps you will be the one to make progress on these conjectures!

    Though perhaps a honycomb over the surface of the Earth is a better analogy:

    Bee Honeycombs

    Looking at the conjectured optimal configuration for four or more bubbles in the plane leads to the following question: what happens if we consider more and more bubbles? You might notice the pattern that the conjectured optimal configuration for more bubbles begins to resemble a hexagonal tiling of the plane. If these configurations are indeed optimal, this leads to the question: do we observe hexagonal tilings in nature?

    Two examples of hexagonal tilings in nature are

    1) soap film patterns formed between two glass plates
    2) bees honeycombs, which are made of wax and are created by many bees working simultaneously in different parts of the honeycomb.
    In the 19th century, Charles Darwin observed that honeycombs were engineering feats “absolutely perfect in economising labour and wax.”

    It’s actually an interesting field of math and physics with unsolved bits in it.

    Which doesn’t make it easier for someone like me trying to do that sort of “bubble coating” of a sphere instead of a flat plate… and hoping for sort of a spherical honeycomb analog…

  34. p.g.sharrow says:

    A Geodesic hemisphere made of triangles to create polygons of pentagons and hexagons.

  35. p.g.sharrow says:

    As long as the bubble from the vertex point will reach it’s neighbors point for information from them and then modify for local observed information if available.

  36. V.P. Elect Smith says:

    I looked at geodesic spheres on the “no good grid” posting. You get 12 pentagons (no matter how many hexagons / triangles) and no straight lines for longer than between 2 pentagons.


    I’d rather have the dots aligned by latitude band (as most of the winds are…) or no polar dot, or both…

    As there is no good grid system, it is a question of choosing the one with the least crap in it and least “patching around the grid effects” needed. I think that’s either the Article Author’s original one, or some improvement on it I might make. Maybe.

    It will always be wrong, whatever the choice, but it might be more useful.

    Just realized I’m referring to the Grids posting but don’t have a link to it in this article. A year from now it will be hard for someone to figure out what I’m talking about and who the “Original Article” is talking about:

  37. V.P. Elect Smith says:

    Staring at this image for a while, it has more alignment discontinuities near the poles and fewer far away (but does have them, of necessity, the sphere not being a flat plane)

    I think I need to learn the plotting package next and make plots of the various options (instead of just tables of numbers…)

  38. Compu Gator says:

    I assume that this is the paper & author that you’ve repeatedly not cited herein:

    https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf (1 pp., 153 KB):
    How to generate equidistributed points on the surface of a sphere

    Markus Deserno
    Max-Planck-Institut für Polymerforschung [Polymer Research], Ackermannweg  10, 55128 Mainz, Germany
    Dated: September 28, 2004

    The publication is not identified in PDF file as downloaded from Carnegie-Mellon U.  Perhaps it was an author’s courtesy preprint?

Anything to say?

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.