diff --git a/source/pysph/base/kernels.h b/source/pysph/base/kernels.h index cfdb18c..a1f48ca 100644 --- a/source/pysph/base/kernels.h +++ b/source/pysph/base/kernels.h @@ -86,6 +86,7 @@ void cubic_spline_gradient(REAL4 pa, REAL4 pb, REAL4* grad, grad[0].s0 = val * rab.s0; grad[0].s1 = val * rab.s1; grad[0].s2 = val * rab.s2; + grad[0].s3 = 0.0F; } diff --git a/source/pysph/sph/funcs/energy_funcs.clt b/source/pysph/sph/funcs/energy_funcs.clt index e69de29..daf88f4 100644 --- a/source/pysph/sph/funcs/energy_funcs.clt +++ b/source/pysph/sph/funcs/energy_funcs.clt @@ -0,0 +1,144 @@ + + +$EnergyEquationNoVisc + +/* Thermal energy equation (eq 10) described in "A shock-captiting SPH + scheme based on adaptive kernel estimation" by Leonardo + Di. G. Sigalotti et al. + + The equation without artificial viscosity and artificial heat is + + \frac{DU_a}{Dt} = \frac{1}{2}\sum_{b=1}^{N}m_b\left[ + \left(\frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2}\right)\,(v_a - + v_b)\right]\,\nabla_a \cdot W_{ab} + +*/ + +#include "cl_common.h" +#include "cl_common.cl" +#include "kernels.h" + +__kernel void EnergyEquationNoVisc(%(kernel_args)s) +{ + %(workgroup_code)s + + // The term `dest_id` will be suitably defined at this point. + + REAL4 ra = (REAL4)( d_x[dest_id], d_y[dest_id], + d_z[dest_id], d_h[dest_id] ); + + REAL rhoa = d_rho[dest_id]; + REAL Pa = d_p[dest_id]; + + REAL4 va = (REAL4)(d_u[dest_id], d_v[dest_id], d_w[dest_id], 0.0F); + + REAL4 grad; + REAL temp; + + %(neighbor_loop_code)s + + { + + // SPH innermost loop code goes here. The index `src_id` will + // be available and looped over, this index. + + REAL4 rb = (REAL4)(s_x[src_id], s_y[src_id], s_z[src_id], s_h[src_id]); + REAL Pb = s_p[src_id]; + REAL mb = s_m[src_id]; + REAL rhob = s_rho[src_id]; + + REAL4 vb = (REAL4)(s_u[src_id], s_v[src_id], s_w[src_id], 0.0F); + + temp = 0.5F*mb*( Pa/(rhoa*rhoa) + Pb/(rhob*rhob) ) ; + + kernel_gradient(ra, rb, &grad, dim, kernel_type); + + tmpx[dest_id] += temp * ( dot(va-vb, grad) ); + + } // neighbor loop + +} // __kernel + +$EnergyEquationNoVisc + + +$EnergyEquationWithVisc + +/* Thermal energy equation (eq 10) described in "A shock-captiting SPH + scheme based on adaptive kernel estimation" by Leonardo + Di. G. Sigalotti et al. + + The equation without artificial heat is + + \frac{DU_a}{Dt} = \frac{1}{2}\sum_{b=1}^{N}m_b\left[ + \left(\frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2} + + \Pi_{ab}\right)\,(v_a - v_b)\right]\,\nabla_a \cdot W_{ab} + +*/ + +#include "cl_common.h" +#include "cl_common.cl" +#include "kernels.h" + +__kernel void EnergyEquationWithVisc(%(kernel_args)s) +{ + %(workgroup_code)s + + // The term `dest_id` will be suitably defined at this point. + + REAL4 ra = (REAL4)( d_x[dest_id], d_y[dest_id], + d_z[dest_id], d_h[dest_id] ); + + REAL4 va = (REAL4)(d_u[dest_id], d_v[dest_id], d_w[dest_id], 0.0F); + + REAL rhoa = d_rho[dest_id]; + REAL Pa = d_p[dest_id]; + REAL4 grad; + REAL temp, Pi_ab; + + %(neighbor_loop_code)s + + { + + // SPH innermost loop code goes here. The index `src_id` will + // be available and looped over, this index. + + REAL4 rb = (REAL4)(s_x[src_id], s_y[src_id], s_z[src_id], s_h[src_id]); + REAL Pb = s_p[src_id]; + REAL rhob = s_rho[src_id]; + + REAL4 vb = (REAL4)(s_u[src_id], s_v[src_id], s_w[src_id], 0.0F); + + REAL dot_product = dot( (va-vb), (ra-rb) ); + + temp = Pa/(rhoa*rhoa) + Pb/(rhob*rhob) ; + + kernel_gradient(ra, rb, &grad, dim, kernel_type); + + Pi_ab = 0.0F; + if ( dot_product < 0.0F ) + { + REAL cab = 0.5F * ( d_cs[dest_id] + s_cs[src_id] ); + REAL rhoab = 0.5F * (rhoa + rhob); + + REAL hab = 0.5F * ( d_h[dest_id] + s_h[src_id] ); + REAL mu = dot_product*hab; + REAL norm = length(ra-rb); + + mu /= ( norm*norm + eta*eta*hab*hab ); + + Pi_ab = -alpha*cab*mu + beta*mu*mu; + Pi_ab /= rhoab; + + } //if + + temp += Pi_ab; + temp *= 0.5F*s_m[src_id]; + + tmpx[dest_id] += temp * dot( va-vb, grad ); + + } //neighbor_loop_code + +} //__kernel + +$EnergyEquationWithVisc diff --git a/source/pysph/sph/funcs/energy_funcs.pyx b/source/pysph/sph/funcs/energy_funcs.pyx index 709ec7f..373271d 100644 --- a/source/pysph/sph/funcs/energy_funcs.pyx +++ b/source/pysph/sph/funcs/energy_funcs.pyx @@ -2,6 +2,8 @@ from pysph.base.point cimport cPoint, cPoint_dot, cPoint_new, cPoint_sub,\ cPoint_norm +from pysph.solver.cl_utils import get_real + cdef extern from "math.h": double sqrt(double) double fabs(double) @@ -31,6 +33,9 @@ cdef class EnergyEquationNoVisc(SPHFunctionParticle): self.src_reads.extend( ['u','v','w','p'] ) self.dst_reads.extend( ['u','v','w','p'] ) + + def _set_extra_cl_args(self): + pass cdef void eval_nbr(self, size_t source_pid, size_t dest_pid, KernelBase kernel, double *nr): @@ -89,7 +94,14 @@ cdef class EnergyEquationNoVisc(SPHFunctionParticle): tmp = 0.5*mb*(pa/(rhoa*rhoa) + pb/(rhob*rhob)) nr[0] += tmp*dot -############################################################################## + + def cl_eval(self, object queue, object context, object kernel): + + self.set_cl_kernel_args() + + self.cl_program.EnergyEquationNoVisc( + queue, self.global_sizes, self.local_sizes, *self.cl_args).wait() + ############################################################################## cdef class EnergyEquationAVisc(SPHFunctionParticle): @@ -232,12 +244,13 @@ cdef class EnergyEquation(SPHFunctionParticle): self.alpha = alpha self.beta = beta self.gamma = gamma + self.eta = eta self.id = 'energyequation' self.tag = "energy" self.cl_kernel_src_file = "energy_funcs.cl" - self.cl_kernel_function_name = "EnergyEquation" + self.cl_kernel_function_name = "EnergyEquationWithVisc" def set_src_dst_reads(self): self.src_reads = [] @@ -247,7 +260,20 @@ cdef class EnergyEquation(SPHFunctionParticle): self.dst_reads.extend( ['x','y','z','h','rho','tag'] ) self.src_reads.extend( ['u','v','w','p','cs'] ) - self.dst_reads.extend( ['u','v','w','p','cs'] ) + self.dst_reads.extend( ['u','v','w','p','cs'] ) + + def _set_extra_cl_args(self): + self.cl_args.append( get_real(self.alpha, self.dest.cl_precision) ) + self.cl_args_name.append( 'REAL const alpha' ) + + self.cl_args.append( get_real(self.beta, self.dest.cl_precision) ) + self.cl_args_name.append( 'REAL const beta' ) + + self.cl_args.append( get_real(self.gamma, self.dest.cl_precision) ) + self.cl_args_name.append( 'REAL const gamma' ) + + self.cl_args.append( get_real(self.eta, self.dest.cl_precision) ) + self.cl_args_name.append( 'REAL const eta' ) cdef void eval_nbr(self, size_t source_pid, size_t dest_pid, KernelBase kernel, double *nr): @@ -332,8 +358,14 @@ cdef class EnergyEquation(SPHFunctionParticle): tmp = cPoint_dot(grad, vab) * tmp nr[0] += 0.5*mb*tmp - + def cl_eval(self, object queue, object context, object kernel): + + self.set_cl_kernel_args() + + self.cl_program.EnergyEquationWithVisc( + queue, self.global_sizes, self.local_sizes, *self.cl_args).wait() + ################################################################################ # `ArtificialHeat` class. ################################################################################ diff --git a/source/pysph/sph/funcs/pressure_funcs.clt b/source/pysph/sph/funcs/pressure_funcs.clt index 86faf08..cd3a05f 100644 --- a/source/pysph/sph/funcs/pressure_funcs.clt +++ b/source/pysph/sph/funcs/pressure_funcs.clt @@ -69,7 +69,6 @@ __kernel void MomentumEquation(%(kernel_args)s) REAL4 grad; REAL temp, Pi_ab; - Pi_ab = 0.0F; %(neighbor_loop_code)s { @@ -80,6 +79,7 @@ __kernel void MomentumEquation(%(kernel_args)s) REAL4 rb = (REAL4)(s_x[src_id], s_y[src_id], s_z[src_id], s_h[src_id]); REAL Pb = s_p[src_id]; REAL rhob = s_rho[src_id]; + REAL mb = s_m[src_id]; REAL4 vb = (REAL4)(s_u[src_id], s_v[src_id], s_w[src_id], 0.0F); @@ -88,7 +88,8 @@ __kernel void MomentumEquation(%(kernel_args)s) temp = Pa/(rhoa*rhoa) + Pb/(rhob*rhob) ; kernel_gradient(ra, rb, &grad, dim, kernel_type); - + + Pi_ab = 0.0F; if ( dot_product < 0.0F ) { REAL cab = 0.5F * ( d_cs[dest_id] + s_cs[src_id] ); @@ -96,8 +97,9 @@ __kernel void MomentumEquation(%(kernel_args)s) REAL hab = 0.5F * ( d_h[dest_id] + s_h[src_id] ); REAL mu = dot_product*hab; - - mu /= ( length(ra-rb) + eta*eta*hab*hab ); + REAL norm2 = length(ra-rb) * length(ra-rb); + + mu /= ( norm2 + eta*eta*hab*hab ); Pi_ab = -alpha*cab*mu + beta*mu*mu; Pi_ab /= rhoab; @@ -105,7 +107,7 @@ __kernel void MomentumEquation(%(kernel_args)s) } //if temp += Pi_ab; - temp *= -s_m[src_id]; + temp = -mb*temp; tmpx[dest_id] += temp*grad.x; tmpy[dest_id] += temp*grad.y; diff --git a/source/pysph/sph/tests/test_energy_funcs.py b/source/pysph/sph/tests/test_energy_funcs.py new file mode 100644 index 0000000..eec9af6 --- /dev/null +++ b/source/pysph/sph/tests/test_energy_funcs.py @@ -0,0 +1,334 @@ +""" Tests for the energy force functions """ + +import pysph.base.api as base +import pysph.solver.api as solver +import pysph.sph.api as sph + +if solver.HAS_CL: + import pyopencl as cl + +import numpy +import unittest +from os import path + +NSquareLocator = base.NeighborLocatorType.NSquareNeighborLocator + +class EnergyFunctionsTestCase(unittest.TestCase): + + def runTest(self): + pass + + def setUp(self): + """ The setup consists of four particles placed at the + vertices of a unit square. + + The function tested is + + ..math:: + + \frac{DU_a}{Dt} = \frac{1}{2}\sum_{b=1}^{N}m_b\left[ + \left(\frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2}\right)\,(v_a - + v_b)\right]\,\nabla_a \cdot W_{ab} + + The mass of each particle is 1 + + """ + + self.precision = "single" + + self.np = 4 + + x = numpy.array([0, 0, 1, 1], numpy.float64) + y = numpy.array([0, 1, 1, 0], numpy.float64) + + z = numpy.zeros_like(x) + m = numpy.ones_like(x) + + u = numpy.array([1, 0, 0, -1], numpy.float64) + p = numpy.array([0, 0, 1, 1], numpy.float64) + + tmpx = numpy.zeros_like(x) + tmpy = numpy.zeros_like(x) + tmpz = numpy.zeros_like(x) + + self.pa = pa = base.get_particle_array(name="test", x=x, y=y, z=z, + m=m, u=u, p=p, + tmpx=tmpx, tmpy=tmpy, tmpz=tmpz, + cl_precision=self.precision) + + env = sph.EnergyEquationNoVisc.withargs() + ewv = sph.EnergyEquation.withargs(alpha=1.0, beta=1.0, + gamma=1.4, eta=0.1) + + + self.env = env.get_func(pa,pa) + self.ewv = ewv.get_func(pa,pa) + + self.env.kernel = base.CubicSplineKernel(dim=2) + self.env.nbr_locator = \ + base.Particles.get_neighbor_particle_locator(pa, + pa) + + self.ewv.kernel = base.CubicSplineKernel(dim=2) + self.ewv.nbr_locator = \ + base.Particles.get_neighbor_particle_locator(pa, + pa) + + self.setup_cl() + + def setup_cl(self): + pass + +class EnergyEquationNoViscTestCase(EnergyFunctionsTestCase): + + def setup_cl(self): + pa = self.pa + + if solver.HAS_CL: + self.ctx = ctx = cl.create_some_context() + self.q = q = cl.CommandQueue(ctx) + + pa.setup_cl(ctx, q) + + pysph_root = solver.get_pysph_root() + + template = solver.cl_read( + path.join(pysph_root, "sph/funcs/energy_funcs.clt"), + function_name=self.env.cl_kernel_function_name, + precision=self.precision) + + prog_src = solver.create_program(template, self.env) + + self.prog=cl.Program(ctx, prog_src).build(solver.get_cl_include()) + + def get_reference_solution(self): + """ Evaluate the force on each particle manually """ + + pa = self.pa + forces = [] + + x,y,z,p,m,h,rho = pa.get('x','y','z','p','m','h','rho') + u,v,w = pa.get('u','v','w') + + kernel = base.CubicSplineKernel(dim=2) + + for i in range(self.np): + + force = base.Point() + xa, ya, za = x[i], y[i], z[i] + ua, va, wa = u[i], v[i], w[i] + + ra = base.Point(xa,ya,za) + Va = base.Point(ua,va,wa) + + Pa, rhoa = p[i], rho[i] + ha = h[i] + + for j in range(self.np): + + grad = base.Point() + xb, yb, zb = x[j], y[j], z[j] + ub, vb, wb = u[j], v[j], w[j] + + Pb, rhob = p[j], rho[j] + hb, mb = m[j], h[j] + + havg = 0.5 * (ha + hb) + + rb = base.Point(xb, yb, zb) + Vb = base.Point(ub, vb, wb) + + tmp = 0.5*mb * ( Pa/(rhoa*rhoa) + Pb/(rhob*rhob) ) + kernel.py_gradient(ra, rb, havg, grad) + + force.x += tmp * grad.dot(Va-Vb) + + forces.append(force) + + return forces + + def test_eval(self): + """ Test the PySPH solution """ + + pa = self.pa + func = self.env + + k = base.CubicSplineKernel(dim=2) + + tmpx = pa.properties['tmpx'] + tmpy = pa.properties['tmpy'] + tmpz = pa.properties['tmpz'] + + func.eval(k, tmpx, tmpy, tmpz) + + reference_solution = self.get_reference_solution() + + for i in range(self.np): + self.assertAlmostEqual(reference_solution[i].x, tmpx[i]) + self.assertAlmostEqual(reference_solution[i].y, tmpy[i]) + self.assertAlmostEqual(reference_solution[i].z, tmpz[i]) + + def test_cl_eval(self): + """ Test the PyOpenCL implementation """ + + if solver.HAS_CL: + + pa = self.pa + func = self.env + + k = base.CubicSplineKernel(dim=2) + + func.setup_cl(self.prog, self.ctx) + + func.cl_eval(self.q, self.ctx, k) + + pa.read_from_buffer() + + reference_solution = self.get_reference_solution() + + for i in range(self.np): + self.assertAlmostEqual(reference_solution[i].x, pa.tmpx[i], 6) + self.assertAlmostEqual(reference_solution[i].y, pa.tmpy[i], 6) + self.assertAlmostEqual(reference_solution[i].z, pa.tmpz[i], 6) + +class EnergyEquationTestCase(EnergyFunctionsTestCase): + + def setup_cl(self): + pa = self.pa + + if solver.HAS_CL: + self.ctx = ctx = cl.create_some_context() + self.q = q = cl.CommandQueue(ctx) + + pa.setup_cl(ctx, q) + + pysph_root = solver.get_pysph_root() + + template = solver.cl_read( + path.join(pysph_root, "sph/funcs/energy_funcs.clt"), + function_name=self.ewv.cl_kernel_function_name, + precision=self.precision) + + self.ewv.set_cl_kernel_args() + prog_src = solver.create_program(template, self.ewv) + + self.prog=cl.Program(ctx, prog_src).build(solver.get_cl_include()) + + def get_reference_solution(self): + """ Evaluate the force on each particle manually """ + + pa = self.pa + forces = [] + + x,y,z,p,m,h,rho = pa.get('x','y','z','p','m','h','rho') + u,v,w,cs = pa.get('u','v','w','cs') + + kernel = base.CubicSplineKernel(dim=2) + + for i in range(self.np): + + force = base.Point() + xa, ya, za = x[i], y[i], z[i] + ua, va, wa = u[i], v[i], w[i] + + ra = base.Point(xa,ya,za) + Va = base.Point(ua,va,wa) + + Pa, rhoa = p[i], rho[i] + ha = h[i] + + for j in range(self.np): + + grad = base.Point() + xb, yb, zb = x[j], y[j], z[j] + Pb, rhob = p[j], rho[j] + hb, mb = h[j], m[j] + + ub, vb, wb = u[j], v[j], w[j] + Vb = base.Point(ub,vb,wb) + + havg = 0.5 * (hb + ha) + + rb = base.Point(xb, yb, zb) + + tmp = Pa/(rhoa*rhoa) + Pb/(rhob*rhob) + kernel.py_gradient(ra, rb, havg, grad) + + vab = Va-Vb + rab = ra-rb + + dot = vab.dot(rab) + piab = 0.0 + + if dot < 0.0: + alpha = 1.0 + beta = 1.0 + gamma = 1.4 + eta = 0.1 + + cab = 0.5 * (cs[i] + cs[j]) + + rhoab = 0.5 * (rhoa + rhob) + muab = havg * dot + + muab /= ( rab.norm() + eta*eta*havg*havg ) + + piab = -alpha*cab*muab + beta*muab*muab + piab /= rhoab + + tmp += piab + tmp *= 0.5*mb + + force.x += tmp * ( vab.dot(grad) ) + + forces.append(force) + + return forces + + def test_eval(self): + """ Test the PySPH solution """ + + pa = self.pa + func = self.ewv + + k = base.CubicSplineKernel(dim=2) + + tmpx = pa.properties['tmpx'] + tmpy = pa.properties['tmpy'] + tmpz = pa.properties['tmpz'] + + func.eval(k, tmpx, tmpy, tmpz) + + reference_solution = self.get_reference_solution() + + for i in range(self.np): + self.assertAlmostEqual(reference_solution[i].x, tmpx[i], 6) + self.assertAlmostEqual(reference_solution[i].y, tmpy[i], 6) + self.assertAlmostEqual(reference_solution[i].z, tmpz[i], 6) + + def test_cl_eval(self): + """ Test the PyOpenCL implementation """ + + if solver.HAS_CL: + + pa = self.pa + func = self.ewv + + k = base.CubicSplineKernel(dim=2) + + func.setup_cl(self.prog, self.ctx) + + func.cl_eval(self.q, self.ctx, k) + + pa.read_from_buffer() + + reference_solution = self.get_reference_solution() + + for i in range(self.np): + self.assertAlmostEqual(reference_solution[i].x, pa.tmpx[i], 6) + self.assertAlmostEqual(reference_solution[i].y, pa.tmpy[i], 6) + self.assertAlmostEqual(reference_solution[i].z, pa.tmpz[i], 6) + + +if __name__ == '__main__': + unittest.main() diff --git a/source/pysph/sph/tests/test_pressure_funcs.py b/source/pysph/sph/tests/test_pressure_funcs.py index 2ed9657..fe6842b 100644 --- a/source/pysph/sph/tests/test_pressure_funcs.py +++ b/source/pysph/sph/tests/test_pressure_funcs.py @@ -38,15 +38,20 @@ def setUp(self): x = numpy.array([0, 0, 1, 1], numpy.float64) y = numpy.array([0, 1, 1, 0], numpy.float64) + z = numpy.zeros_like(x) m = numpy.ones_like(x) + + u = numpy.array([1, 0, 0, -1], numpy.float64) + p = numpy.array([0, 0, 1, 1], numpy.float64) + tmpx = numpy.zeros_like(x) tmpy = numpy.zeros_like(x) tmpz = numpy.zeros_like(x) self.pa = pa = base.get_particle_array(name="test", x=x, y=y, z=z, - m=m, tmpx=tmpx, tmpy=tmpy, - tmpz=tmpz, + m=m, u=u, p=p, + tmpx=tmpx, tmpy=tmpy, tmpz=tmpz, cl_precision=self.precision) grad_func = sph.SPHPressureGradient.withargs() @@ -221,7 +226,7 @@ def get_reference_solution(self): ui, vi, wi = u[i], v[i], w[i] ri = base.Point(xi,yi,zi) - va = base.Point(ui,vi,wi) + Va = base.Point(ui,vi,wi) Pi, rhoi = p[i], rho[i] hi = h[i] @@ -231,10 +236,10 @@ def get_reference_solution(self): grad = base.Point() xj, yj, zj = x[j], y[j], z[j] Pj, rhoj = p[j], rho[j] - hj, mj = m[j], h[j] + hj, mj = h[j], m[j] uj, vj, wj = u[j], v[j], w[j] - vb = base.Point(uj,vj,wj) + Vb = base.Point(uj,vj,wj) havg = 0.5 * (hi + hj) @@ -243,7 +248,7 @@ def get_reference_solution(self): tmp = Pi/(rhoi*rhoi) + Pj/(rhoj*rhoj) kernel.py_gradient(ri, rj, havg, grad) - vab = va-vb + vab = Va-Vb rab = ri-rj dot = vab.dot(rab) @@ -260,9 +265,9 @@ def get_reference_solution(self): rhoab = 0.5 * (rhoi + rhoj) muab = havg * dot - mu /= ( rab.length() + eta*eta*havg*havg ) + muab /= ( rab.norm() + eta*eta*havg*havg ) - piab -alpha*cab*mu + beta*mu*mu + piab = -alpha*cab*muab + beta*muab*muab piab /= rhoab tmp += piab