From f5f993fcc32de90f676150e60235071c5abcb120 Mon Sep 17 00:00:00 2001 From: "kunal.r.puri@gmail.com" Date: Thu, 5 May 2011 10:22:05 +0530 Subject: [PATCH] Changing ADKEConduction coefficient function to handle multiple arrays --- examples/shock-tube/sigalotti_shock_tube.py | 11 ++++- source/pysph/sph/funcs/adke_funcs.pxd | 3 +- source/pysph/sph/funcs/adke_funcs.pyx | 51 ++++++++------------- 3 files changed, 31 insertions(+), 34 deletions(-) diff --git a/examples/shock-tube/sigalotti_shock_tube.py b/examples/shock-tube/sigalotti_shock_tube.py index 8d524d2..cc91c49 100644 --- a/examples/shock-tube/sigalotti_shock_tube.py +++ b/examples/shock-tube/sigalotti_shock_tube.py @@ -31,6 +31,7 @@ pa = particles.get_named_particle_array('fluid') pa.add_property({'name':'rhop'}) +pa.add_property({'name':'div'}) # ensure that the array 'q' is available @@ -57,12 +58,19 @@ # add the update conduction coefficient after the density calculation +s.add_operation(solver.SPHOperation( + + sph.VelocityDivergence.withargs(), + on_types=[Fluid], from_types=[Fluid], updates=['div'], id='vdivergence'), + + before=False, id='density') + s.add_operation(solver.SPHOperation( sph.ADKEConductionCoeffUpdate.withargs(g1=g1, g2=g2), on_types=[Fluid], from_types=[Fluid], updates=['q'], id='qcoeff'), - before=False, id='density') + before=False, id='vdivergence') # add the artificial heat term after the energy equation @@ -74,7 +82,6 @@ before=False, id="enr") - s.set_final_time(0.15) s.set_time_step(3e-4) diff --git a/source/pysph/sph/funcs/adke_funcs.pxd b/source/pysph/sph/funcs/adke_funcs.pxd index 8cb353a..805cf8a 100644 --- a/source/pysph/sph/funcs/adke_funcs.pxd +++ b/source/pysph/sph/funcs/adke_funcs.pxd @@ -37,6 +37,7 @@ cdef class SPHVelocityDivergence(SPHFunctionParticle): ############################################################################### # `ADKEConductionCoeffUpdate` class. ############################################################################### -cdef class ADKEConductionCoeffUpdate(SPHVelocityDivergence): +cdef class ADKEConductionCoeffUpdate(SPHFunction): """ Compute the new smoothing length for the ADKE algorithm """ cdef double g1, g2 + cdef DoubleArray d_div diff --git a/source/pysph/sph/funcs/adke_funcs.pyx b/source/pysph/sph/funcs/adke_funcs.pyx index 4479b22..c6d4840 100644 --- a/source/pysph/sph/funcs/adke_funcs.pyx +++ b/source/pysph/sph/funcs/adke_funcs.pyx @@ -98,7 +98,7 @@ cdef class ADKESmoothingUpdate(SPHFunction): cpdef setup_arrays(self): """ """ - SPHFunctionParticle.setup_arrays(self) + SPHFunction.setup_arrays(self) self.d_rhop = self.dest.get_carray('rhop') @@ -220,13 +220,13 @@ cdef class SPHVelocityDivergence(SPHFunctionParticle): ############################################################################### # `ADKEConductionCoeffUpdate` class. ############################################################################### -cdef class ADKEConductionCoeffUpdate(SPHVelocityDivergence): +cdef class ADKEConductionCoeffUpdate(SPHFunction): """ Compute the pilot estimate of density for the ADKE algorithm """ def __init__(self, ParticleArray source, ParticleArray dest, bint setup_arrays=True, g1=0.0, g2=0.0, **kwargs): - SPHVelocityDivergence.__init__(self,source,dest,setup_arrays,**kwargs) + SPHFunction.__init__(self,source,dest,setup_arrays,**kwargs) self.g1 = g1 self.g2 = g2 @@ -236,8 +236,6 @@ cdef class ADKEConductionCoeffUpdate(SPHVelocityDivergence): self.cl_kernel_src_file = "adke_funcs.cl" - self.dst_reads.append('cs') - self.cl_kernel_function_name = "ADKEConductionCoeffUpdate" def set_src_dst_reads(self): @@ -246,6 +244,13 @@ cdef class ADKEConductionCoeffUpdate(SPHVelocityDivergence): self.src_reads.extend( ['x','y','z','h','u','v','w','m'] ) self.dst_reads.extend( ['x','y','z','h','u','v','w','rho','cs','tag'] ) + + cpdef setup_arrays(self): + """ + """ + SPHFunction.setup_arrays(self) + + self.d_div = self.dest.get_carray('div') cpdef eval(self, KernelBase kernel, DoubleArray output1, DoubleArray output2, DoubleArray output3): @@ -258,42 +263,26 @@ cdef class ADKEConductionCoeffUpdate(SPHVelocityDivergence): """ - cdef double div, g1, g2, ca, ha + cdef double g1, g2, ca, ha cdef int i self.setup_iter_data() cdef size_t np = self.dest.num_real_particles cdef LongArray tag_arr = self.dest.get_carray('tag') + cdef DoubleArray div = self.d_div g1 = self.g1 g2 = self.g2 for i in range(np): - self.eval_single(i, kernel, &div) - - ca = self.d_cs.data[i] - ha = self.d_h.data[i] - - abs_div = fabs(div) - - # set q_a = g1 h_a c_a + g2 h_a^2 [abs(div_a) - div_a] - - output1.data[i] = g1 * ca + ( g2 * ha * (abs_div - div) ) - output1.data[i] *= ha - - cdef void eval_single(self, size_t dest_pid, KernelBase kernel, - double * result): - """ Computes contribution of all neighbors on particle at dest_pid """ - - cdef LongArray nbrs = self.nbr_locator.get_nearest_particles(dest_pid) - cdef size_t nnbrs = nbrs.length + if tag_arr.data[i] == LocalReal: + ca = self.d_cs.data[i] + ha = self.d_h.data[i] - if self.exclude_self: - if self.src is self.dest: - nnbrs -= 1 + abs_div = fabs( div.data[i] ) - result[0] = 0.0 - for j in range(nnbrs): - self.eval_nbr(nbrs.data[j], dest_pid, kernel, result) - + # set q_a = g1 h_a c_a + g2 h_a^2 [abs(div_a) - div_a] + + output1.data[i] = g1*ca + ( g2 * ha * (abs_div - div.data[i]) ) + output1.data[i] *= ha