Skip to content
Snippets Groups Projects
Commit 6e30b66d authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

adds temporary alternate field convergence test

parent 696e7669
No related branches found
No related tags found
No related merge requests found
This commit is part of merge request !35. Comments created here will be created in the context of that merge request.
Show whitespace changes
Inline Side-by-side
examples/convergence_only_time.py 0 → 100644
+ 281
− 0
View file @ 6e30b66d
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()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment