Skip to content
Snippets Groups Projects

Feature/convergence test

Merged Feature/convergence test
feature/convergence_test into develop
Merged Cristian Lalescu requested to merge feature/convergence_test into develop
Viewing commit 6e30b66d
Show latest version
1 file
+ 281
− 0
Compare changes
  • Side-by-side
  • Inline
examples/convergence_only_time.py 0 → 100644
+ 281
− 0
 
import numpy as np
 
import h5py
 
import matplotlib.pyplot as plt
 
 
import TurTLE
 
from TurTLE import DNS, PP
 
 
base_niterations = 512
 
test_niterations = 32
 
 
def main():
 
generate_initial_conditions()
 
 
run_simulations_field()
 
plot_error_field2()
 
return None
 
 
def generate_initial_conditions():
 
# change these two values as neded.
 
# number of MPI processes to use
 
nprocesses = 8
 
# number of OpenMP threads per MPI process to use
 
nthreads_per_process = 1
 
 
# 1. Generate quasistationary state to use for initial conditions.
 
# create a dns object
 
c0 = DNS()
 
# launch the simulation
 
c0.launch([
 
'NSVE',
 
'-n', '64',
 
'--np', '{0}'.format(nprocesses),
 
'--ntpp', '{0}'.format(nthreads_per_process),
 
'--precision', 'double',
 
'--src-simname', 'B32p1e4',
 
'--src-wd', TurTLE.data_dir,
 
'--src-iteration', '0',
 
'--simname', 'base',
 
'--niter_todo', '{0}'.format(base_niterations),
 
'--niter_out', '{0}'.format(base_niterations),
 
'--overwrite-src-parameters',
 
'--kMeta', '1',
 
'--niter_stat', '16',
 
'--checkpoints_per_file', '{0}'.format(16)])
 
 
# 3. Generate initial conditions for NSE runs at 1x, 2x and 4x the
 
# base resolution.
 
# create postprocess object
 
cc = PP()
 
# compute velocity field at last iteration
 
cc.launch([
 
'get_velocity',
 
'--np', '{0}'.format(nprocesses),
 
'--ntpp', '{0}'.format(nthreads_per_process),
 
'--simname', 'base',
 
'--precision', 'double',
 
'--iter0', '{0}'.format(base_niterations),
 
'--iter1', '{0}'.format(base_niterations)])
 
# we don't need to update "checkpoint" file after "get_velocity"
 
return None
 
 
def run_simulations_field():
 
# change these two values as neded.
 
# number of MPI processes to use
 
nprocesses = 8
 
# number of OpenMP threads per MPI process to use
 
nthreads_per_process = 1
 
 
# 1. Run NSVE for the three resolutions.
 
for factor in [1, 2, 4]:
 
# create dns object
 
cc = DNS()
 
# launch simulation
 
cc.launch([
 
'NSVE',
 
'-n', '{0}'.format(64),
 
'--np', '{0}'.format(nprocesses),
 
'--ntpp', '{0}'.format(nthreads_per_process),
 
'--src-simname', 'base',
 
'--src-iteration', '{0}'.format(base_niterations),
 
'--simname', 'nsve_{0}x'.format(factor),
 
'--precision', 'double',
 
'--dtfactor', '{0}'.format(1.0 / factor),
 
'--niter_todo', '{0}'.format(test_niterations*factor),
 
'--niter_out', '{0}'.format(test_niterations*factor),
 
'--niter_stat', '{0}'.format(4*factor)])
 
 
# 2. Run NSE for the three resolutions.
 
for factor in [1, 2, 4]:
 
# create dns object
 
cc = DNS()
 
# launch simulation
 
cc.launch([
 
'NSE',
 
'-n', '{0}'.format(64),
 
'--np', '{0}'.format(nprocesses),
 
'--ntpp', '{0}'.format(nthreads_per_process),
 
'--src-simname', 'base',
 
'--src-iteration', '{0}'.format(base_niterations),
 
'--simname', 'nse_{0}x'.format(factor),
 
'--precision', 'double',
 
'--dtfactor', '{0}'.format(1.0 / factor),
 
'--niter_todo', '{0}'.format(test_niterations*factor),
 
'--niter_out', '{0}'.format(test_niterations*factor),
 
'--niter_stat', '{0}'.format(4*factor)])
 
return None
 
 
def plot_error_field2():
 
c0_list = [DNS(simname = 'nsve_{0}x'.format(factor))
 
for factor in [1, 2, 4]]
 
c1_list = [DNS(simname = 'nse_{0}x'.format(factor))
 
for factor in [1, 2, 4]]
 
cl = [c0_list, c1_list]
 
for cc in c0_list + c1_list:
 
cc.compute_statistics()
 
 
factor_list = [1, 2]
 
 
error_list = [[], []]
 
for ii in range(len(factor_list)):
 
factor = factor_list[ii]
 
for jj in [0, 1]:
 
c0 = cl[jj][ii]
 
c1 = cl[jj][ii+1]
 
df0 = h5py.File(c0.get_checkpoint_fname(iteration = test_niterations*factor), 'r')
 
df1 = h5py.File(c1.get_checkpoint_fname(iteration = test_niterations*2*factor), 'r')
 
vel0 = get_velocity(c0, iteration = test_niterations*factor)
 
vel1 = get_velocity(c1, iteration = test_niterations*2*factor)
 
diff = vel1 - vel0
 
vv = np.sum(vel0**2, axis = 3)
 
dd = np.sum(diff**2, axis = 3)
 
vrms = np.sqrt(np.mean(vv))
 
err = np.mean(dd)**0.5 / vrms
 
error_list[jj].append(err)
 
print('maximum error for factor {0} is {1}'.format(
 
factor, np.max(dd**0.5) / vrms))
 
df0.close()
 
df1.close()
 
 
f = plt.figure(figsize = (4, 4))
 
a = f.add_subplot(111)
 
fl = np.array(factor_list).astype(np.float)
 
a.plot(fl,
 
error_list[0],
 
marker = '.',
 
label = 'NSVE')
 
a.plot(fl,
 
error_list[1],
 
marker = '.',
 
label = 'NSE')
 
a.plot(fl, 1e-3 * fl**(-1),
 
dashes = (1, 1),
 
color = 'black',
 
label = '$\propto f^{-1}$')
 
a.plot(fl, 1e-3 * fl**(-2),
 
dashes = (2, 2),
 
color = 'black',
 
label = '$\propto f^{-2}$')
 
a.plot(fl, 1e-3 * fl**(-3),
 
dashes = (3, 3),
 
color = 'black',
 
label = '$\propto f^{-3}$')
 
a.set_ylabel('relative error')
 
a.set_xlabel('resolution factor $f$')
 
a.set_xscale('log')
 
a.set_yscale('log')
 
a.legend(loc = 'best')
 
f.tight_layout()
 
f.savefig('err_vs_dt_field2.pdf')
 
f.savefig('err_vs_dt_field2.svg')
 
plt.close(f)
 
return None
 
 
def plot_error_field():
 
c0_list = [DNS(simname = 'nsve_{0}x'.format(factor))
 
for factor in [1, 2, 4]]
 
c1_list = [DNS(simname = 'nse_{0}x'.format(factor))
 
for factor in [1, 2, 4]]
 
for cc in c0_list + c1_list:
 
cc.compute_statistics()
 
 
factor_list = [1, 2, 4]
 
 
error_list = []
 
for ii in range(len(factor_list)):
 
factor = factor_list[ii]
 
c0 = c0_list[ii]
 
c1 = c1_list[ii]
 
df0 = h5py.File(c0.get_checkpoint_fname(iteration = test_niterations*factor), 'r')
 
df1 = h5py.File(c1.get_checkpoint_fname(iteration = test_niterations*factor), 'r')
 
vel0 = get_velocity(c0, iteration = test_niterations*factor)
 
vel1 = get_velocity(c1, iteration = test_niterations*factor)
 
diff = vel1 - vel0
 
vv = np.sum(vel0**2, axis = 3)
 
dd = np.sum(diff**2, axis = 3)
 
vrms = np.sqrt(np.mean(vv))
 
err = np.mean(dd)**0.5 / vrms
 
error_list.append(err)
 
print('maximum error for factor {0} is {1}'.format(
 
factor, np.max(dd**0.5) / vrms))
 
df0.close()
 
df1.close()
 
 
f = plt.figure(figsize = (4, 4))
 
a = f.add_subplot(111)
 
fl = np.array(factor_list).astype(np.float)
 
a.plot(fl,
 
error_list,
 
marker = '.',
 
label = 'NSVE vs NSE')
 
a.plot(fl, 1e-2 * fl**(-1),
 
dashes = (1, 1),
 
color = 'black',
 
label = '$\propto f^{-1}$')
 
a.plot(fl, 1e-2 * fl**(-2),
 
dashes = (2, 2),
 
color = 'black',
 
label = '$\propto f^{-2}$')
 
a.plot(fl, 1e-2 * fl**(-3),
 
dashes = (3, 3),
 
color = 'black',
 
label = '$\propto f^{-3}$')
 
a.set_ylabel('relative error')
 
a.set_xlabel('resolution factor $f$')
 
a.set_xscale('log')
 
a.set_yscale('log')
 
a.legend(loc = 'best')
 
f.tight_layout()
 
f.savefig('err_vs_dt_field.pdf')
 
f.savefig('err_vs_dt_field.svg')
 
plt.close(f)
 
return None
 
 
def get_velocity(
 
cc,
 
iteration = 0):
 
data_file = h5py.File(
 
cc.get_checkpoint_fname(iteration = iteration),
 
'r')
 
if 'velocity' in data_file.keys():
 
if 'real' in data_file['velocity'].keys():
 
vel = data_file['velocity/real/{0}'.format(iteration)][()]
 
else:
 
cvel = data_file['velocity/complex/{0}'.format(iteration)][()]
 
vel = ift(cvel)
 
#vel *= cc.parameters['nx']*cc.parameters['ny']*cc.parameters['nz']
 
else:
 
assert('complex' in data_file['vorticity'].keys())
 
cvort = data_file['vorticity/complex/{0}'.format(iteration)][()]
 
cvel = compute_curl(cc, cvort, inverse_curl = True)
 
vel = ift(cvel)
 
data_file.close()
 
return vel
 
 
def ift(cfield):
 
ff = np.fft.irfftn(cfield, axes = (0, 1, 2))
 
field = np.transpose(ff, (1, 0, 2, 3)).copy()
 
return field
 
 
def compute_curl(
 
cc,
 
cfield,
 
inverse_curl = False):
 
kx = cc.get_data_file()['kspace/kx'][()]
 
ky = cc.get_data_file()['kspace/ky'][()]
 
kz = cc.get_data_file()['kspace/kz'][()]
 
ff = cfield.copy()
 
ff[..., 0] = 1j*(ky[:, None, None]*cfield[..., 2] - kz[None, :, None]*cfield[..., 1])
 
ff[..., 1] = 1j*(kz[None, :, None]*cfield[..., 0] - kx[None, None, :]*cfield[..., 2])
 
ff[..., 2] = 1j*(kx[None, None, :]*cfield[..., 1] - ky[:, None, None]*cfield[..., 0])
 
if inverse_curl:
 
k2 = (kx[None, None, :]**2 +
 
ky[ :, None, None]**2 +
 
kz[None, :, None]**2)
 
k2[0, 0, 0] = 1.
 
ff /= k2[:, :, :, None]
 
return ff
 
 
if __name__ == '__main__':
 
main()