diff --git a/.gitmodules b/.gitmodules index 391503e..e69de29 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +0,0 @@ -[submodule "JYU-LB"] - path = JYU-LB - url = https://github.com/simphony/JYU-LB diff --git a/.travis.yml b/.travis.yml index 5bb2bc0..c8a1d8a 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,20 +1,24 @@ language: python env: - - SIMPHONY_VERSION=0.4.0 + - SIMPHONY_VERSION=0.6.0 - SIMPHONY_VERSION=master matrix: allow_failures: - env: SIMPHONY_VERSION=master install: - pip install -e git+https://github.com/simphony/simphony-common.git@${SIMPHONY_VERSION}#egg=simphony + - pip install Cython + - pip install -e git+https://github.com/simphony/JYU-LB.git#egg=jyulb - pip install -r dev_requirements.txt - - sudo ./install_external.sh +# - sudo ./install_external.sh - pip install coveralls - python setup.py develop script: - - flake8 . - - coverage run -m unittest discover + - flake8 wrapper/ + - flake8 examples/ +# - flake8 . + - coverage run -m unittest discover -s wrapper - coverage report after_success: diff --git a/JYU-LB b/JYU-LB deleted file mode 160000 index c7d07ba..0000000 --- a/JYU-LB +++ /dev/null @@ -1 +0,0 @@ -Subproject commit c7d07ba1666966f7c1cf5a8202d078240cf43ab4 diff --git a/dev_requirements.txt b/dev_requirements.txt index 4518fe6..0a42a47 100644 --- a/dev_requirements.txt +++ b/dev_requirements.txt @@ -1,4 +1,4 @@ -numpy==1.9.1 +numpy>=1.9.1 +cython>=0.25.2 coverage flake8 -cython diff --git a/examples/couette.py b/examples/couette.py new file mode 100644 index 0000000..3c953a1 --- /dev/null +++ b/examples/couette.py @@ -0,0 +1,226 @@ +"""Coeutte fluid flow simulation (flow between two moving, parallel plates). + +Details +------- +- 3D simulation (with two "dummy" direction) +- flow driven by moving walls +- periodic bc in the "dummy" directions +""" +import time +import numpy as np +from simphony.cuds.meta import api +from simphony.core.cuba import CUBA +from simphony.api import CUDS, Simulation +from simphony.engine import EngineInterface +from simphony.cuds.lattice import make_cubic_lattice +from jyulb.flow_field import couette_plate_steady_state + +# CUDS +cuds = CUDS(name='couette') + +# Physics model +cfd = api.Cfd(name='fluid flow') + +# Submodels (in fact, these are the default submodels in CFD) +cfd.thermal_model = api.IsothermalModel(name='isothermal') +cfd.turbulence_model = api.LaminarFlowModel(name='laminar') +cfd.rheology_model = api.NewtonianFluidModel(name='newtonian') +cfd.multiphase_model = api.SinglePhaseModel(name='singlephase') +cfd.compressibility_model = api.IncompressibleFluidModel(name='incompressible') +cuds.add([cfd]) + +# Materials +fluid = api.Material(name='fluid') +solid = api.Material(name='wall') +fluid.data[CUBA.DENSITY] = 1.0 +fluid.data[CUBA.KINEMATIC_VISCOSITY] = 1.0/6.0 +cuds.add([fluid, solid]) + +# Dataset (lattice) +lat_size = (10, 1, 10) # number of lattice nodes +lat_h = 1.0 # lattice spacing +lat = make_cubic_lattice('channel', lat_h, lat_size) + +# Create a channel geometry, set the initial flow field +lower_wall_vel1 = -1.0e-3 +lower_wall_vel2 = 0.0e-3 +upper_wall_vel1 = 1.0e-3 +upper_wall_vel2 = -2.0e-3 + +for node in lat.iter(): + ijk = node.index + # The two plates are separated in the z-direction + if ijk[2] == 0: + node.data[CUBA.MATERIAL] = solid + node.data[CUBA.VELOCITY] = (lower_wall_vel1, lower_wall_vel2, 0) + elif ijk[2] == (lat_size[2]-1): + node.data[CUBA.MATERIAL] = solid + node.data[CUBA.VELOCITY] = (upper_wall_vel1, upper_wall_vel2, 0) + else: + node.data[CUBA.MATERIAL] = fluid + node.data[CUBA.VELOCITY] = (0, 0, 0) + + lat.update([node]) + +cuds.add([lat]) + +# Boundary conditions +# for the exterior boundary faces of the simulation domain +# (the directions refer to unit outward normals) +face1y = api.Boundary(name='face1y', condition=[api.Periodic()]) +face1y.data[CUBA.DIRECTION] = (0, -1, 0) + +facey1 = api.Boundary(name='facey1', condition=[api.Periodic()]) +facey1.data[CUBA.DIRECTION] = (0, 1, 0) + +face1z = api.Boundary(name='face1z', condition=[api.Periodic()]) +face1z.data[CUBA.DIRECTION] = (0, 0, -1) + +facez1 = api.Boundary(name='facez1', condition=[api.Periodic()]) +facez1.data[CUBA.DIRECTION] = (0, 0, 1) + +cuds.add([face1y, facey1, face1z, facez1]) + +# Solver parameters (time integration, collision operator) +ti = api.IntegrationTime(name='simulation time', + current=0.0, final=1000.0, size=1.0) + +ti.data[CUBA.COLLISION_OPERATOR] = 'TRT' +cuds.add([ti]) + +# Simulate +print '-'*77 +sim = Simulation(cuds, 'JYU-LB', engine_interface=EngineInterface.Internal) +lat = cuds.get_by_name(lat.name) + +# Total fluid mass +tot_fmass = 0.0 +for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_fmass += node.data[CUBA.DENSITY] + +print '='*77 +print 'Couette flow simulation' +print '='*77 +print "Total fluid mass before simulation: {0:.4e}".format(tot_fmass) +print '-'*77 + +start = time.time() + +for e in range(10): + sim.run() + + # Total velocity + tot_vx = 0.0 + tot_vy = 0.0 + tot_vz = 0.0 + for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_vx += node.data[CUBA.VELOCITY][0] + tot_vy += node.data[CUBA.VELOCITY][1] + tot_vz += node.data[CUBA.VELOCITY][2] + + sf = 'Time = {0:f}, tot.vel = ({1:11.4e}, {2:11.4e}, {3:11.4e})' + print sf.format(ti.current, tot_vx, tot_vy, tot_vz) + +end = time.time() + +# Total fluid mass +tot_fmass = 0.0 +for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_fmass += node.data[CUBA.DENSITY] + +print '-'*77 +print "Total fluid mass after simulation: {0:.4e}".format(tot_fmass) +print '-'*77 +print "Time spend in run (s): ", end-start +print '-'*77 + +# --------------------------------------------------------------------------- +# Post-processing: +# Analytical solution for the Couette flow, error, and output +# --------------------------------------------------------------------------- +# GLE: flow speed in y-data format +# --------------------------------------------------------------------------- +H = lat_size[2]-2 # plates separated in the z-direction +centrel_dist = np.zeros((H), dtype=np.float64) +ana_vel = np.zeros((H, 2), dtype=np.float64) + +couette_plate_steady_state(H, lower_wall_vel1, lower_wall_vel2, + upper_wall_vel1, upper_wall_vel2, + centrel_dist, ana_vel) + +sim_loc_vel = np.zeros(3, dtype=np.float64) +coord = np.zeros(3, dtype=np.int32) +coord[0] = int(lat_size[0]/2) +coord[1] = int(lat_size[1]/2) +coord[2] = int(lat_size[2]/2) + +ana2 = 0.0 +sim_ana_diff2 = 0.0 +sim_loc_vel = np.zeros((3), dtype=np.float64) + +with open('flow_prof.txt', 'w') as f: + for h in range(H): + coord[2] = h+1 + node = lat.get(coord) + dist = centrel_dist[h] + + if node.data[CUBA.MATERIAL] == solid: + continue + + ana_vel1 = ana_vel[h, 0] + ana_vel2 = ana_vel[h, 1] + + sim_loc_p = node.data[CUBA.PRESSURE] + sim_loc_den = node.data[CUBA.DENSITY] + sim_loc_vel[0] = node.data[CUBA.VELOCITY][0] + sim_loc_vel[1] = node.data[CUBA.VELOCITY][1] + sim_loc_vel[2] = node.data[CUBA.VELOCITY][2] + + sim_loc_speed = np.sqrt(sim_loc_vel[0]*sim_loc_vel[0] + + sim_loc_vel[1]*sim_loc_vel[1] + + sim_loc_vel[2]*sim_loc_vel[2]) + + ana_loc_speed = np.sqrt(ana_vel1*ana_vel1 + ana_vel2*ana_vel2) + + sim_ana_diff = sim_loc_speed - ana_loc_speed + sim_ana_diff2 += sim_ana_diff*sim_ana_diff + ana2 += ana_loc_speed*ana_loc_speed + + f.write('{0:f} {1:e} {2:e} {3:e} {4:e} {5:e} {6:e}\n'.format(dist, + sim_loc_den, sim_loc_vel[0], sim_loc_vel[1], sim_loc_vel[2], + sim_loc_speed, ana_loc_speed)) + +rel_l2_err_norm_vel = np.sqrt(sim_ana_diff2/ana2) + +print 'Relative L2-error norm: {0:.4e}'.format(rel_l2_err_norm_vel) +print '-'*77 + +# --------------------------------------------------------------------------- +# Matlab: 2D flow field +# --------------------------------------------------------------------------- +coord[0] = int(lat_size[0]/2) +coord[1] = int(lat_size[1]/2) +coord[2] = int(lat_size[2]/2) + +with open('flow_2D.field', 'w') as f: + for x1 in range(lat_size[0]): + for x2 in range(lat_size[2]): + coord[0] = x1 + coord[2] = x2 + + node = lat.get(coord) + if node.data[CUBA.MATERIAL] == solid: + continue + + p = node.data[CUBA.PRESSURE] + dn = node.data[CUBA.DENSITY] + v1 = node.data[CUBA.VELOCITY][0] + v2 = node.data[CUBA.VELOCITY][1] + + sf = '{0:d} {1:d} {2:e} {3:e} {4:e}\n' + f.write(sf.format(x1, x2, dn, v1, v2)) + +# --------------------------------------------------------------------------- diff --git a/examples/lid_driven_cavity.py b/examples/lid_driven_cavity.py new file mode 100644 index 0000000..43bb197 --- /dev/null +++ b/examples/lid_driven_cavity.py @@ -0,0 +1,169 @@ +"""Lid-driven cavity flow (2D setup). + +Details +------- +- 3D simulation (2D setup, i.e. periodic in one direction) +- flow driven by a moving wall +""" +import time +import numpy as np +from simphony.cuds.meta import api +from simphony.core.cuba import CUBA +from simphony.api import CUDS, Simulation +from simphony.engine import EngineInterface +from simphony.cuds.lattice import make_cubic_lattice + +# CUDS +cuds = CUDS(name='lid-driven cavity') + +# Physics model +cfd = api.Cfd(name='fluid flow') + +# Submodels (in fact, these are the default submodels in CFD) +cfd.thermal_model = api.IsothermalModel(name='isothermal') +cfd.turbulence_model = api.LaminarFlowModel(name='laminar') +cfd.rheology_model = api.NewtonianFluidModel(name='newtonian') +cfd.multiphase_model = api.SinglePhaseModel(name='singlephase') +cfd.compressibility_model = api.IncompressibleFluidModel(name='incompressible') +cuds.add([cfd]) + +# Materials +fluid = api.Material(name='fluid') +solid = api.Material(name='wall') +fluid.data[CUBA.DENSITY] = 1.0 +fluid.data[CUBA.KINEMATIC_VISCOSITY] = 1.0/120.0 +cuds.add([fluid, solid]) + +# Dataset (lattice) +H = 100 # Distance between the lower (not moving) and upper (moving) wall +W = 100 # Distance between the side walls (not moving) +lat_size = (W+2, 1, H+2) # number of lattice nodes +lat_h = 1.0 # lattice spacing +lat = make_cubic_lattice('cavity', lat_h, lat_size) + +# Create a cavity geometry, set the initial flow field +# The moving wall is located at the upper exterior boundary face (z-dir.) +# No walls in the "dummy" direction (y-dir.) +mwall_vel = 5.0e-2 +for node in lat.iter(): + ijk = node.index + + if ijk[2] == (lat_size[2]-1): + node.data[CUBA.MATERIAL] = solid + node.data[CUBA.VELOCITY] = (mwall_vel, 0, 0) + elif ijk[0] == 0 or ijk[0] == (lat_size[0]-1) or ijk[2] == 0: + node.data[CUBA.MATERIAL] = solid + node.data[CUBA.VELOCITY] = (0, 0, 0) + else: + node.data[CUBA.MATERIAL] = fluid + node.data[CUBA.VELOCITY] = (0, 0, 0) + + lat.update([node]) + +cuds.add([lat]) + +# Reynolds number +Re = H*mwall_vel/fluid.data[CUBA.KINEMATIC_VISCOSITY] + +# Boundary conditions +# for the exterior boundary faces of the simulation domain +# (the directions refer to unit outward normals) +face1y = api.Boundary(name='face1y', condition=[api.Periodic()]) +face1y.data[CUBA.DIRECTION] = (0, -1, 0) + +facey1 = api.Boundary(name='facey1', condition=[api.Periodic()]) +facey1.data[CUBA.DIRECTION] = (0, 1, 0) + +face1z = api.Boundary(name='face1z', condition=[api.Periodic()]) +face1z.data[CUBA.DIRECTION] = (0, 0, -1) + +facez1 = api.Boundary(name='facez1', condition=[api.Periodic()]) +facez1.data[CUBA.DIRECTION] = (0, 0, 1) + +cuds.add([face1y, facey1, face1z, facez1]) + +# Solver parameters (time integration, collision operator) +ti = api.IntegrationTime(name='simulation time', + current=0.0, final=5000.0, size=1.0) + +ti.data[CUBA.COLLISION_OPERATOR] = 'Regularization' +cuds.add([ti]) + +# Simulate +print '-'*77 +sim = Simulation(cuds, 'JYU-LB', engine_interface=EngineInterface.Internal) +lat = cuds.get_by_name(lat.name) + +# Total fluid mass +tot_fmass = 0.0 +for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_fmass += node.data[CUBA.DENSITY] + +print '-'*77 +print 'Lid-Driven Cavity flow simulation (Re = {0:f})'.format(Re) +print '-'*77 +print "Total fluid mass before simulation: {0:.4e}".format(tot_fmass) +print '-'*77 + +start = time.time() + +for e in range(25): + sim.run() + + # Total velocity + tot_vx = 0.0 + tot_vy = 0.0 + tot_vz = 0.0 + for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_vx += node.data[CUBA.VELOCITY][0] + tot_vy += node.data[CUBA.VELOCITY][1] + tot_vz += node.data[CUBA.VELOCITY][2] + + sf = 'Time = {0:f}, tot.vel = ({1:11.4e}, {2:11.4e}, {3:11.4e})' + print sf.format(ti.current, tot_vx, tot_vy, tot_vz) + +end = time.time() + +# Total fluid mass +tot_fmass = 0.0 +for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_fmass += node.data[CUBA.DENSITY] + +print '-'*77 +print "Total fluid mass after simulation: {0:.4e}".format(tot_fmass) +print '-'*77 +print "Time spend in run (s): ", end-start +print '-'*77 + +# --------------------------------------------------------------------------- +# Post-processing: output +# --------------------------------------------------------------------------- +# Matlab: 2D flow field +# --------------------------------------------------------------------------- +coord = np.zeros(3, dtype=np.int32) +coord[0] = int(lat_size[0]/2) +coord[1] = int(lat_size[1]/2) +coord[2] = int(lat_size[2]/2) + +with open('flow_2D.field', 'w') as f: + for x1 in range(lat_size[0]): + for x2 in range(lat_size[2]): + coord[0] = x1 + coord[2] = x2 + + node = lat.get(coord) + if node.data[CUBA.MATERIAL] == solid: + continue + + p = node.data[CUBA.PRESSURE] + dn = node.data[CUBA.DENSITY] + v1 = node.data[CUBA.VELOCITY][0] + v2 = node.data[CUBA.VELOCITY][2] + + sf = '{0:d} {1:d} {2:e} {3:e} {4:e}\n' + f.write(sf.format(x1, x2, dn, v1, v2)) + +# --------------------------------------------------------------------------- diff --git a/examples/plot_field_2D.m b/examples/plot_field_2D.m new file mode 100644 index 0000000..c11ff30 --- /dev/null +++ b/examples/plot_field_2D.m @@ -0,0 +1,109 @@ +function plot_field_2D(fname) + +field2D = load(fname); +ncnt = length(field2D); + +% Assume integer coordinates +min_i = min(field2D(:,1)); +min_j = min(field2D(:,2)); +max_i = max(field2D(:,1)); +max_j = max(field2D(:,2)); + +ni = max_i - min_i + 1 +nj = max_j - min_j + 1 + +% Allocate 2D arrays for the flow field +is = zeros(nj,ni); +js = zeros(nj,ni); +dns = zeros(nj,ni); +uxs = zeros(nj,ni); +uys = zeros(nj,ni); +us = zeros(nj,ni); +vorts = zeros(nj,ni); +ekins = zeros(nj,ni); +%is = zeros(ni,nj); +%js = zeros(ni,nj); +%dns = zeros(ni,nj); +%uxs = zeros(ni,nj); +%uys = zeros(ni,nj); +%us = zeros(ni,nj); +%vorts = zeros(ni,nj); +%ekins = zeros(ni,nj); + +% Read and store flow field into 2D arrays +for n=1:ncnt + i = field2D(n,1) - min_i + 1; + j = field2D(n,2) - min_j + 1; + dn = field2D(n,3); + ux = field2D(n,4); + uy = field2D(n,5); + u2 = ux^2 + uy^2; + u = sqrt(u2); + ekin = 0.5*dn*u2; + + is(j,i) = i; + js(j,i) = j; + dns(j,i) = dn; + uxs(j,i) = ux; + uys(j,i) = uy; + us(j,i) = u; + ekins(j,i) = ekin; +end + +% Compute vorticity +for i=2:(ni-1) + for j=2:(nj-1) + dx_uy = 0.5*(uys(j,i+1) - uys(j,i-1)); + dy_ux = 0.5*(uxs(j+1,i) - uxs(j-1,i)); + vorts(j,i) = dx_uy - dy_ux; + end + end + +figure; +surface(dns); +%contourf(dns); +title('Simulated den'); + +figure; +surface(uxs); +title('Simulated ux'); + +figure; +contourf(uxs); +title('Simulated ux'); + +figure; +surface(uys); +title('Simulated uy'); + +figure; +contourf(uys); +title('Simulated uy'); + +figure; +%surface(us); +contourf(us); +title('Simulated flow speed'); + +%figure; +%surface(ekins); +%contourf(ekins); +%title('Simulated kinetic energy'); + +%figure; +%surface(vorts); +%contourf(vorts); +%title('Simulated vorticity'); + +figure; +title('Simulated velocity'); +%step = 2 +%[X,Y] = meshgrid(1:step:nj, 1:step:ni); +%quiver(X,Y,uxs(1:step:nj,1:step:ni),uys(1:step:nj,1:step:ni),6); +quiver(uxs,uys,6); + +figure; +title('Simulated velocity'); +%[X,Y] = meshgrid(1:(nj/5):nj, 1:(ni/5):ni); +%streamline(uxs,uys,X,Y); +streamslice(uxs,uys,2.5,'cubic'); diff --git a/examples/plot_flow_prof.gle b/examples/plot_flow_prof.gle new file mode 100644 index 0000000..fdfda5f --- /dev/null +++ b/examples/plot_flow_prof.gle @@ -0,0 +1,22 @@ +include "graphutil.gle" + +size 12 12 + +begin graph + size 12 12 + nobox + +! yaxis min -2E-05 max 3E-05 +! xaxis min -1E-02 max 1E-02 + + title "Flow profile" hei 0.37 + xtitle "Distance from the centerline" hei 0.37 + ytitle "Flow speed" hei 0.37 + + data flow_prof.txt d1=c1,c6 + data flow_prof.txt d2=c1,c7 + + key pos bc compact offset 0.0 1.5 hei 0.35 + d1 line lwidth 0.03 marker square msize 0.35 color black key "Simulated" + d2 line lwidth 0.03 marker circle msize 0.31 color red key "Analytical" +end graph diff --git a/examples/poiseuille_grav.py b/examples/poiseuille_grav.py new file mode 100644 index 0000000..1c3b533 --- /dev/null +++ b/examples/poiseuille_grav.py @@ -0,0 +1,219 @@ +"""Poiseuille fluid flow simulation (flow between two parallel plates). + +Details +------- +- 3D simulation (with one "dummy" direction) +- flow driven by gravity (effective pressure gradient) +- periodic bc in the flow direction (as well as in the "dummy" direction) +""" +import time +import numpy as np +from simphony.cuds.meta import api +from simphony.core.cuba import CUBA +from simphony.api import CUDS, Simulation +from simphony.engine import EngineInterface +from simphony.cuds.lattice import make_cubic_lattice +from jyulb.flow_field import poise_plate_steady_state_pgrad + +# CUDS +cuds = CUDS(name='poiseuille') + +# Physics model +cfd = api.Cfd(name='fluid flow') + +# Submodels (in fact, these are the default submodels in CFD) +cfd.thermal_model = api.IsothermalModel(name='isothermal') +cfd.turbulence_model = api.LaminarFlowModel(name='laminar') +cfd.rheology_model = api.NewtonianFluidModel(name='newtonian') +cfd.multiphase_model = api.SinglePhaseModel(name='singlephase') +cfd.compressibility_model = api.IncompressibleFluidModel(name='incompressible') + +# Gravity in the z-direction +cfd.gravity_model.acceleration = (0, 0, 1e-5) +cuds.add([cfd]) + +# Materials +fluid = api.Material(name='fluid') +solid = api.Material(name='wall') +fluid.data[CUBA.DENSITY] = 1.0 +fluid.data[CUBA.KINEMATIC_VISCOSITY] = 1.0/6.0 +cuds.add([fluid, solid]) + +# Dataset (lattice) +lat_size = (10, 1, 10) # number of lattice nodes +lat_h = 1.0 # lattice spacing +lat = make_cubic_lattice('channel', lat_h, lat_size) + +# Create a channel geometry, set the initial flow field +for node in lat.iter(): + ijk = node.index + # The two plates are separated in the x-direction + if ijk[0] == 0 or ijk[0] == (lat_size[0]-1): + node.data[CUBA.MATERIAL] = solid + node.data[CUBA.VELOCITY] = (0, 0, 0) + else: + node.data[CUBA.MATERIAL] = fluid + node.data[CUBA.VELOCITY] = (0, 0, 0) + + lat.update([node]) + +cuds.add([lat]) + +# Boundary conditions +# for the exterior boundary faces of the simulation domain +# (the directions refer to unit outward normals) +face1y = api.Boundary(name='face1y', condition=[api.Periodic()]) +face1y.data[CUBA.DIRECTION] = (0, -1, 0) + +facey1 = api.Boundary(name='facey1', condition=[api.Periodic()]) +facey1.data[CUBA.DIRECTION] = (0, 1, 0) + +face1z = api.Boundary(name='face1z', condition=[api.Periodic()]) +face1z.data[CUBA.DIRECTION] = (0, 0, -1) + +facez1 = api.Boundary(name='facez1', condition=[api.Periodic()]) +facez1.data[CUBA.DIRECTION] = (0, 0, 1) + +cuds.add([face1y, facey1, face1z, facez1]) + +# Solver parameters (time integration, collision operator) +ti = api.IntegrationTime(name='simulation time', + current=0.0, final=1000.0, size=1.0) + +ti.data[CUBA.COLLISION_OPERATOR] = 'TRT' +cuds.add([ti]) + +# Simulate +print '-'*77 +sim = Simulation(cuds, 'JYU-LB', engine_interface=EngineInterface.Internal) +lat = cuds.get_by_name(lat.name) + +# Total fluid mass +tot_fmass = 0.0 +for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_fmass += node.data[CUBA.DENSITY] + +print '='*77 +print 'Poiseuille flow simulation' +print '='*77 +print "Total fluid mass before simulation: {0:.4e}".format(tot_fmass) +print '-'*77 + +start = time.time() + +for e in range(10): + sim.run() + + # Total velocity + tot_vx = 0.0 + tot_vy = 0.0 + tot_vz = 0.0 + for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_vx += node.data[CUBA.VELOCITY][0] + tot_vy += node.data[CUBA.VELOCITY][1] + tot_vz += node.data[CUBA.VELOCITY][2] + + sf = 'Time = {0:f}, tot.vel = ({1:11.4e}, {2:11.4e}, {3:11.4e})' + print sf.format(ti.current, tot_vx, tot_vy, tot_vz) + +end = time.time() + +# Total fluid mass +tot_fmass = 0.0 +for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_fmass += node.data[CUBA.DENSITY] + +print '-'*77 +print "Total fluid mass after simulation: {0:.4e}".format(tot_fmass) +print '-'*77 +print "Time spend in run (s): ", end-start +print '-'*77 + +# --------------------------------------------------------------------------- +# Post-processing: +# Analytical solution for the Poiseuille flow, error, and output +# --------------------------------------------------------------------------- +# GLE: flow speed in y-data format +# --------------------------------------------------------------------------- +H = lat_size[0]-2 # plates separated in the x-direction +dvisc = fluid.data[CUBA.DENSITY]*fluid.data[CUBA.KINEMATIC_VISCOSITY] +eff_pgrad = -fluid.data[CUBA.DENSITY]*cfd.gravity_model.acceleration[2] + +centrel_dist = np.zeros((H), dtype=np.float64) +ana_vel = np.zeros((H), dtype=np.float64) + +poise_plate_steady_state_pgrad(H, dvisc, eff_pgrad, centrel_dist, ana_vel) + +sim_loc_vel = np.zeros(3, dtype=np.float64) +coord = np.zeros(3, dtype=np.int32) +coord[0] = int(lat_size[0]/2) +coord[1] = int(lat_size[1]/2) +coord[2] = int(lat_size[2]/2) + +ana2 = 0.0 +sim_ana_diff2 = 0.0 +sim_loc_vel = np.zeros((3), dtype=np.float64) + +with open('flow_prof.txt', 'w') as f: + for h in range(H): + coord[0] = h+1 + node = lat.get(coord) + dist = centrel_dist[h] + + if node.data[CUBA.MATERIAL] == solid: + continue + + sim_loc_p = node.data[CUBA.PRESSURE] + sim_loc_den = node.data[CUBA.DENSITY] + sim_loc_vel[0] = node.data[CUBA.VELOCITY][0] + sim_loc_vel[1] = node.data[CUBA.VELOCITY][1] + sim_loc_vel[2] = node.data[CUBA.VELOCITY][2] + + sim_loc_speed = np.sqrt(sim_loc_vel[0]*sim_loc_vel[0] + + sim_loc_vel[1]*sim_loc_vel[1] + + sim_loc_vel[2]*sim_loc_vel[2]) + + ana_loc_vel = ana_vel[h] + + sim_ana_diff = sim_loc_vel[2] - ana_loc_vel + sim_ana_diff2 += sim_ana_diff*sim_ana_diff + ana2 += ana_loc_vel*ana_loc_vel + + f.write('{0:f} {1:e} {2:e} {3:e} {4:e} {5:e} {6:e}\n'.format(dist, + sim_loc_den, sim_loc_vel[0], sim_loc_vel[1], sim_loc_vel[2], + sim_loc_speed, ana_loc_vel)) + +rel_l2_err_norm_vel = np.sqrt(sim_ana_diff2/ana2) + +print 'Relative L2-error norm: {0:.4e}'.format(rel_l2_err_norm_vel) +print '-'*77 + +# --------------------------------------------------------------------------- +# Matlab: 2D flow field +# --------------------------------------------------------------------------- +coord[0] = int(lat_size[0]/2) +coord[1] = int(lat_size[1]/2) +coord[2] = int(lat_size[2]/2) + +with open('flow_2D.field', 'w') as f: + for x1 in range(lat_size[2]): + for x2 in range(lat_size[0]): + coord[2] = x1 + coord[0] = x2 + + node = lat.get(coord) + if node.data[CUBA.MATERIAL] == solid: + continue + + p = node.data[CUBA.PRESSURE] + dn = node.data[CUBA.DENSITY] + v1 = node.data[CUBA.VELOCITY][2] + v2 = node.data[CUBA.VELOCITY][0] + + sf = '{0:d} {1:d} {2:e} {3:e} {4:e}\n' + f.write(sf.format(x1, x2, dn, v1, v2)) + +# --------------------------------------------------------------------------- diff --git a/examples/poiseuille_invel.py b/examples/poiseuille_invel.py new file mode 100644 index 0000000..5963c50 --- /dev/null +++ b/examples/poiseuille_invel.py @@ -0,0 +1,234 @@ +"""Poiseuille fluid flow simulation (flow between two parallel plates). + +Details +------- +- 3D simulation (with one "dummy" direction) +- flow driven by pressure gradient +- periodic bc in the flow direction (as well as in the "dummy" direction) +- wall bc (immobile & adiabatic) for representing the plates +""" +import time +import numpy as np +from simphony.cuds.meta import api +from simphony.core.cuba import CUBA +from simphony.api import CUDS, Simulation +from simphony.engine import EngineInterface +from simphony.cuds.lattice import make_cubic_lattice +from jyulb.flow_field import poise_plate_steady_state_umax + +# CUDS +cuds = CUDS(name='poiseuille') + +# Physics model +cfd = api.Cfd(name='fluid flow') + +# Submodels (in fact, these are the default submodels in CFD) +cfd.thermal_model = api.IsothermalModel(name='isothermal') +cfd.turbulence_model = api.LaminarFlowModel(name='laminar') +cfd.rheology_model = api.NewtonianFluidModel(name='newtonian') +cfd.multiphase_model = api.SinglePhaseModel(name='singlephase') +cfd.compressibility_model = api.IncompressibleFluidModel(name='incompressible') +cuds.add([cfd]) + +# Materials +fluid = api.Material(name='fluid') +solid = api.Material(name='wall') +fluid.data[CUBA.DENSITY] = 1.0 +fluid.data[CUBA.KINEMATIC_VISCOSITY] = 1.0/6.0 +cuds.add([fluid, solid]) + +# Dataset (lattice) +L = 10 +lat_size = (10, 1, L) # number of lattice nodes +lat_h = 1.0 # lattice spacing +lat = make_cubic_lattice('channel', lat_h, lat_size) + +# Create a channel geometry +for node in lat.iter(): + ijk = node.index + # The two plates are separated in the x-direction + if ijk[0] == 0 or ijk[0] == (lat_size[0]-1): + node.data[CUBA.MATERIAL] = solid + else: + node.data[CUBA.MATERIAL] = fluid + + lat.update([node]) + +cuds.add([lat]) + +# Boundary conditions +# for the exterior boundary faces of the simulation domain +# (the directions refer to unit outward normals) +face1y = api.Boundary(name='face1y', condition=[api.Periodic()]) +face1y.data[CUBA.DIRECTION] = (0, -1, 0) + +facey1 = api.Boundary(name='facey1', condition=[api.Periodic()]) +facey1.data[CUBA.DIRECTION] = (0, 1, 0) + +in_bc = api.Dirichlet(name='in_vel', material=fluid) +in_bc.data[CUBA.VARIABLE] = CUBA.VELOCITY + +face1z = api.Boundary(name='face1z', condition=[in_bc]) +face1z.data[CUBA.DIRECTION] = (0, 0, -1) + +out_bc = api.Dirichlet(name='out_p', material=fluid) +out_bc.data[CUBA.VARIABLE] = CUBA.PRESSURE + +facez1 = api.Boundary(name='facez1', condition=[out_bc]) +facez1.data[CUBA.DIRECTION] = (0, 0, 1) + +cuds.add([face1y, facey1, face1z, facez1]) + +# Solver parameters (time integration, collision operator) +ti = api.IntegrationTime(name='simulation time', + current=0.0, final=1000.0, size=1.0) + +ti.data[CUBA.COLLISION_OPERATOR] = 'TRT' +cuds.add([ti]) + +# Simulate +print '-'*77 +sim = Simulation(cuds, 'JYU-LB', engine_interface=EngineInterface.Internal) +lat = cuds.get_by_name(lat.name) + +# --------------------------------------------------------------------------- +# Analytical solution for the Poiseuille flow +# --------------------------------------------------------------------------- +max_vel = 1e-3 +H = lat_size[0]-2 # plates separated in the x-direction + +centrel_dist = np.zeros((H), dtype=np.float64) +ana_vel = np.zeros((H), dtype=np.float64) + +poise_plate_steady_state_umax(H, max_vel, centrel_dist, ana_vel) + +# Set the initial velocity field +for node in lat.iter(): + ijk = node.index + + if node.data[CUBA.MATERIAL] == fluid and ijk[2] == 0: + node.data[CUBA.VELOCITY] = (0, 0, ana_vel[ijk[0]-1]) + else: + node.data[CUBA.VELOCITY] = (0, 0, 0) + + lat.update([node]) + +# Total fluid mass +tot_fmass = 0.0 +for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_fmass += node.data[CUBA.DENSITY] + +print '='*77 +print 'Poiseuille flow simulation' +print '='*77 +print "Total fluid mass before simulation: {0:.4e}".format(tot_fmass) +print '-'*77 + +start = time.time() + +for e in range(10): + sim.run() + + # Total velocity + tot_vx = 0.0 + tot_vy = 0.0 + tot_vz = 0.0 + for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_vx += node.data[CUBA.VELOCITY][0] + tot_vy += node.data[CUBA.VELOCITY][1] + tot_vz += node.data[CUBA.VELOCITY][2] + + sf = 'Time = {0:f}, tot.vel = ({1:11.4e}, {2:11.4e}, {3:11.4e})' + print sf.format(ti.current, tot_vx, tot_vy, tot_vz) + +end = time.time() + +# Total fluid mass +tot_fmass = 0.0 +for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_fmass += node.data[CUBA.DENSITY] + +print '-'*77 +print "Total fluid mass after simulation: {0:.4e}".format(tot_fmass) +print '-'*77 +print "Time spend in run (s): ", end-start +print '-'*77 + +# --------------------------------------------------------------------------- +# Post-processing: error and output +# --------------------------------------------------------------------------- +# GLE: flow speed in y-data format +# --------------------------------------------------------------------------- +sim_loc_vel = np.zeros(3, dtype=np.float64) +coord = np.zeros(3, dtype=np.int32) +coord[0] = int(lat_size[0]/2) +coord[1] = int(lat_size[1]/2) +coord[2] = int(lat_size[2]/2) + +ana2 = 0.0 +sim_ana_diff2 = 0.0 +sim_loc_vel = np.zeros((3), dtype=np.float64) + +with open('flow_prof.txt', 'w') as f: + for h in range(H): + coord[0] = h+1 + node = lat.get(coord) + dist = centrel_dist[h] + + if node.data[CUBA.MATERIAL] == solid: + continue + + sim_loc_p = node.data[CUBA.PRESSURE] + sim_loc_den = node.data[CUBA.DENSITY] + sim_loc_vel[0] = node.data[CUBA.VELOCITY][0] + sim_loc_vel[1] = node.data[CUBA.VELOCITY][1] + sim_loc_vel[2] = node.data[CUBA.VELOCITY][2] + + sim_loc_speed = np.sqrt(sim_loc_vel[0]*sim_loc_vel[0] + + sim_loc_vel[1]*sim_loc_vel[1] + + sim_loc_vel[2]*sim_loc_vel[2]) + + ana_loc_vel = ana_vel[h] + + sim_ana_diff = sim_loc_vel[2] - ana_loc_vel + sim_ana_diff2 += sim_ana_diff*sim_ana_diff + ana2 += ana_loc_vel*ana_loc_vel + + f.write('{0:f} {1:e} {2:e} {3:e} {4:e} {5:e} {6:e}\n'.format(dist, + sim_loc_den, sim_loc_vel[0], sim_loc_vel[1], sim_loc_vel[2], + sim_loc_speed, ana_loc_vel)) + +rel_l2_err_norm_vel = np.sqrt(sim_ana_diff2/ana2) + +print 'Relative L2-error norm: {0:.4e}'.format(rel_l2_err_norm_vel) +print '-'*77 + +# --------------------------------------------------------------------------- +# Matlab: 2D flow field +# --------------------------------------------------------------------------- +coord[0] = int(lat_size[0]/2) +coord[1] = int(lat_size[1]/2) +coord[2] = int(lat_size[2]/2) + +with open('flow_2D.field', 'w') as f: + for x1 in range(lat_size[2]): + for x2 in range(lat_size[0]): + coord[2] = x1 + coord[0] = x2 + + node = lat.get(coord) + if node.data[CUBA.MATERIAL] == solid: + continue + + p = node.data[CUBA.PRESSURE] + dn = node.data[CUBA.DENSITY] + v1 = node.data[CUBA.VELOCITY][2] + v2 = node.data[CUBA.VELOCITY][0] + + sf = '{0:d} {1:d} {2:e} {3:e} {4:e}\n' + f.write(sf.format(x1, x2, dn, v1, v2)) + +# --------------------------------------------------------------------------- diff --git a/examples/poiseuille_pdiff.py b/examples/poiseuille_pdiff.py new file mode 100644 index 0000000..b3e4719 --- /dev/null +++ b/examples/poiseuille_pdiff.py @@ -0,0 +1,237 @@ +"""Poiseuille fluid flow simulation (flow between two parallel plates). + +Details +------- +- 3D simulation (with one "dummy" direction) +- flow driven by pressure gradient +- periodic bc in the flow direction (as well as in the "dummy" direction) +- wall bc (immobile & adiabatic) for representing the plates +""" +import time +import numpy as np +from simphony.cuds.meta import api +from simphony.core.cuba import CUBA +from simphony.api import CUDS, Simulation +from simphony.engine import EngineInterface +from simphony.cuds.lattice import make_cubic_lattice +from jyulb.flow_field import poise_plate_steady_state_pgrad + +# CUDS +cuds = CUDS(name='poiseuille') + +# Physics model +cfd = api.Cfd(name='fluid flow') + +# Submodels (in fact, these are the default submodels in CFD) +cfd.thermal_model = api.IsothermalModel(name='isothermal') +cfd.turbulence_model = api.LaminarFlowModel(name='laminar') +cfd.rheology_model = api.NewtonianFluidModel(name='newtonian') +cfd.multiphase_model = api.SinglePhaseModel(name='singlephase') +cfd.compressibility_model = api.IncompressibleFluidModel(name='incompressible') +cuds.add([cfd]) + +# Materials +fluid = api.Material(name='fluid') +solid = api.Material(name='wall') +fluid.data[CUBA.DENSITY] = 1.0 +fluid.data[CUBA.KINEMATIC_VISCOSITY] = 1.0/6.0 +cuds.add([fluid, solid]) + +# Dataset (lattice) +L = 10 +lat_size = (10, 1, L) # number of lattice nodes +lat_h = 1.0 # lattice spacing +lat = make_cubic_lattice('channel', lat_h, lat_size) + +# Create a channel geometry +for node in lat.iter(): + ijk = node.index + # The two plates are separated in the x-direction + if ijk[0] == 0 or ijk[0] == (lat_size[0]-1): + node.data[CUBA.MATERIAL] = solid + else: + node.data[CUBA.MATERIAL] = fluid + + lat.update([node]) + +cuds.add([lat]) + +# Boundary conditions +# for the exterior boundary faces of the simulation domain +# (the directions refer to unit outward normals) +face1y = api.Boundary(name='face1y', condition=[api.Periodic()]) +face1y.data[CUBA.DIRECTION] = (0, -1, 0) + +facey1 = api.Boundary(name='facey1', condition=[api.Periodic()]) +facey1.data[CUBA.DIRECTION] = (0, 1, 0) + +in_bc = api.Dirichlet(name='in_p', material=fluid) +in_bc.data[CUBA.VARIABLE] = CUBA.PRESSURE + +face1z = api.Boundary(name='face1z', condition=[in_bc]) +face1z.data[CUBA.DIRECTION] = (0, 0, -1) + +out_bc = api.Dirichlet(name='out_p', material=fluid) +out_bc.data[CUBA.VARIABLE] = CUBA.PRESSURE + +facez1 = api.Boundary(name='facez1', condition=[out_bc]) +facez1.data[CUBA.DIRECTION] = (0, 0, 1) + +cuds.add([face1y, facey1, face1z, facez1]) + +# Solver parameters (time integration, collision operator) +ti = api.IntegrationTime(name='simulation time', + current=0.0, final=1000.0, size=1.0) + +ti.data[CUBA.COLLISION_OPERATOR] = 'TRT' +cuds.add([ti]) + +# Simulate +print '-'*77 +sim = Simulation(cuds, 'JYU-LB', engine_interface=EngineInterface.Internal) +lat = cuds.get_by_name(lat.name) + +# Pressure gradient +kinem_pgrad = -1e-5 # gradient of kinematic pressure +pgrad = fluid.data[CUBA.DENSITY]*kinem_pgrad # gradient of pressure +pdiff = (L-1.0)*pgrad # pressure difference between outlet and inlet +pin_off = -0.5*pdiff # inlet pressure (offset to the reference pressure) +pout_off = 0.5*pdiff # outlet pressure (offset to the reference pressure) + +# Set the initial pressure and velocity field +for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + ijk = node.index + p = node.data[CUBA.PRESSURE] + node.data[CUBA.PRESSURE] = p + pin_off + ijk[2]*pgrad + node.data[CUBA.VELOCITY] = (0, 0, 0) + lat.update([node]) + +# Total fluid mass +tot_fmass = 0.0 +for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_fmass += node.data[CUBA.DENSITY] + +print '='*77 +print 'Poiseuille flow simulation' +print '='*77 +print "Total fluid mass before simulation: {0:.4e}".format(tot_fmass) +print '-'*77 + +start = time.time() + +for e in range(10): + sim.run() + + # Total velocity + tot_vx = 0.0 + tot_vy = 0.0 + tot_vz = 0.0 + for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_vx += node.data[CUBA.VELOCITY][0] + tot_vy += node.data[CUBA.VELOCITY][1] + tot_vz += node.data[CUBA.VELOCITY][2] + + sf = 'Time = {0:f}, tot.vel = ({1:11.4e}, {2:11.4e}, {3:11.4e})' + print sf.format(ti.current, tot_vx, tot_vy, tot_vz) + +end = time.time() + +# Total fluid mass +tot_fmass = 0.0 +for node in lat.iter(): + if node.data[CUBA.MATERIAL] == fluid: + tot_fmass += node.data[CUBA.DENSITY] + +print '-'*77 +print "Total fluid mass after simulation: {0:.4e}".format(tot_fmass) +print '-'*77 +print "Time spend in run (s): ", end-start +print '-'*77 + +# --------------------------------------------------------------------------- +# Post-processing: +# Analytical solution for the Poiseuille flow, error, and output +# --------------------------------------------------------------------------- +# GLE: flow speed in y-data format +# --------------------------------------------------------------------------- +H = lat_size[0]-2 # plates separated in the x-direction +dvisc = fluid.data[CUBA.DENSITY]*fluid.data[CUBA.KINEMATIC_VISCOSITY] + +centrel_dist = np.zeros((H), dtype=np.float64) +ana_vel = np.zeros((H), dtype=np.float64) + +poise_plate_steady_state_pgrad(H, dvisc, pgrad, centrel_dist, ana_vel) + +sim_loc_vel = np.zeros(3, dtype=np.float64) +coord = np.zeros(3, dtype=np.int32) +coord[0] = int(lat_size[0]/2) +coord[1] = int(lat_size[1]/2) +coord[2] = int(lat_size[2]/2) + +ana2 = 0.0 +sim_ana_diff2 = 0.0 +sim_loc_vel = np.zeros((3), dtype=np.float64) + +with open('flow_prof.txt', 'w') as f: + for h in range(H): + coord[0] = h+1 + node = lat.get(coord) + dist = centrel_dist[h] + + if node.data[CUBA.MATERIAL] == solid: + continue + + sim_loc_p = node.data[CUBA.PRESSURE] + sim_loc_den = node.data[CUBA.DENSITY] + sim_loc_vel[0] = node.data[CUBA.VELOCITY][0] + sim_loc_vel[1] = node.data[CUBA.VELOCITY][1] + sim_loc_vel[2] = node.data[CUBA.VELOCITY][2] + + sim_loc_speed = np.sqrt(sim_loc_vel[0]*sim_loc_vel[0] + + sim_loc_vel[1]*sim_loc_vel[1] + + sim_loc_vel[2]*sim_loc_vel[2]) + + ana_loc_vel = ana_vel[h] + + sim_ana_diff = sim_loc_vel[2] - ana_loc_vel + sim_ana_diff2 += sim_ana_diff*sim_ana_diff + ana2 += ana_loc_vel*ana_loc_vel + + f.write('{0:f} {1:e} {2:e} {3:e} {4:e} {5:e} {6:e}\n'.format(dist, + sim_loc_den, sim_loc_vel[0], sim_loc_vel[1], sim_loc_vel[2], + sim_loc_speed, ana_loc_vel)) + +rel_l2_err_norm_vel = np.sqrt(sim_ana_diff2/ana2) + +print 'Relative L2-error norm: {0:.4e}'.format(rel_l2_err_norm_vel) +print '-'*77 + +# --------------------------------------------------------------------------- +# Matlab: 2D flow field +# --------------------------------------------------------------------------- +coord[0] = int(lat_size[0]/2) +coord[1] = int(lat_size[1]/2) +coord[2] = int(lat_size[2]/2) + +with open('flow_2D.field', 'w') as f: + for x1 in range(lat_size[2]): + for x2 in range(lat_size[0]): + coord[2] = x1 + coord[0] = x2 + + node = lat.get(coord) + if node.data[CUBA.MATERIAL] == solid: + continue + + p = node.data[CUBA.PRESSURE] + dn = node.data[CUBA.DENSITY] + v1 = node.data[CUBA.VELOCITY][2] + v2 = node.data[CUBA.VELOCITY][0] + + sf = '{0:d} {1:d} {2:e} {3:e} {4:e}\n' + f.write(sf.format(x1, x2, dn, v1, v2)) + +# --------------------------------------------------------------------------- diff --git a/jyulb/__init__.py b/jyulb_cpp/__init__.py similarity index 100% rename from jyulb/__init__.py rename to jyulb_cpp/__init__.py diff --git a/jyulb/cuba_extension.py b/jyulb_cpp/cuba_extension.py similarity index 100% rename from jyulb/cuba_extension.py rename to jyulb_cpp/cuba_extension.py diff --git a/jyulb/fileio/__init__.py b/jyulb_cpp/fileio/__init__.py similarity index 100% rename from jyulb/fileio/__init__.py rename to jyulb_cpp/fileio/__init__.py diff --git a/jyulb/fileio/common/__init__.py b/jyulb_cpp/fileio/common/__init__.py similarity index 100% rename from jyulb/fileio/common/__init__.py rename to jyulb_cpp/fileio/common/__init__.py diff --git a/jyulb/fileio/common/proxy_lattice.py b/jyulb_cpp/fileio/common/proxy_lattice.py similarity index 100% rename from jyulb/fileio/common/proxy_lattice.py rename to jyulb_cpp/fileio/common/proxy_lattice.py diff --git a/jyulb/fileio/common/tests/__init__.py b/jyulb_cpp/fileio/common/tests/__init__.py similarity index 100% rename from jyulb/fileio/common/tests/__init__.py rename to jyulb_cpp/fileio/common/tests/__init__.py diff --git a/jyulb/fileio/common/tests/test_proxy_lattice.py b/jyulb_cpp/fileio/common/tests/test_proxy_lattice.py similarity index 100% rename from jyulb/fileio/common/tests/test_proxy_lattice.py rename to jyulb_cpp/fileio/common/tests/test_proxy_lattice.py diff --git a/jyulb/fileio/isothermal/__init__.py b/jyulb_cpp/fileio/isothermal/__init__.py similarity index 100% rename from jyulb/fileio/isothermal/__init__.py rename to jyulb_cpp/fileio/isothermal/__init__.py diff --git a/jyulb/fileio/isothermal/jyulb_engine.py b/jyulb_cpp/fileio/isothermal/jyulb_engine.py similarity index 100% rename from jyulb/fileio/isothermal/jyulb_engine.py rename to jyulb_cpp/fileio/isothermal/jyulb_engine.py diff --git a/jyulb/fileio/isothermal/tests/__init__.py b/jyulb_cpp/fileio/isothermal/tests/__init__.py similarity index 100% rename from jyulb/fileio/isothermal/tests/__init__.py rename to jyulb_cpp/fileio/isothermal/tests/__init__.py diff --git a/jyulb/fileio/isothermal/tests/test_jyulb_engine.py b/jyulb_cpp/fileio/isothermal/tests/test_jyulb_engine.py similarity index 100% rename from jyulb/fileio/isothermal/tests/test_jyulb_engine.py rename to jyulb_cpp/fileio/isothermal/tests/test_jyulb_engine.py diff --git a/jyulb/fileio/isothermal/tests/test_plugin_integration.py b/jyulb_cpp/fileio/isothermal/tests/test_plugin_integration.py similarity index 100% rename from jyulb/fileio/isothermal/tests/test_plugin_integration.py rename to jyulb_cpp/fileio/isothermal/tests/test_plugin_integration.py diff --git a/jyulb/internal/__init__.py b/jyulb_cpp/internal/__init__.py similarity index 100% rename from jyulb/internal/__init__.py rename to jyulb_cpp/internal/__init__.py diff --git a/jyulb/internal/common/__init__.py b/jyulb_cpp/internal/common/__init__.py similarity index 100% rename from jyulb/internal/common/__init__.py rename to jyulb_cpp/internal/common/__init__.py diff --git a/jyulb/internal/common/domain.pxd b/jyulb_cpp/internal/common/domain.pxd similarity index 100% rename from jyulb/internal/common/domain.pxd rename to jyulb_cpp/internal/common/domain.pxd diff --git a/jyulb/internal/common/domain.pyx b/jyulb_cpp/internal/common/domain.pyx similarity index 100% rename from jyulb/internal/common/domain.pyx rename to jyulb_cpp/internal/common/domain.pyx diff --git a/jyulb_cpp/internal/common/domain_setup.py b/jyulb_cpp/internal/common/domain_setup.py new file mode 100644 index 0000000..13e72c9 --- /dev/null +++ b/jyulb_cpp/internal/common/domain_setup.py @@ -0,0 +1,14 @@ +from distutils.core import setup, Extension +from Cython.Build import cythonize +import numpy + +setup(ext_modules = cythonize(Extension( + "domain", # the extesion name + sources=["domain.pyx", # the Cython and C++ source files + "../../../../JYU-LB/include/common/node.cpp"], + include_dirs=['../../../../JYU-LB/include/common/', # include paths + numpy.get_include()], + extra_compile_args=['-fopenmp','-O3'], # compiler flags + extra_link_args=['-fopenmp','-O3'], # linker flags + language="c++", # generate and compile C++ code +))) diff --git a/jyulb/internal/common/proxy_lattice.py b/jyulb_cpp/internal/common/proxy_lattice.py similarity index 100% rename from jyulb/internal/common/proxy_lattice.py rename to jyulb_cpp/internal/common/proxy_lattice.py diff --git a/jyulb_cpp/internal/common/test_domain.py b/jyulb_cpp/internal/common/test_domain.py new file mode 100644 index 0000000..348c2cc --- /dev/null +++ b/jyulb_cpp/internal/common/test_domain.py @@ -0,0 +1,93 @@ +import domain +import time +import numpy as np + +print 'Solid material id = {}'.format(domain.SOLID_ENUM) +print 'Fluid material id = {}'.format(domain.FLUID_ENUM) + +lat1 = domain.PyLattice(np.array((4, 2, 3), dtype=np.uint32), + np.array((0.0, 1.0, 2.0), dtype=np.float64)) + +size1 = np.zeros(3, dtype=np.uint32) +origin1 = np.zeros(3, dtype=np.float64) + +lat1.get_size(size1) +lat1.get_origin(origin1) +ncount1 = lat1.get_node_count() + +print type(size1), size1, ncount1 +print type(origin1), origin1 + +lat2 = domain.PyLattice.fromlattice(lat1) + +size2 = np.zeros(3, dtype=np.uint32) +origin2 = np.zeros(3, dtype=np.float64) + +lat2.get_size(size2) +lat2.get_origin(origin2) +ncount2 = lat2.get_node_count() + +print type(size2), size2, ncount2 +print type(origin2), origin2 + +for index in np.ndindex(tuple(size1)): + print index, lat1.get_n(np.array(index, dtype=np.uint32)) + +for n in np.arange(ncount1): + ijk = np.zeros(3, dtype=np.uint32) + lat1.get_ijk(n, ijk) + print n, ijk + +geom1 = domain.PyGeometry(lat1) + +for index in np.ndindex(tuple(size1)): + if index[0] > 0 and index[0] < size1[0]-1: + geom1.set_material_ijk(np.array(index, dtype=np.uint32), domain.FLUID_ENUM) + +fluid_ncount = 0 +solid_ncount = 0 +for index in np.ndindex(tuple(size1)): + material = geom1.get_material_ijk(np.array(index, dtype=np.uint32)) + print 'Node {} material id = {}'.format(index, material) + + if material == domain.SOLID_ENUM: + solid_ncount = solid_ncount + 1 + else: + fluid_ncount = fluid_ncount + 1 + +print 'Solid node count = {}'.format(solid_ncount) +print 'Fluid node count = {}'.format(fluid_ncount) + +ivel = np.array((0.0, 1.0, 2.0), dtype=np.float64) +ifrc = np.array((0.0, 0.0, 0.0), dtype=np.float64) +field_data = domain.PyIsothermalData(lat1, 1.0, ivel, ifrc) + +#data_count = field_data.get_node_set().get_node_count() +data_count = field_data.get_data_count() +print 'Isothermal data count = {}'.format(data_count) + +vel = np.zeros(3, dtype=np.float64) +frc = np.ones(3, dtype=np.float64) +ijk = np.zeros(3, dtype=np.uint32) + +for n in np.arange(data_count): +# field_data.get_node_set().get_ijk(n, ijk) + field_data.get_ijk(n, ijk) + + vel[0] = n*1.0 + vel[1] = 0.0 + vel[2] = 0.0 + + frc[0] = 0.0 + frc[1] = 0.0 + frc[2] = n*1.0 + + field_data.set_vel_ijk(ijk, vel) + field_data.set_frc_ijk(ijk, frc) + + den = field_data.get_den_n(n) + field_data.get_vel_n(n, vel) + field_data.get_frc_n(n, frc) + + print n, ijk, den, vel, frc + diff --git a/jyulb/internal/common/tests/__init__.py b/jyulb_cpp/internal/common/tests/__init__.py similarity index 100% rename from jyulb/internal/common/tests/__init__.py rename to jyulb_cpp/internal/common/tests/__init__.py diff --git a/jyulb_cpp/internal/common/tests/test_lattice.py b/jyulb_cpp/internal/common/tests/test_lattice.py new file mode 100644 index 0000000..56b17b0 --- /dev/null +++ b/jyulb_cpp/internal/common/tests/test_lattice.py @@ -0,0 +1,153 @@ +""" + Testing module for jyu-lb proxy-lattice data class. +""" +import unittest +import numpy as np +import numpy.testing as np_test + +from simphony.core.cuba import CUBA +from simphony.cuds.lattice import LatticeNode +from simphony.cuds.abstractlattice import ABCLattice +from simphony.testing.abc_check_lattice import ABCCheckLattice +from simphony.core.data_container import DataContainer +from jyulb.internal.common.proxy_lattice import ProxyLattice +from jyulb.internal.common.domain import PyLattice, PyGeometry +from jyulb.internal.common.domain import PyIsothermalData + + +class TestProxyLattice(ABCCheckLattice, unittest.TestCase): + +# def setUp(self): +# self.hl_data = DataContainer() +# self.hl_data[CUBA.STATUS] = 0 + + def container_factory(self, name, type_, base_vect, size, origin): + size_n = len(size) + orig_n = len(origin) +# np_size = np.array(size, dtype=np.uint32) +# np_orig = np.array(origin, dtype=np.float64) + np_size = np.ones(3, dtype=np.uint32) + np_orig = np.ones(3, dtype=np.float64) + + np_size[0:len(size)] = size[:] + np_orig[0:len(origin)] = origin[:] + +# if size_n > 0: +# np_size[0] = size[0] +# if size_n > 1: +# np_size[1] = size[1] +# if size_n > 2: +# np_size[2] = size[2] + +# if orig_n > 0: +# np_orig[0] = origin[0] +# if orig_n > 1: +# np_orig[1] = origin[1] +# if orig_n > 2: +# np_orig[2] = origin[2] + + self.lat = PyLattice(np_size, np_orig) + self.geom = PyGeometry(self.lat) + + for index in np.ndindex(tuple(np_size)): + if index[0] > 0 and index[0] < np_size[0]-1: + ijk = np.array(index, dtype=np.uint32) + self.geom.set_material_ijk(ijk, ProxyLattice.FLUID_ENUM) + + ivel = np.array((0.0, 0.0, 0.0), dtype=np.float64) + ifrc = np.array((0.0, 0.0, 0.0), dtype=np.float64) + self.fdata = PyIsothermalData(self.lat, 0.0, ivel, ifrc) + + return ProxyLattice(name, type_, base_vect, self.geom, self.fdata) + + def supported_cuba(self): + return set((CUBA.MATERIAL_ID,CUBA.DENSITY,CUBA.VELOCITY,CUBA.FORCE)) + + def test_construct_lattice(self): + """Construction of a lattice.""" + proxy = self.container + +# proxy = self.container_factory('test_lat', 'Cubic', +# (0.1, 0.1, 0.1), (10, 20, 30), +# (0, 0, 0)) +# proxy.data = self.hl_data + + self.assertIsInstance(proxy, ABCLattice, "Error: not a ABCLattice!") + + self.assertEqual(proxy._geom, self.geom) + self.assertEqual(proxy._fdata, self.fdata) +# self.assertEqual(proxy.data, self.hl_data) + + try: + proxy.get_node((-1, 0, 0)) + except IndexError: + pass + else: + raise AssertionError('Negative indices incorrectly accepted.') + + test_node = LatticeNode((0, -1, 0)) + try: + proxy.update_node(test_node) + except IndexError: + pass + else: + raise AssertionError('Negative indices incorrectly accepted.') + + def test_set_get_iter_lattice_nodes(self): + """Creation of lattices using the factory functions.""" + proxy = ProxyLattice(self.name, self.type, + self.bvec, self.geom, self.fdata) + + iden = 1.0 + ivel = np.array((0.0, 0.0, 1.0), dtype=np.float64) + ifrc = np.array((0.0, 1.0, 0.0), dtype=np.float64) + + for node in proxy.iter_nodes(): + node.data[CUBA.DENSITY] = iden + node.data[CUBA.VELOCITY] = ivel + node.data[CUBA.FORCE] = ifrc + proxy.update_node(node) + + check_iden = 0.0 + check_ivelx = 0.0 + check_ively = 0.0 + check_ivelz = 0.0 + check_ifrcx = 0.0 + check_ifrcy = 0.0 + check_ifrcz = 0.0 + check_sum = (self.nx-2)*self.ny*self.nz + + for node in proxy.iter_nodes(): + if node.data[CUBA.MATERIAL_ID] == ProxyLattice.FLUID_ENUM: + check_iden = check_iden + node.data[CUBA.DENSITY] + check_ivelx = check_ivelx + node.data[CUBA.VELOCITY][0] + check_ively = check_ively + node.data[CUBA.VELOCITY][1] + check_ivelz = check_ivelz + node.data[CUBA.VELOCITY][2] + check_ifrcx = check_ifrcx + node.data[CUBA.FORCE][0] + check_ifrcy = check_ifrcy + node.data[CUBA.FORCE][1] + check_ifrcz = check_ifrcz + node.data[CUBA.FORCE][2] + + self.assertEqual(check_sum, check_iden) + self.assertEqual(0.0, check_ivelx) + self.assertEqual(0.0, check_ively) + self.assertEqual(check_sum, check_ivelz) + self.assertEqual(0.0, check_ifrcx) + self.assertEqual(check_sum, check_ifrcy) + self.assertEqual(0.0, check_ifrcz) + + for node in proxy.iter_nodes(np.ndindex(1, self.ny, self.nz)): + try: + node.data[CUBA.DENSITY] + node.data[CUBA.VELOCITY] + node.data[CUBA.FORCE] + except: + pass + else: + raise AssertionError('Field data stored for solid nodes.') + + np_test.assert_array_equal((1*self.h, 2*self.h, 3*self.h), + proxy.get_coordinate((1, 2, 3))) + + +if __name__ == '__main__': + unittest.main() diff --git a/jyulb/internal/common/tests/test_proxy_lattice.py b/jyulb_cpp/internal/common/tests/test_proxy_lattice.py similarity index 100% rename from jyulb/internal/common/tests/test_proxy_lattice.py rename to jyulb_cpp/internal/common/tests/test_proxy_lattice.py diff --git a/jyulb/internal/isothermal/__init__.py b/jyulb_cpp/internal/isothermal/__init__.py similarity index 98% rename from jyulb/internal/isothermal/__init__.py rename to jyulb_cpp/internal/isothermal/__init__.py index 7b15fcf..46e7496 100644 --- a/jyulb/internal/isothermal/__init__.py +++ b/jyulb_cpp/internal/isothermal/__init__.py @@ -1,2 +1,3 @@ from jyulb.internal.isothermal import jyulb_engine + __all__ = ['jyulb_engine'] diff --git a/jyulb/internal/isothermal/jyulb_engine.py b/jyulb_cpp/internal/isothermal/jyulb_engine.py similarity index 100% rename from jyulb/internal/isothermal/jyulb_engine.py rename to jyulb_cpp/internal/isothermal/jyulb_engine.py diff --git a/jyulb/internal/isothermal/solver.pxd b/jyulb_cpp/internal/isothermal/solver.pxd similarity index 100% rename from jyulb/internal/isothermal/solver.pxd rename to jyulb_cpp/internal/isothermal/solver.pxd diff --git a/jyulb/internal/isothermal/solver.pyx b/jyulb_cpp/internal/isothermal/solver.pyx similarity index 100% rename from jyulb/internal/isothermal/solver.pyx rename to jyulb_cpp/internal/isothermal/solver.pyx diff --git a/jyulb_cpp/internal/isothermal/solver_setup.py b/jyulb_cpp/internal/isothermal/solver_setup.py new file mode 100644 index 0000000..e80ed8f --- /dev/null +++ b/jyulb_cpp/internal/isothermal/solver_setup.py @@ -0,0 +1,22 @@ +from distutils.core import setup, Extension +from Cython.Build import cythonize +import numpy + +setup(ext_modules = cythonize(Extension( + "solver", # the extesion name + sources=["solver.pyx", # the Cython and C++ source files + "../../../../JYU-LB/include/common/node.cpp", + "../../../../JYU-LB/include/common/filter.cpp", + "../../../../JYU-LB/include/collision/collision.cpp", + "../../../../JYU-LB/include/kernel/kernel.cpp", + "../../../../JYU-LB/include/solver/solver.cpp"], + include_dirs=['../../../../JYU-LB/include/common/', # include paths + '../../../../JYU-LB/include/dvs/', + '../../../../JYU-LB/include/collision/', + '../../../../JYU-LB/include/kernel/', + '../../../../JYU-LB/include/solver/', + numpy.get_include()], + extra_compile_args=['-fopenmp','-O3'], # compiler flags + extra_link_args=['-fopenmp','-O3'], # linker flags + language="c++", # generate and compile C++ code +))) diff --git a/jyulb/internal/isothermal/tests/__init__.py b/jyulb_cpp/internal/isothermal/tests/__init__.py similarity index 100% rename from jyulb/internal/isothermal/tests/__init__.py rename to jyulb_cpp/internal/isothermal/tests/__init__.py diff --git a/jyulb/internal/isothermal/tests/test_jyulb_engine.py b/jyulb_cpp/internal/isothermal/tests/test_jyulb_engine.py similarity index 100% rename from jyulb/internal/isothermal/tests/test_jyulb_engine.py rename to jyulb_cpp/internal/isothermal/tests/test_jyulb_engine.py diff --git a/jyulb/internal/isothermal/tests/test_plugin_integration.py b/jyulb_cpp/internal/isothermal/tests/test_plugin_integration.py similarity index 100% rename from jyulb/internal/isothermal/tests/test_plugin_integration.py rename to jyulb_cpp/internal/isothermal/tests/test_plugin_integration.py diff --git a/jyulb/testing/__init__.py b/jyulb_cpp/testing/__init__.py similarity index 100% rename from jyulb/testing/__init__.py rename to jyulb_cpp/testing/__init__.py diff --git a/jyulb/testing/jyulb_check_engine.py b/jyulb_cpp/testing/jyulb_check_engine.py similarity index 100% rename from jyulb/testing/jyulb_check_engine.py rename to jyulb_cpp/testing/jyulb_check_engine.py diff --git a/jyulb/testing/jyulb_check_proxy_lattice.py b/jyulb_cpp/testing/jyulb_check_proxy_lattice.py similarity index 100% rename from jyulb/testing/jyulb_check_proxy_lattice.py rename to jyulb_cpp/testing/jyulb_check_proxy_lattice.py diff --git a/jyulb_cpp/version.py b/jyulb_cpp/version.py new file mode 100644 index 0000000..01673f6 --- /dev/null +++ b/jyulb_cpp/version.py @@ -0,0 +1 @@ +version = '0.2.1' diff --git a/requirements.txt b/requirements.txt index 1ab95b6..eb6999a 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,2 +1,3 @@ numpy>=1.9.1 -cython +cython>=0.25.2 +jyulb>=0.2.0 \ No newline at end of file diff --git a/setup.py b/setup.py index 632d84f..0cef935 100644 --- a/setup.py +++ b/setup.py @@ -1,72 +1,29 @@ import os -from setuptools import setup, find_packages, Extension -from Cython.Build import cythonize -import numpy - -VERSION = '0.2.1' +from setuptools import setup, find_packages +VERSION = '0.3.0' def write_version_py(filename=None): if filename is None: filename = os.path.join( - os.path.dirname(__file__), 'jyulb', 'version.py') - ver = """\ -version = '%s' -""" + os.path.dirname(__file__), 'wrapper', 'version.py') + + ver = """version = '%s'\n""" fh = open(filename, 'wb') try: fh.write(ver % VERSION) finally: fh.close() - write_version_py() -extensions = [ - Extension( - name="jyulb.internal.isothermal.solver", - sources=["jyulb/internal/isothermal/solver.pyx", - "JYU-LB/include/common/node.cpp", - "JYU-LB/include/common/filter.cpp", - "JYU-LB/include/collision/collision.cpp", - "JYU-LB/include/kernel/kernel.cpp", - "JYU-LB/include/solver/solver.cpp"], - include_dirs=['JYU-LB/include/common/', - 'JYU-LB/include/dvs/', - 'JYU-LB/include/collision/', - 'JYU-LB/include/kernel/', - 'JYU-LB/include/solver/', - numpy.get_include()], - extra_compile_args=['-fopenmp', '-O3'], - extra_link_args=['-fopenmp', '-O3'], - language="c++", - ), - Extension( - name="jyulb.internal.common.domain", - sources=["jyulb/internal/common/domain.pyx", - "JYU-LB/include/common/node.cpp"], - include_dirs=['JYU-LB/include/common/', - numpy.get_include()], - extra_compile_args=['-fopenmp', '-O3'], - extra_link_args=['-fopenmp', '-O3'], - language="c++", - ) -] - setup( - name='jyulb_engine', + name='simphony_jyulb', version=VERSION, - author='SimPhoNy FP7 European Project', - description='Implementation of JYU-LB wrappers', + author='SimPhoNy, EU FP7 Project (Nr. 604005) www.simphony-project.eu', + description='JYU-LB wrapper for the SimPhoNy framework', packages=find_packages(), - install_requires=['simphony~=0.4'], - ext_modules=cythonize(extensions), - entry_points={ - 'simphony.engine': [ - 'jyulb_fileio_isothermal =' + - 'jyulb.fileio.isothermal.jyulb_engine', - 'jyulb_internal_isothermal =' + - 'jyulb.internal.isothermal.jyulb_engine' - ] - }, +# install_requires=['numpy>=1.9.1','simphony>=0.6','jyulb>=0.2.0'], + install_requires=['numpy>=1.9.1','cython>=0.25.2','simphony>=0.6'], + entry_points={'simphony.engine': ['jyulb = wrapper']} ) diff --git a/wrapper/__init__.py b/wrapper/__init__.py new file mode 100644 index 0000000..24a2b21 --- /dev/null +++ b/wrapper/__init__.py @@ -0,0 +1,65 @@ +from simphony.engine import ABCEngineExtension +from simphony.engine import EngineInterface +from simphony.engine.decorators import register + +from .jyulb_engine import Wrapper as InternalWrapper + +__all__ = ["InternalWrapper"] + + +@register +class JYULBExtension(ABCEngineExtension): + + """JYU-LB extension. + + This extension provides support for JYU-LB engine. + """ + + def get_supported_engines(self): + """Get metadata about supported engines. + + Returns + ------- + list: a list of EngineMetadata objects + """ + # TODO: Add proper features as soon as the metadata classes are ready. + # Flow type, relaxation model etc. + # jyulb_features =\ + # self.create_engine_metadata_feature(LAMINAR_FLOW, TRT) + + jyulb_features = None + + jyulb = self.create_engine_metadata('JYU-LB', + jyulb_features, + [EngineInterface.Internal]) + return [jyulb] + + def create_wrapper(self, cuds, engine_name, engine_interface): + """Create a wrapper to the requested engine. + + Parameters + ---------- + cuds: CUDS + CUDS computational model data + engine_name: str + name of the engine, must be supported by this extension + engine_interface: EngineInterface + the interface to interact with engine + + Returns + ------- + ABCEngineExtension: A wrapper configured with cuds and ready to run + """ + use_internal_interface = True + if engine_interface == EngineInterface.FileIO: + use_internal_interface = False + + if engine_name == 'JYU-LB': + if use_internal_interface: + return InternalWrapper(cuds=cuds) + else: + raise Exception('File-IO wrapper for the JYU-LB engine ' + 'is not supported.') + else: + raise Exception('Only JYU-LB engine is supported. ' + 'Unsupported engine: %s', engine_name) diff --git a/wrapper/jyulb_engine.py b/wrapper/jyulb_engine.py new file mode 100644 index 0000000..0d5c14e --- /dev/null +++ b/wrapper/jyulb_engine.py @@ -0,0 +1,688 @@ +"""JYU-LB wrapper for the SimPhoNy framework. + +The wrapper relies on the jyulb-extension module +(i.e. an internal interface wrapper). +""" +import numpy as np + +from simphony.cuds.meta import api +from simphony.core.cuba import CUBA +from simphony.cuds.abc_lattice import ABCLattice +from simphony.cuds.primitive_cell import BravaisLattice +from simphony.cuds.abc_modeling_engine import ABCModelingEngine + +from jyulb.isothermal import Solver +from jyulb.defs import FACE_1X, FACE_X1, FACE_1Y, FACE_Y1, FACE_1Z, FACE_Z1 +from jyulb.defs import PERIODIC, FIXED_DEN, FIXED_VEL +from jyulb.defs import BGK, TRT, REGULARIZATION +from jyulb.defs import FLUID, SOLID +from jyulb.defs import X, Y, Z + +from .jyulb_lattice import ProxyLattice + + +class Wrapper(ABCModelingEngine): + + """JYU-LB wrapper (internal interface) for the SimPhoNy framework. + + CUDS configuration + ================== + Physics model + ------------- + - Expecting CFD (exactly one) with submodels + incompressible, + single phase, + isothermal, + newtonian, + laminar + - Acceleration due to gravity is read from the Gravity submodel of CFD + (i.e. Gravity models not associated with the CFD are ignored); + acceleration due to gravity is updated before each run + - The electrostatic submodel of CFD is ignored + + Materials + --------- + - Assume 1 (fluid) or 2 (fluid and solid/wall) materials + - If 2 materials, one material must be named 'solid' or 'wall' + - The non-solid/non-wall material is interpreted as fluid with data + DENSITY or REFERENCE_DENSITY and + KINEMATIC_VISCOSITY or DYNAMIC_VISCOSITY or + VISCOSITY (interpreted as DYNAMIC_VISCOSITY) + + Dataset + ------- + - Cubic Bravais lattice (exactly one), + length of a primitive vector defines the lattice spacing + - Lattice node data acknowledged includes + MATERIAL (the default material is fluid if not specified), + PRESSURE (or alternatively DENSITY), + VELOCITY (or alternatively MOMENTUM), + ACCELERATION (or alternatively FORCE); + see also the documentation of ProxyLattice + + Solver parameters + ----------------- + - TimeIntegration (exactly one), + the attribute size defines the discrete time step; + attributes final and current are updated by the wrapper + after each run + - The collision operator can be specified by providing + a SolverParameter with data COLLISION_OPERATOR + (value 'bgk' or 'trt' or 'regularization'); + e.g. the TimeIntegration can be used for this purpose + + Boundaries + ---------- + - In order to configure inlet/outlet/etc. boundary conditions, + exterior/outer boundary faces of the simulation domain + must be specified + - Exterior boundary faces of the simulation domain are specified + by providing Boundaries with particular unit outward normals, + i.e. expecting Boundaries with data DIRECTION = (-1,0,0) or + (1,0,0) or (0,-1,0) or (0,1,0) or (0,0,-1) or (0,0,1) + - Two or more Boundaries with equal DIRECTION value (see above) + is not allowed + + Conditions + ---------- + - Each Boundary (see above) can be associated with conditions + (i.e. conditions not associated with Boundaries are ignored) + - The default boundary condition is periodic + - Inlet boundary condition (non-periodic) can be assigned for + exactly one exterior boundary face, and the opposite boundary face + is then expected to be the outlet (with non-periodic condition) + - Allowed inlet/outlet boundary conditions are + Dirichlet with data VARIABLE = PRESSURE + (where also the in-plane velocity components are fixed) or + Dirichlet with data VARIABLE = VELOCITY + - Note that + Dirichlet conditions must be assigned for the fluid material, + only one Dirichlet condition per boundary is acknowledged, and + the boundary condition data is read directly from the flow field + (i.e. from the data of associated lattice nodes) + """ + + def __init__(self, **kwargs): + """Constructor.""" + # Solver and simulation data + self._solver = None + self._proxy_lat = None + self._is_init_finalized = False + + self._phase = None + self._den = None + self._vel = None + self._acc = None + + # CFD model + self._is_incompressible = True + self._is_singlephase = True + self._is_isothermal = True + self._is_newtonian = True + self._is_laminar = True + self._cfd = None + + # Gravity + self._g = np.zeros((3), dtype=np.float64) + + # Materials + self._fluid_mat = None + self._solid_mat = None + self._material_cnt = 0 + self._kvisc = 1.0/6.0 + self._ref_den = 1.0 + self._inv_ref_den = 1.0 + + # Lattice + self._lat = None + self._lat_size = np.zeros((3), dtype=np.int32) + self._lat_fnode_cnt = 0 # number of fluid lattice sites + self._dr = 0 # lattice spacing + self._inv_dr = 0 + self._V = 0 # volume of the lattice cell + self._inv_V = 0 + + # Solver parameters + self._time_integr = None + self._coll_oper = BGK + self._dt = 0 # discrete time step + self._inv_dt = 0 + self._cr = 0 # lattice reference velocity cr = dr/dt + self._inv_cr = 0.0 + + # Boundary conditions + self._bc = np.full((6), PERIODIC, dtype=np.int32) + + # Call the base class in order to load CUDS + super(Wrapper, self).__init__(**kwargs) + + def _update_gravity(self): + """Read acceleration due to gravity from the CFD model.""" + self._g[0] = self._cfd.gravity_model.acceleration[0] + self._g[1] = self._cfd.gravity_model.acceleration[1] + self._g[2] = self._cfd.gravity_model.acceleration[2] + + def _check_cuds_model(self, cuds): + """Check physics model in the given cuds.""" + print 'Physics model' + print '-------------' + + # General CFD model + cfd_cnt = cuds.count_of(CUBA.CFD) + self._cfd = next(cuds.iter(item_type=CUBA.CFD), None) + + if cfd_cnt != 1: + err_str = 'Expecting exactly one CFD physics model ' + err_str += '(not {0:d})!'.format(cfd_cnt) + raise ValueError(err_str) + + # Electrostatic submodel (ignored) + print '* Ignoring CFD electrostatic model' + + # Gravity submodel + self._update_gravity() + + # Compressibility submodel + self._is_incompressible = True + if not isinstance(self._cfd.compressibility_model, + api.IncompressibleFluidModel): + err_str = 'Compressibility models are not yet supported!' + raise NotImplementedError(err_str) + + # Multiphase submodel + self._is_singlephase = True + if not isinstance(self._cfd.multiphase_model, + api.SinglePhaseModel): + err_str = 'Multiphase models are not yet supported!' + raise NotImplementedError(err_str) + + # Thermal submodel + self._is_isothermal = True + if not isinstance(self._cfd.thermal_model, + api.IsothermalModel): + err_str = 'Thermal models are not yet supported!' + raise NotImplementedError(err_str) + + # Rheology submodel + self._is_newtonian = True + if not isinstance(self._cfd.rheology_model, + api.NewtonianFluidModel): + err_str = 'Rheology models are not yet supported!' + raise NotImplementedError(err_str) + + # Turbulence submodel + self._is_laminar = True + if not isinstance(self._cfd.turbulence_model, + api.LaminarFlowModel): + err_str = 'Turbulence models are not yet supported!' + raise NotImplementedError(err_str) + + print '* Single phase, incompressible, isothermal, '\ + 'laminar, newtonian' + print '* Gravity: ({0:.4e}, {1:.4e}, {2:.4e})'.format(self._g[X], + self._g[Y], + self._g[Z]) + + def _check_cuds_material(self, cuds): + """Check materials in the given cuds.""" + print '---------' + print 'Materials' + print '---------' + + # Fluid and solid (or wall) material + self._fluid_mat = None + self._solid_mat = None + self._material_cnt = cuds.count_of(CUBA.MATERIAL) + + if self._material_cnt < 1 or self._material_cnt > 2: + err_str = 'Expecting exactly 1 or 2 materials ' + err_str += '(not {0:d})!'.format(self._material_cnt) + raise ValueError(err_str) + + mat_ind = 1 + for imat in cuds.iter(item_type=CUBA.MATERIAL): + mat_name = imat.name.lower() + if mat_name == 'solid' or mat_name == 'wall': + self._solid_mat = imat + else: + self._fluid_mat = imat + + print '* Material {0:d}: {1:s}'.format(mat_ind, mat_name) + mat_ind += 1 + + if self._fluid_mat is None: + err_str = 'Expecting one material which is ' + err_str += 'not named as solid or wall!' + raise ValueError(err_str) + + if self._material_cnt == 2 and self._solid_mat is None: + err_str = 'With 2 materials the other must be ' + err_str += 'named as solid or wall!' + raise ValueError(err_str) + + # Density and kinematic viscosity for the fluid material + if CUBA.DENSITY in self._fluid_mat.data: + self._ref_den = self._fluid_mat.data[CUBA.DENSITY] + else: # try + self._ref_den = self._fluid_mat.data[CUBA.REFERENCE_DENSITY] + + self._inv_ref_den = 1.0/self._ref_den + + if CUBA.KINEMATIC_VISCOSITY in self._fluid_mat.data: + self._kvisc = self._fluid_mat.data[CUBA.KINEMATIC_VISCOSITY] + elif CUBA.DYNAMIC_VISCOSITY in self._fluid_mat.data: + visc = self._fluid_mat.data[CUBA.DYNAMIC_VISCOSITY] + self._kvisc = visc/self._ref_den + else: # try + print '* Trying to interpret VISCOSITY as dynamic viscosity' + + visc = self._fluid_mat.data[CUBA.VISCOSITY] + self._kvisc = visc/self._ref_den + + print '* Fluid density: {0:.4e}'.format(self._ref_den) + print '* Fluid kinematic viscosity: {0:.4e}'.format(self._kvisc) + + def _check_cuds_dataset(self, cuds): + """Check dataset in the given cuds.""" + print '-------' + print 'Dataset' + print '-------' + + # Cubic Bravais lattice + cubic_lat_cnt = 0 + self._lat = None + + for ic in cuds.iter(): + if isinstance(ic, ABCLattice): + lat_type = ic.primitive_cell.bravais_lattice + if lat_type is BravaisLattice.CUBIC: + cubic_lat_cnt += 1 + self._lat = ic + + if cubic_lat_cnt != 1: + err_str = 'Expecting exactly one cubic Bravais lattice dataset ' + err_str += '(not {0:d})!'.format(cubic_lat_cnt) + raise ValueError(err_str) + + self._lat_size[0] = self._lat.size[0] + self._lat_size[1] = self._lat.size[1] + self._lat_size[2] = self._lat.size[2] + + p1x = self._lat.primitive_cell.p1[X] + p1y = self._lat.primitive_cell.p1[Y] + p1z = self._lat.primitive_cell.p1[Z] + + dr = np.sqrt(p1x*p1x + p1y*p1y + p1z*p1z) + self._dr = dr + self._inv_dr = 1.0/self._dr + + self._V = self._lat.primitive_cell.volume + self._inv_V = 1.0/self._V + + print '* Cubic lattice: {0:s}'.format(self._lat.name) + print '* Lattice size: {0:}'.format(self._lat.size) + print '* Lattice spacing: {0:.4e}'.format(self._dr) + + def _check_cuds_solver_parameter(self, cuds): + """Check solver parameters in the given cuds.""" + print '-----------------' + print 'Solver parameters' + print '-----------------' + + # Time integration + tint_cnt = cuds.count_of(CUBA.INTEGRATION_TIME) + self._time_integr = next(cuds.iter(item_type=CUBA.INTEGRATION_TIME), + None) + + if tint_cnt != 1: + err_str = 'JYU-LB wrapper (load cuds): must configure exactly ' + err_str += 'one time integration (not {0:d})!'.format(tint_cnt) + raise ValueError(err_str) + + self._dt = self._time_integr.size + self._inv_dt = 1.0/self._dt + self._cr = self._dr/self._dt + self._inv_cr = 1.0/self._cr + + print '* Discrete time step: {0:.4e}'.format(self._dt) + print '* Lattice reference velocity: {0:.4e}'.format(self._cr) + + # Check collision operator; + # related to time integration from the CFD point of view, + # i.e. here NOT part of the (kinetic) model equation specification + coll_oper_cnt = 0 + coll_oper_str = 'bgk' + self._coll_oper = BGK # default + + for sp in cuds.iter(item_type=CUBA.SOLVER_PARAMETER): + if CUBA.COLLISION_OPERATOR in sp.data: + aux_str = sp.data[CUBA.COLLISION_OPERATOR].lower() + coll_oper_cnt += 1 + + if aux_str == 'bgk': + self._coll_oper = BGK + coll_oper_str = aux_str + elif aux_str == 'trt': + self._coll_oper = TRT + coll_oper_str = aux_str + elif aux_str == 'regularization': + self._coll_oper = REGULARIZATION + coll_oper_str = aux_str + else: + wrn_str = 'Ignoring the unknown collision operator ' + wrn_str += '({0:s})'.format(aux_str) + print wrn_str + + if coll_oper_cnt > 1: + wrn_str = ' Several ({0:d}) '.format(coll_oper_cnt) + wrn_str += 'collision operators specified' + print wrn_str + + print '* Using collision operator: {0:s}'.format(coll_oper_str) + + def _check_cuds_boundary_condition(self, cuds): + """Check boundary conditions in the given cuds.""" + print '-------------------' + print 'Boundary conditions' + print '-------------------' + + # Unit outward normals for + # the exterior boundary faces of the simulation domain + bface_uon_str = np.empty((6), dtype="S10") + bface_uon_str[FACE_1X] = '(-1,0,0)' + bface_uon_str[FACE_X1] = '(1,0,0)' + bface_uon_str[FACE_1Y] = '(0,-1,0)' + bface_uon_str[FACE_Y1] = '(0,1,0)' + bface_uon_str[FACE_1Z] = '(0,0,-1)' + bface_uon_str[FACE_Z1] = '(0,0,1)' + + bface_cnt = np.zeros((6), dtype=np.int32) + bface = np.full((6), None, dtype=object) + + # Check the boundaries + for b in cuds.iter(item_type=CUBA.BOUNDARY): + # Check the unit outward normal (aligned with cart.coord.axes) + if CUBA.DIRECTION in b.data: + nx = b.data[CUBA.DIRECTION][0] + ny = b.data[CUBA.DIRECTION][1] + nz = b.data[CUBA.DIRECTION][2] + nx2 = nx*nx + ny2 = ny*ny + nz2 = nz*nz + n2 = nx2 + ny2 + nz2 + + # Not a unit vector + if np.fabs(n2-1.0) > 1e-06: + continue + + # Not aligned with cart.coord.axes + if (np.fabs(nx2-1.0) > 1e-06 and np.fabs(ny2-1.0) > 1e-06 and + np.fabs(nz2-1.0) > 1e-06): + continue + + if nx < 0.0: + bface_cnt[FACE_1X] += 1 + bface[FACE_1X] = b + elif nx > 0.0: + bface_cnt[FACE_X1] += 1 + bface[FACE_X1] = b + elif ny < 0.0: + bface_cnt[FACE_1Y] += 1 + bface[FACE_1Y] = b + elif ny > 0.0: + bface_cnt[FACE_Y1] += 1 + bface[FACE_Y1] = b + elif nz < 0.0: + bface_cnt[FACE_1Z] += 1 + bface[FACE_1Z] = b + else: # nz > 0.0 + bface_cnt[FACE_Z1] += 1 + bface[FACE_Z1] = b + + for f in range(6): + msg_str = '* Exterior boundary face with unit outward normal ' + msg_str += '{0:s}'.format(bface_uon_str[f]) + print msg_str + + bc_str = 'Periodic' + bc = PERIODIC + + if bface_cnt[f] == 0: + print ' Not configured' + elif bface_cnt[f] > 1: + err_str = 'Found more than one ({0:d})!'.format(bface_cnt[f]) + raise ValueError(err_str) + else: + b = bface[f] + bc_cnt = len(b.condition) + + if bc_cnt == 0: + print ' No conditions configured' + elif bc_cnt > 1: + wrn_str = ' Expected 0 or 1 conditions ' + wrn_str += '(not {0:d})'.format(bc_cnt) + print wrn_str + + for c in b.condition: + if isinstance(c, api.Periodic): + bc_str = 'Periodic' + bc = PERIODIC + elif isinstance(c, api.Dirichlet): + c_mat = c.data[CUBA.MATERIAL] + if c_mat != self._fluid_mat: + err_str = 'Dirichlet condition must be specified' + err_str += ' for the fluid material (not for ' + err_str += '{0:s})!'.format(c_mat.name) + raise ValueError(err_str) + + if c.data[CUBA.VARIABLE] == CUBA.PRESSURE: + bc = FIXED_DEN + bc_str = 'Dirichlet (pressure ' + bc_str += '+ in-plane velocity components)' + elif c.data[CUBA.VARIABLE] == CUBA.VELOCITY: + bc = FIXED_VEL + bc_str = 'Dirichlet (velocity)' + else: + wrn_str = ' Dirichlet bc supported only for ' + wrn_str += 'pressure or velocity (not for ' + wrn_str += '{0:})'.format(c.data[CUBA.VARIABLE]) + print wrn_str + wrn_str = ' Ignoring Dirichlet bc (' + wrn_str += 'name: {0:s})'.format(c.name) + print wrn_str + else: + wrn_str = ' Found unsupported bc (' + wrn_str += 'name: {0:s}, '.format(c.name) + wrn_str += 'cuba key: {0:}), '.format(c.cuba_key) + wrn_str += 'ignored' + print wrn_str + + self._bc[f] = bc + print ' Using {0:s} bc'.format(bc_str) + + def _load_cuds(self): + """Load CUDS data for the configuration of the JYU-LB engine.""" + print 'JYU-LB wrapper, load cuds' + print '-------------------------' + + cuds = self.get_cuds() + if cuds is None: + raise ValueError('Unable to load cuds!') + + # Check the cuds for consistency and extract simulation parameters + self._check_cuds_model(cuds) + self._check_cuds_material(cuds) + self._check_cuds_dataset(cuds) + self._check_cuds_solver_parameter(cuds) + self._check_cuds_boundary_condition(cuds) + + # Set up the solver + print '-------------' + print 'Solver set-up' + print '-------------' + kvisc_redu = self._inv_dr*self._inv_cr*self._kvisc + g_redu = np.zeros((3), dtype=np.float64) + g_redu[X] = self._dt*self._inv_cr*self._g[X] + g_redu[Y] = self._dt*self._inv_cr*self._g[Y] + g_redu[Z] = self._dt*self._inv_cr*self._g[Z] + + msg_str = '* Kinematic viscosity (reduced units): ' + msg_str += '{0:.4e}'.format(kvisc_redu) + print msg_str + msg_str = '* Gravity (reduced units): ' + msg_str += '({0:.4e}, {1:.4e}, {2:.4e})'.format(g_redu[X], + g_redu[Y], + g_redu[Z]) + print msg_str + + self._solver = Solver(self._lat_size, self._coll_oper, + kvisc_redu, g_redu, self._bc) + + if self._solver.error_occured(): + raise RuntimeError(self._solver.error_msg) + + # Simulation data + self._den = self._solver.get_den() # density field updated by solver + self._vel = self._solver.get_vel() # vel.field updated by solver + self._acc = self._solver.get_acc() # accel.field utilizied by solver + self._phase = self._solver.get_phase() # phase field util. by solver + + # Initialize simulation data for the lattice nodes + # First initialize the phase information (geometry) + fluid_node_cnt = 0 + for node in self._lat.iter(): + ijk = node.index + + if CUBA.MATERIAL in node.data: + mat = node.data[CUBA.MATERIAL] + if mat.uid == self._fluid_mat.uid: + self._phase[ijk] = FLUID + fluid_node_cnt += 1 + elif mat.uid == self._solid_mat.uid: + self._phase[ijk] = SOLID + else: # unidentified material + err_str = 'JYU-LB wrapper (load cuds): unidentified ' + err_str += 'material (name {0:s}) '.format(mat.name) + err_str += 'associated with the lattice node ' + err_str += '{0}!'.format(ijk) + raise ValueError(err_str) + else: # default material is fluid + self._phase[ijk] = FLUID + fluid_node_cnt += 1 + + print '* Fluid lattice nodes: {0:d}'.format(fluid_node_cnt) + + self._proxy_lat = ProxyLattice(self._lat, self) + self._proxy_lat.update(self._lat.iter()) # Then init. the flow field + + # Replace the original input lattice + # with the newly created proxy lattice + cuds.remove([self._lat.uid]) + cuds.add([self._proxy_lat]) + + self._is_init_finalized = False + self._lat = None # the input dataset/lattice not needed anymore + + def run(self): + """Run modeling engine using the configured settings.""" + # Solver initialization finalized before the first run + if not self._is_init_finalized: + self._lat_fnode_cnt = self._solver.finalize_init() + self._is_init_finalized = True + + # Update particuler simulation parameters before each run + # (bc data updated internally by the solver) + self._update_gravity() + gx_redu = self._dt*self._inv_cr*self._g[X] + gy_redu = self._dt*self._inv_cr*self._g[Y] + gz_redu = self._dt*self._inv_cr*self._g[Z] + self._solver.gravity = (gx_redu, gy_redu, gz_redu) + + tperiod = self._time_integr.final - self._time_integr.current + tsteps = int(tperiod/self._dt) + tsimulated = tsteps*self._dt + + self._solver.evolve(tsteps) + + # Advance the integration times after each run + self._time_integr.current += tsimulated + self._time_integr.final += tsimulated + + def add_dataset(self, container): + """Add a CUDS Lattice container. + + Parameters + ---------- + container : {ABCLattice} + The CUDS Lattice container to add to the engine. + + Raises + ------ + NotImplementedError: + Allways (no support for this feature). + """ + err_str = 'JYU-LB wrapper (add dataset): adding datasets ' + err_str += 'directly to the wrapper not supported!' + raise NotImplementedError(err_str) + + def remove_dataset(self, name): + """Delete a lattice. + + Parameters + ---------- + name : str + name of the lattice to be deleted. + + Raises + ------ + NotImplementedError: + Allways (no support for this feature). + """ + err_str = 'JYU-LB wrapper (add dataset): removing datasets ' + err_str += 'from the wrapper not supported!' + raise NotImplementedError(err_str) + + def get_dataset(self, name): + """ Get a lattice. + + The returned lattice can be used to query and update the state of the + lattice inside the modeling engine. + + Parameters + ---------- + name : str + name of the lattice to be retrieved + (not used, the proxy lattice is always returned). + + Returns + ------- + ABCLattice + """ + return self._proxy_lattice + + def get_dataset_names(self): + """ Return the names of all datasets in the engine workspace.""" + return [self._proxy_lattice.name] + + def iter_datasets(self, names=None): + """Iterate over a subset or all of the lattices. + + Parameters + ---------- + names : sequence of str, optional + names of specific lattices to be iterated over. If names is not + given, then all lattices will be iterated over. + + Returns + ------- + A generator of ABCLattice objects + + Raises + ------ + NotImplementedError: + Allways (no support for this feature; forces to use get_dataset). + """ + err_str = 'JYU-LB wrapper (iter.dataset): iterating over datasets is' + err_str += ' not supported (use get_dataset to access the single ' + err_str += 'dataset)!' + raise NotImplementedError(err_str) diff --git a/wrapper/jyulb_lattice.py b/wrapper/jyulb_lattice.py new file mode 100644 index 0000000..3398294 --- /dev/null +++ b/wrapper/jyulb_lattice.py @@ -0,0 +1,269 @@ +"""A proxy lattice for accessing simulation data of a JYU-LB engine.""" +import uuid +import numpy as np + +from simphony.core.data_container import DataContainer +from simphony.core.cuba import CUBA + +from simphony.cuds.abc_lattice import ABCLattice +from simphony.cuds.lattice_items import LatticeNode +from simphony.cuds.primitive_cell import PrimitiveCell + +from jyulb.defs import X, Y, Z, SOLID + + +class ProxyLattice(ABCLattice): + + """Access to the simulation data of a JYU-LB engine. + + Lattice node data + ================= + Provided by get + --------------- + MATERIAL, PRESSURE, DENSITY, VELOCITY, ACCELERATION + + Acknowledged by update + ---------------------- + PRESSURE (or alternatively DENSITY), + VELOCITY (or alternatively MOMENTUM), + ACCELERATION (or alternatively FORCE) + """ + + cuba_key = CUBA.LATTICE + + def __init__(self, src_lat, jyulb_wrapper): + """Constructor.""" + # Copy data from the source lattice + self.name = src_lat.name + + pc = src_lat._primitive_cell + self._primitive_cell = PrimitiveCell(pc.p1, pc.p2, pc.p3, + pc.bravais_lattice) + + self._size = np.array(src_lat.size, dtype=np.int32) + self._origin = np.array(src_lat.origin, dtype=np.float) + + self._data = DataContainer(src_lat.data) + + self._items_count = { + CUBA.NODE: lambda: self._size + } + self._uid = uuid.uuid4() + + # References to the simulation data + self._phase = jyulb_wrapper._phase + self._den = jyulb_wrapper._den + self._vel = jyulb_wrapper._vel + self._acc = jyulb_wrapper._acc + + # Parameters for unit transformations + # (SI <-> reduced unit system) + self._ref_den = jyulb_wrapper._ref_den + self._inv_ref_den = 1.0/self._ref_den + + self._dr = jyulb_wrapper._dr + self._inv_dr = 1.0/self._dr + + self._dt = jyulb_wrapper._dt + self._inv_dt = 1.0/self._dt + + self._cr = self._dr/self._dt + self._inv_cr = 1.0/self._cr + + # References to the fluid and solid materials + self._fluid_mat = jyulb_wrapper._fluid_mat + self._solid_mat = jyulb_wrapper._solid_mat + + # Reference to a JYU-LB engine + self._solver = jyulb_wrapper._solver + + @property + def uid(self): + """Universal identifier.""" + return self._uid + + def count_of(self, item_type): + """ Return the count of item_type in the container. + + Parameters + ---------- + item_type : CUBA + The CUBA enum of the type of the items to return + the count of. + + Returns + ------- + count : int + The number of items of item_type in the container. + + Raises + ------ + ValueError : + If the type of the item is not supported in the current + container. + + """ + try: + return np.prod(self._items_count[item_type]()) + except KeyError: + error_str = "Trying to obtain count a of non-supported item: {}" + raise ValueError(error_str.format(item_type)) + + @property + def size(self): + """Access (get).""" + return self._size + + @property + def origin(self): + """Access (get).""" + return self._origin + + @property + def data(self): + """Access (get).""" + return self._data + + @data.setter + def data(self, value): + """Access (set).""" + self._data = DataContainer(value) + + def _get_node(self, index): + """Get a copy of the node corresponding to the given index. + + Parameters + ---------- + index : int[3] + node index coordinate + + Returns + ------- + A reference to a LatticeNode object + + """ + ijk = tuple(index) + ijkx = ijk + (X,) + ijky = ijk + (Y,) + ijkz = ijk + (Z,) + + if ijk[0] < 0 or ijk[1] < 0 or ijk[2] < 0: + raise IndexError('invalid index: {}'.format(ijk)) + + node = LatticeNode(index) + + if self._phase[ijk] == SOLID: + node.data[CUBA.MATERIAL] = self._solid_mat + else: + node.data[CUBA.MATERIAL] = self._fluid_mat + + # Transformation from the reduced unit system to SI-units + den = self._den[ijk] # reduced units + p = self._solver.eos_den2p(den) # reduced units + p *= self._ref_den*self._cr*self._cr # SI + node.data[CUBA.PRESSURE] = p + node.data[CUBA.DENSITY] = self._ref_den*den + + vx = self._vel[ijkx] # reduced units + vy = self._vel[ijky] # reduced units + vz = self._vel[ijkz] # reduced units + vx *= self._cr # SI + vy *= self._cr # SI + vz *= self._cr # SI + node.data[CUBA.VELOCITY] = (vx, vy, vz) + + ax = self._acc[ijkx] # reduced units + ay = self._acc[ijky] # reduced units + az = self._acc[ijkz] # reduced units + ax *= self._cr*self._inv_dt # SI + ay *= self._cr*self._inv_dt # SI + az *= self._cr*self._inv_dt # SI + node.data[CUBA.ACCELERATION] = (ax, ay, az) + + return node + + def _update_nodes(self, nodes): + """Update the corresponding lattice nodes (data copied). + + Parameters + ---------- + nodes : iterable of LatticeNode objects + reference to LatticeNode objects from where the data is copied + to the Lattice + """ + for node in nodes: + ijk = tuple(node.index) + ijkx = ijk + (X,) + ijky = ijk + (Y,) + ijkz = ijk + (Z,) + + if any(value < 0 for value in ijk): + raise IndexError('invalid index: {}'.format(ijk)) + + # Transformation from SI-units to the reduced unit system + if CUBA.PRESSURE in node.data: + p = node.data[CUBA.PRESSURE] # in SI-units + p *= self._inv_ref_den*self._inv_cr*self._inv_cr + self._den[ijk] = self._solver.eos_p2den(p) + elif CUBA.DENSITY in node.data: + self._den[ijk] = self._inv_ref_den*node.data[CUBA.DENSITY] + + if CUBA.VELOCITY in node.data: + self._vel[ijkx] = self._inv_cr*node.data[CUBA.VELOCITY][X] + self._vel[ijky] = self._inv_cr*node.data[CUBA.VELOCITY][Y] + self._vel[ijkz] = self._inv_cr*node.data[CUBA.VELOCITY][Z] + elif CUBA.MOMENTUM in node.data: + # interpreted as momentum density + den = self._ref_den*self._den[ijk] # in SI-units + inv_den = 1.0/den # in SI-units + vx = inv_den*node.data[CUBA.MOMENTUM][X] # in SI-units + vy = inv_den*node.data[CUBA.MOMENTUM][Y] # in SI-units + vz = inv_den*node.data[CUBA.MOMENTUM][Z] # in SI-units + self._vel[ijkx] = self._inv_cr*vx + self._vel[ijky] = self._inv_cr*vy + self._vel[ijkz] = self._inv_cr*vz + + if CUBA.ACCELERATION in node.data: + ax = node.data[CUBA.ACCELERATION][X] # in SI-units + ay = node.data[CUBA.ACCELERATION][Y] # in SI-units + az = node.data[CUBA.ACCELERATION][Z] # in SI-units + self._acc[ijkx] = self._inv_cr*self._dt*ax + self._acc[ijky] = self._inv_cr*self._dt*ay + self._acc[ijkz] = self._inv_cr*self._dt*az + elif CUBA.FORCE in node.data: + # interpreted as body force (i.e. force density) + den = self._ref_den*self._den[ijk] # in SI-units + inv_den = 1.0/den # in SI-units + ax = inv_den*node.data[CUBA.FORCE][X] # in SI-units + ay = inv_den*node.data[CUBA.FORCE][Y] # in SI-units + az = inv_den*node.data[CUBA.FORCE][Z] # in SI-units + self._acc[ijkx] = self._inv_cr*self._dt*ax + self._acc[ijky] = self._inv_cr*self._dt*ay + self._acc[ijkz] = self._inv_cr*self._dt*az + + def _iter_nodes(self, indices=None): + """Get an iterator over the LatticeNodes described by the indices. + + Parameters + ---------- + indices : iterable set of int[3], optional + When indices (i.e. node index coordinates) are provided, then + nodes are returned in the same order of the provided indices. + If indices is None, there is no restriction on the order of the + returned nodes. + + Returns + ------- + A generator for LatticeNode objects + + """ + sx = self.size[X] + sy = self.size[Y] + sz = self.size[Z] + + if indices is None: + for index in np.ndindex((sx, sy, sz)): + yield self.get(index) + else: + for index in indices: + yield self.get(index) diff --git a/wrapper/tests/__init__.py b/wrapper/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/wrapper/tests/test_jyulb_engine.py b/wrapper/tests/test_jyulb_engine.py new file mode 100644 index 0000000..e3d2988 --- /dev/null +++ b/wrapper/tests/test_jyulb_engine.py @@ -0,0 +1,219 @@ +"""Testing module for the JYU-LB wrapper.""" +import os +import shutil +import tempfile +import unittest +import numpy as np + +from simphony.cuds.meta import api +from simphony.core.cuba import CUBA +from simphony.api import CUDS, Simulation +from simphony.engine import EngineInterface +from simphony.cuds.lattice import make_cubic_lattice +from jyulb.flow_field import poise_plate_steady_state_pgrad + + +class JYULBEngineTestCaseSimulate(unittest.TestCase): + + """Simulation test for checking that installation is ok.""" + + def setUp(self): + """Initialize the test case.""" + self.temp_dir = tempfile.mkdtemp() + self.saved_path = os.getcwd() + + os.chdir(self.temp_dir) + self.addCleanup(self.cleanup) + + def cleanup(self): + """Clean files after the test case.""" + os.chdir(self.saved_path) + shutil.rmtree(self.temp_dir) + + def cuds_factory(self): + """Create CUDS for poiseuille fluid flow simulation.""" + # CUDS + cuds = CUDS(name='poiseuille') + + # Physics model + cfd = api.Cfd(name='fluid flow') + + # Submodels (in fact, these are the default submodels in CFD) + cfd.thermal_model = api.IsothermalModel(name='isothermal') + cfd.turbulence_model = api.LaminarFlowModel(name='laminar') + cfd.rheology_model = api.NewtonianFluidModel(name='newtonian') + cfd.multiphase_model = api.SinglePhaseModel(name='singlephase') + compr_model = api.IncompressibleFluidModel(name='incompressible') + cfd.compressibility_model = compr_model + + # Gravity in the z-direction + cfd.gravity_model.acceleration = (0, 0, 1e-5) + cuds.add([cfd]) + + # Materials + fluid = api.Material(name='fluid') + solid = api.Material(name='wall') + fluid.data[CUBA.DENSITY] = 1.0 + fluid.data[CUBA.KINEMATIC_VISCOSITY] = 2.0/6.0 + cuds.add([fluid, solid]) + + # Dataset (lattice) + lat_size = (5, 1, 1) # number of lattice nodes + lat_h = 1.0 # lattice spacing + lat = make_cubic_lattice('channel', lat_h, lat_size) + + # Create a channel geometry, set the initial flow field + for node in lat.iter(): + ijk = node.index + # The two plates are separated in the x-direction + if ijk[0] == 0 or ijk[0] == (lat_size[0]-1): + node.data[CUBA.MATERIAL] = solid + node.data[CUBA.VELOCITY] = (0, 0, 0) + else: + node.data[CUBA.MATERIAL] = fluid + node.data[CUBA.VELOCITY] = (0, 0, 0) + + lat.update([node]) + + cuds.add([lat]) + + # Boundary conditions + # for the exterior boundary faces of the simulation domain + # (the directions refer to unit outward normals) + face1y = api.Boundary(name='face1y', condition=[api.Periodic()]) + face1y.data[CUBA.DIRECTION] = (0, -1, 0) + + facey1 = api.Boundary(name='facey1', condition=[api.Periodic()]) + facey1.data[CUBA.DIRECTION] = (0, 1, 0) + + face1z = api.Boundary(name='face1z', condition=[api.Periodic()]) + face1z.data[CUBA.DIRECTION] = (0, 0, -1) + + facez1 = api.Boundary(name='facez1', condition=[api.Periodic()]) + facez1.data[CUBA.DIRECTION] = (0, 0, 1) + + cuds.add([face1y, facey1, face1z, facez1]) + + # Solver parameters (time integration, collision operator) + ti = api.IntegrationTime(name='simulation time', + current=0.0, final=100.0, size=1.0) + + ti.data[CUBA.COLLISION_OPERATOR] = 'TRT' + cuds.add([ti]) + + self._fluid = fluid + self._solid = solid + self._cuds = cuds + self._cfd = cfd + self._lat = lat + self._ti = ti + + return cuds + + def analyse_results(self): + """Compute the relative L2-error norm.""" + H = self._lat.size[0]-2 # plates separated in the x-direction + + den = self._fluid.data[CUBA.DENSITY] + dvisc = den*self._fluid.data[CUBA.KINEMATIC_VISCOSITY] + eff_pgrad = -den*self._cfd.gravity_model.acceleration[2] + + centrel_dist = np.zeros((H), dtype=np.float64) + ana_vel = np.zeros((H), dtype=np.float64) + + poise_plate_steady_state_pgrad(H, dvisc, eff_pgrad, + centrel_dist, ana_vel) + + sim_loc_vel = np.zeros(3, dtype=np.float64) + coord = np.zeros(3, dtype=np.int32) + coord[0] = int(self._lat.size[0]/2) + coord[1] = int(self._lat.size[1]/2) + coord[2] = int(self._lat.size[2]/2) + + ana2 = 0.0 + sim_ana_diff2 = 0.0 + sim_loc_vel = np.zeros((3), dtype=np.float64) + + for h in range(H): + coord[0] = h+1 + node = self._lat.get(coord) + + if node.data[CUBA.MATERIAL] == self._solid: + continue + + sim_loc_vel[0] = node.data[CUBA.VELOCITY][0] + sim_loc_vel[1] = node.data[CUBA.VELOCITY][1] + sim_loc_vel[2] = node.data[CUBA.VELOCITY][2] + + sim_loc_speed = np.sqrt(sim_loc_vel[0]*sim_loc_vel[0] + + sim_loc_vel[1]*sim_loc_vel[1] + + sim_loc_vel[2]*sim_loc_vel[2]) + + ana_loc_vel = ana_vel[h] + + sim_ana_diff = sim_loc_speed - ana_loc_vel + sim_ana_diff2 += sim_ana_diff*sim_ana_diff + ana2 += ana_loc_vel*ana_loc_vel + + rel_l2_err_norm_vel = np.sqrt(sim_ana_diff2/ana2) + + print 'Relative L2-error norm: {0:.4e}'.format(rel_l2_err_norm_vel) + print '-'*77 + + return rel_l2_err_norm_vel + + def test_run_engine(self): + """Running the JYU-LB engine.""" + cuds = self.cuds_factory() + + sim = Simulation(cuds, 'JYU-LB', + engine_interface=EngineInterface.Internal) + + lat = cuds.get_by_name(self._lat.name) + self._lat = lat + + # Total fluid mass before simulation + tot_fmass = 0.0 + for node in lat.iter(): + if node.data[CUBA.MATERIAL] == self._fluid: + tot_fmass += node.data[CUBA.DENSITY] + + print '='*77 + print 'Poiseuille flow simulation' + print '='*77 + print "Total fluid mass before simulation: {0:.4e}".format(tot_fmass) + print '-'*77 + + # Simulate + sim.run() + + # Total velocity after simulation + tot_vx = 0.0 + tot_vy = 0.0 + tot_vz = 0.0 + for node in lat.iter(): + if node.data[CUBA.MATERIAL] == self._fluid: + tot_vx += node.data[CUBA.VELOCITY][0] + tot_vy += node.data[CUBA.VELOCITY][1] + tot_vz += node.data[CUBA.VELOCITY][2] + + sf = 'Time = {0:f}, tot.vel = ({1:11.4e}, {2:11.4e}, {3:11.4e})' + print sf.format(self._ti.current, tot_vx, tot_vy, tot_vz) + + # Total fluid mass after simulation + tot_fmass = 0.0 + for node in lat.iter(): + if node.data[CUBA.MATERIAL] == self._fluid: + tot_fmass += node.data[CUBA.DENSITY] + + print '-'*77 + print "Total fluid mass after simulation: {0:.4e}".format(tot_fmass) + print '-'*77 + + rel_l2_error = self.analyse_results() + + self.assertTrue(rel_l2_error < 1.0e-10) + + +if __name__ == '__main__': + unittest.main()