Fun with OpenFOAM and meandering channel

This post is the result of my TA activities. So this time I've decided to show students how to make curved meshes and how to run free surface simulations with OpenFOAM.

Cause all students are from Civil and Environmental Engineering department typical curved mesh is necessary for meandering channels simulation. Surely there're another types of flumes with curved parts (rising sector weir gate for ex.) but with meandering channels we'll have more fun.

Scheme of the channel are shown on Fig. 1. We've got 270 degree meandering part with centerline radius of 1.5 m which is connected to straight 1 m inlet and outlet channels with two small curved parts of smaller radius.

Figure 1

Progress on mesh creation is shown on Figs. 2-5. First mesh without curved edges is built, then arc edges are added, then density of mesh increased, and at last upper part (where air is located) of the channel added.

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Boundary conditions are rather usual for this type of simulations. We have constant discharge at the inlet and outlet (though maybe it'll be better to set zero gradient BC at the outlet). Then it is possible to run simulation. On the video first 200 s from the start are shown.

Most of the time when comparing simulation results with experiment it is necessary to get some more of less steady flow in the flume. Usually it is achieved by running simulation for rather long time (as an estimation we can take the time necessary to refill the channel, i.e. channel length / inlet velocity). With initially still fluid in the channel it'll take rather long time to settle all the oscillations of the free surface. So it's better to get some initial estimation of the velocity distribution in the channel. For the straight channel this first estimation (constant velocity along the channel) can be set with setFields utility, in case of more complex channel shapes we can use pimpleFoam where fluid free surface do not move and actually is wall with slip BC. So in this case I took mesh from Fig. 5, ran pimpleFoam with the same inlet and outlet BCs and then used mapFields for initial flow field setup.

Video of the free surface evolution in this new simulation is shown below.

As in the previous case in the video around 200 seconds from the start of simulation are shown. And this time we get almost steady flow during this time.

OpenFOAM files for the cases can be found at https://bitbucket.org/mrklein/openfoam-course/src/9f9c41fdb851125a6a4bedb56dd3b8fd67c37330/meandering?at=master.

Fun with OpenFOAM and turbulence models

Creation of these movies was a part of sudden teaching assistance duties. In the beginning the idea was to just reproduce Kármán vortex street. And then it evolved to the comparison of vortex street simulation results with different turbulence models. Case files for the LES simulations can be found at https://bitbucket.org/mrklein/openfoam-ravings.

Channels and meshes

We have simple straight channel with two bodies of different shape put inside the channels. Simple schemes of the channels are shown on the figures below.

"Circle" channel scheme

This saved my evening ;)

Initially I've got points from OpenFOAM sample in the plane, then I collected all Xes and Ys, and then (cause there was also fluid with channel involved and velocities was not available at every point and when data is not available then it's zero)...

>> vels = importdata('vels.dat');
>> I = ~(vels(:, 3) == 0 & vels(:, 4) == 0);
>> quiver(vels(I, 1), vels(I, 2), vels(I, 3), vels(I, 4))

As all these points for zero values gone, figures look much clearer.

Strange number 123456

To exclude I/O from this post I've modified code. So now it not only calculates factorials but also sums their digits and at last sums the sums of digits and as a result there's on one I/O operation. Lisp as usual behaving quite good. And a result of this code:

	(defun factorial (n)
	(defun i_factorial (k acc)
	(if (= 0 k)
	acc
	(i_factorial (- k 1) (* k acc))))
	(i_factorial n 1))
	(defun digits (n)
	(defun i_digits (k acc)
	(if (= 0 k)
	acc
	(i_digits (truncate k 10) (cons (mod k 10) acc))))
	(i_digits n '()))
	(defun digits_s (n)
	(mapcar
	#'(lambda (k) (parse-integer (string k)))
	(coerce (write-to-string n) 'list)))
	(format T "~D~%"
	(time (reduce
	'+
	(mapcar
	#'(lambda (n)
	(reduce
	'+
	(digits_s (factorial n))))
	'(1 12 123 1234 12345 123456)))))

view raw test.lisp hosted with ❤ by GitHub

is

Evaluation took:
  7.762 seconds of real time
  7.565354 seconds of total run time (7.363660 user, 0.201694 system)
  [ Run times consist of 0.399 seconds GC time, and 7.167 seconds non-GC time. ]
  97.46% CPU
  17,812,436,719 processor cycles
  177 page faults
  15,624,347,568 bytes consed
  
2652733

Almost similar results I've got from Go program:

	package main
	import (
	"fmt"
	"math/big"
	"strings"
	"strconv"
	"time"
	)
	func Factorial(n *big.Int) *big.Int {
	var unit = big.NewInt(1)
	var acc *big.Int = big.NewInt(1)
	var t = n
	for t.Cmp(unit) == 1 {
	acc.Mul(acc, t)
	t.Sub(t, unit)
	}
	return acc
	}
	func SumDigits(n string) int64 {
	var r int64
	l := strings.Split(n, "")
	for i := 0; i < len(l); i++ {
	d, _ := strconv.ParseInt(l[i], 10, 0)
	r += d
	}
	return r
	}
	func main() {
	var s int64 = 0
	t1 := time.Now()
	n := []int64{1, 12, 123, 1234, 12345, 123456}
	for i := 0; i < len(n); i++ {
	r := Factorial(big.NewInt(n[i]))
	s += SumDigits(r.String())
	}
	t2 := time.Now()
	fmt.Printf("Elapsed time: %v\n", t2.Sub(t1))
	fmt.Println(s)
	}

view raw test.go hosted with ❤ by GitHub

Elapsed time: 4.509443s
2652733

But strangely Clojure, Haskell and Julia have got some very strange behavior. Till 12345! everything seems to be ok (though for Clojure it takes almost 1.5 sec to calculate while Lisp done it in 0.091 seconds of real time). But when we get 123456! Clojure needs almost 1 minute, Haskell and Julia just plainly refused to calculate. Even in case of Julia it was not so plainly cause it decided to eat almost all memory.

Euler #125

The palindromic number 595 is interesting because it can be written as the sum of consecutive squares: 62 + 72 + 82 + 92 + 102 + 112 + 122...

	#!/usr/bin/env python
	from math import sqrt, floor
	N = 100000000
	def is_palindrome(n):
	return str(n) == str(n)[::-1]
	def sum_list(n):
	res = set()
	for i in xrange(n - 1, 0, -1):
	l = [x*x for x in xrange(i, n + 1)]
	t = sum(l)
	if t > N:
	break
	if is_palindrome(t):
	res.add(t)
	return list(res)
	if __name__ == '__main__':
	t = int(sqrt(N/2)) + 1
	l = [item for sublist in map(sum_list, xrange(2, t + 1)) for item in sublist]
	print len(set(l)), sum(set(l))

view raw euler125.py hosted with ❤ by GitHub

What's wrong with the code?

Decided to take naive implementation of factorial in Racket and Lisp:

	#lang racket
	(define (factorial n)
	(define (i_factorial n acc)
	(if (= 0 n)
	acc
	(i_factorial (- n 1) (* n acc))))
	(i_factorial n 1))
	(map (lambda (n)
	(printf "~a~n" (factorial n))) (list 1 12 123 1234 12345))

view raw factorial.rkt hosted with ❤ by GitHub

and

	(defun factorial (n)
	(defun i_factorial (n acc)
	(if (= 0 n)
	acc
	(i_factorial (- n 1) (* acc n))))
	(i_factorial n 1))
	(defun print-factorial (n)
	(format T "~D~%" (factorial n)))
	(mapcar 'print-factorial '(1 12 123 1234 12345))

view raw factorial.lisp hosted with ❤ by GitHub

And even after

raco exe test.rkt

timing is as follows:

raco exe test.rkt
time ./test

real 0m0.697s
user 0m0.581s
sys 0m0.092s

time sbcl --script test.lisp

real 0m0.139s
user 0m0.085s
sys 0m0.049s

Maybe I/O in Racket is rather slow?

Euler #87

How many numbers below fifty million can be expressed as the sum of a prime square, prime cube, and prime fourth power?

	#!/usr/bin/env python
	from libeuler import primes
	if __name__ == '__main__':
	p = primes(7100)
	res = set()
	for s in p:
	for c in p:
	for q in p:
	n = s*s + c*c*c + q*q*q*q
	if n < 50000000:
	res.add(n)
	print len(res)

view raw euler87.py hosted with ❤ by GitHub

Choosing flash cards application

After 8 months in Korea at last I’ve decided to start learning language. I’ve got Korean language textbook which I bought as soon as I came here. With my colleague we arranged language exchange lessons. English for Korean. So the time for learning words came.
I the very beginning of my life here I bought paper flip-cards to learn alphabet. This time I’ve decided to be more technological and use iPad and iPhone. Surely there are lots of flip card sets on the Internet but they are just abstract flip cards: korean numbers, colors, verbs etc. I needed special sets of flip cards linked to the lessons in my textbook. So there appeared necessity for creation of flip cards. For this activity iPad with its onscreen keyboard was perfect. And for word learning I planned to use iPhone. And final procedure for learning became:

1. Create flip cards for the lesson on iPad
2. Export it somewhere
3. Import them on iPhone
4. Learn words

Euler #64

All square roots are periodic when written as continued fractions and can be written in the form...

Method for calculation of square root continued fraction is taken from  http://web.math.princeton.edu/mathlab/jr02fall/Periodicity/mariusjp.pdf

	#!/usr/bin/env python
	def period(n):
	a_0 = int(n ** 0.5)
	if a_0 * a_0 == n:
	return 0
	b, b_0 = a_0, a_0
	c, c_0 = n - a_0*a_0, n - a_0*a_0
	res = 0
	while True:
	a = (a_0 + b) / c
	b = a * c - b
	c = (n - b * b) / c
	res += 1
	if b == b_0 and c == c_0:
	break
	return res
	if __name__ == '__main__':
	print sum(
	map(lambda n: period(n) % 2,
	xrange(2, 10001)))

view raw 64.py hosted with ❤ by GitHub

Brute-forcing Eurler #70

Euler's Totient function, φ(n) [sometimes called the phi function], is used to determine the number of positive numbers less than or equal to n which are relatively prime to n. For example, as 1, 2, 4, 5, 7, and 8, are all less than nine and relatively prime to nine, φ(9)=6....

First I was trying to  brute-force the solution using FT of Euler totient function (http://en.wikipedia.org/wiki/Euler%27s_totient_function#Fourier_transform), than through the definition but all these solutions was too long to wait. In the end finished with product formula.

In the process of solution found quite faster recipe from ActiveState on generating prime numbers (http://code.activestate.com/recipes/366178-a-fast-prime-number-list-generator/).

Also used multiprocessing module to speed up solution.

	#!/usr/bin/python
	from operator import mul
	from multiprocessing import Pool, Manager
	manager = Manager()
	primes_cache = manager.list()
	def primes(n):
	if n == 2:
	return [2]
	elif n < 2:
	return []
	s = range(3, n + 1, 2)
	mroot = n ** 0.5
	half = (n + 1) / 2 - 1
	i = 0
	m = 3
	while m <= mroot:
	if s[i]:
	j = (m * m - 3) / 2
	s[j] = 0
	while j < half:
	s[j] = 0
	j += m
	i = i + 1
	m = 2 * i + 3
	return [2] + [x for x in s if x]
	primes_cache = primes(10000000)
	def factorize(n):
	for prime in primes_cache:
	if prime > n:
	return
	exponent = 0
	while n % prime == 0:
	exponent, n = exponent + 1, n / prime
	if exponent != 0:
	yield prime, exponent
	def totient(n):
	return reduce(
	mul,
	((p - 1) * p ** (exp - 1)
	for p, exp in factorize(n)),
	1)
	def euler70(n):
	t = totient(n)
	if sorted(str(n)) == sorted(str(t)):
	return (n, t, float(n) / t)
	else:
	return (n, t, 1e7)
	if __name__ == '__main__':
	pool = Pool(processes=8)
	print sorted(
	pool.map(
	euler70,
	xrange(2, 10000000)),
	key=lambda n: n[2])[0]

view raw 70.py hosted with ❤ by GitHub

And after...

$ time ./70.py 
(8319823, 8313928, 1.0007090511248113)

real    235m23.873s
user    1644m0.217s
sys     0m36.206s

Euler #80

It is well known that if the square root of a natural number is not an integer, then it is irrational.

from decimal import *
from math import sqrt

getcontext().prec = 105

reduce(int.__add__, 
       map(lambda a: \
           reduce(lambda x,y: int(x)+int(y), 
                  str(a.sqrt()).replace('.', '')[:100], 
                  0), 
           map(lambda n: Decimal(str(n)), 
               map(lambda k: \
                   0 if int(sqrt(k)) ** 2 == k else k, 
                   range(100))
              )
          ), 
       0)

Surely it's cheating ;)