diff --git a/examples/shock-tube/sigalotti_shock_tube.py b/examples/shock-tube/sigalotti_shock_tube.py index 3506a18..cc91c49 100644 --- a/examples/shock-tube/sigalotti_shock_tube.py +++ b/examples/shock-tube/sigalotti_shock_tube.py @@ -30,6 +30,8 @@ name='fluid', type=Fluid) pa = particles.get_named_particle_array('fluid') +pa.add_property({'name':'rhop'}) +pa.add_property({'name':'div'}) # ensure that the array 'q' is available @@ -42,19 +44,33 @@ s.add_operation(solver.SPHOperation( - sph.ADKESmoothingUpdate.withargs(h0=h0, k=k, eps=eps), - on_types=[Fluid], from_types=[Fluid], updates=['h'], id='adke'), + sph.ADKEPilotRho.withargs(h0=h0), + on_types=[Fluid], from_types=[Fluid], updates=['rhop'], id='adke_rho'), + + before=True, id="density") - before=True, id="density") +s.add_operation(solver.SPHOperation( + sph.ADKESmoothingUpdate.withargs(h0=h0, k=k, eps=eps), + on_types=[Fluid], updates=['h'], id='adke'), + + before=True, id="density") + # 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 @@ -66,7 +82,6 @@ before=False, id="enr") - s.set_final_time(0.15) s.set_time_step(3e-4) diff --git a/source/pysph/base/kernels.pyx b/source/pysph/base/kernels.pyx index f6290b3..221fdbe 100755 --- a/source/pysph/base/kernels.pyx +++ b/source/pysph/base/kernels.pyx @@ -337,7 +337,7 @@ cdef class DummyKernel(KernelBase): return 1.0 cpdef double radius(self): - return 1.0 + return 2.0 ############################################################################## #`CubicSplineKernel` diff --git a/source/pysph/solver/integrator.py b/source/pysph/solver/integrator.py index 6d0f140..51dd61f 100644 --- a/source/pysph/solver/integrator.py +++ b/source/pysph/solver/integrator.py @@ -310,6 +310,8 @@ def setup_integrator(self): self.initial_props[calc.id].append(prop_initial) + # set the step properties for non integrating calcs + for l in range(self.nsteps): k_num = 'k'+str(l+1) @@ -442,66 +444,31 @@ def eval(self, calcs): for i in range(ncalcs): calc = calcs[i] - updates = calc.updates - nupdates = calc.nupdates - - # get the destination particle array for this calc - - pa = self.arrays[calc.dnum] - if logger.level < 30: logger.info("Integrator:eval: operating on calc %d, %s"%( i, calc.id)) # Evaluate the calc + if calc.integrates: + calc.sph(*calc.dst_writes[k_num]) - calc.sph(*self.step_props) - - for j in range(nupdates): - update_prop = updates[j] - step_prop = self.step_props[j] - - step_array = pa.get(step_prop) - - if not calc.integrates: - - #set the evaluated property - if logger.level < 30: - logger.info("""Integrator:eval: setting the prop - %s for calc %d, %s""" - %(update_prop, i, calc.id)) - - pa.set(**{update_prop:step_array.copy()}) - - # ensure that all processes have reached this point - - particles.barrier() - - # update neighbor information if 'h' has been updated + else: + calc.sph(*calc.updates) - if calc.tag == "h": - particles.update() + # ensure all processes have reached this point + particles.barrier() - # update the remote particle properties + # update the remote particle properties - self.rupdate_list[calc.dnum] = [update_prop] + self.rupdate_list[calc.dnum] = [calc.updates] - if logger.level < 30: - logger.info("""Integrator:eval: updating remote particle - properties %s"""%(self.rupdate_list)) + if logger.level < 30: + logger.info("""Integrator:eval: updating remote particle + properties %s"""%(self.rupdate_list)) - particles.update_remote_particle_properties( - self.rupdate_list) - - else: - k_prop = self.k_props[calc.id][k_num][j] - - pa.set(**{k_prop:step_array.copy()}) - - pass - - #ensure that the eval phase is completed for all processes - + particles.update_remote_particle_properties( + self.rupdate_list) + particles.barrier() def step(self, calcs, dt): diff --git a/source/pysph/sph/funcs/adke_funcs.pxd b/source/pysph/sph/funcs/adke_funcs.pxd index ae7b8ec..805cf8a 100644 --- a/source/pysph/sph/funcs/adke_funcs.pxd +++ b/source/pysph/sph/funcs/adke_funcs.pxd @@ -1,7 +1,8 @@ """ Declarations for the adke functions """ #sph imports -from pysph.sph.sph_func cimport SPHFunctionParticle, CSPHFunctionParticle +from pysph.sph.sph_func cimport SPHFunctionParticle, CSPHFunctionParticle,\ + SPHFunction #base imports from pysph.base.particle_array cimport ParticleArray @@ -20,9 +21,10 @@ cdef class ADKEPilotRho(CSPHFunctionParticle): ############################################################################### # `ADKESmoothingUpdate` class. ############################################################################### -cdef class ADKESmoothingUpdate(ADKEPilotRho): +cdef class ADKESmoothingUpdate(SPHFunction): """ Compute the new smoothing length for the ADKE algorithm """ - cdef double k, eps + cdef double k, eps, h0 + cdef DoubleArray d_rhop ############################################################################### @@ -35,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 400646e..c6d4840 100644 --- a/source/pysph/sph/funcs/adke_funcs.pyx +++ b/source/pysph/sph/funcs/adke_funcs.pyx @@ -4,6 +4,7 @@ from libc.math cimport log, exp from pysph.base.point cimport cPoint, cPoint_new, cPoint_sub, cPoint_dot +from pysph.base.particle_array cimport LocalReal cdef extern from "math.h": double fabs (double) @@ -27,6 +28,14 @@ cdef class ADKEPilotRho(CSPHFunctionParticle): self.cl_kernel_src_file = "adke_funcs.cl" self.cl_kernel_function_name = "ADKEPilotRho" + def set_src_dst_reads(self): + + self.src_reads = [] + self.dst_reads = [] + + self.src_reads.extend( ['x','y','z','m','rho'] ) + self.dst_reads.extend( ['x','y','z','tag'] ) + cdef void eval_nbr_csph(self, size_t source_pid, size_t dest_pid, KernelBase kernel, double *nr, double* dnr): """ Compute the contribution from source_pid on dest_pid. @@ -67,78 +76,66 @@ cdef class ADKEPilotRho(CSPHFunctionParticle): ############################################################################### # `ADKESmoothingUpdate` class. ############################################################################### -cdef class ADKESmoothingUpdate(ADKEPilotRho): +cdef class ADKESmoothingUpdate(SPHFunction): """ Compute the pilot estimate of density for the ADKE algorithm """ def __init__(self, ParticleArray source, ParticleArray dest, bint setup_arrays=True, k=1.0, eps=0.0, double h0=1.0, **kwargs): - ADKEPilotRho.__init__(self, source, dest, - setup_arrays=setup_arrays, h0=h0, **kwargs) + SPHFunction.__init__(self, source, dest, setup_arrays, **kwargs) + self.k = k self.eps = eps + self.h0 = h0 - self.id = "adke_smoothing" + self.id = "adke_smoothing_update" self.tag = "h" self.cl_kernel_src_file = "adke_funcs.cl" self.cl_kernel_function_name = "ADKESmoothingUpdate" + cpdef setup_arrays(self): + """ + """ + SPHFunction.setup_arrays(self) + + self.d_rhop = self.dest.get_carray('rhop') + def set_src_dst_reads(self): self.src_reads = [] - self.dst_reads = [] - - self.src_reads.extend( ['x','y','z','m','rho'] ) - self.dst_reads.extend( ['x','y','z','tag'] ) + self.dst_reads = ['rhop'] cpdef eval(self, KernelBase kernel, DoubleArray output1, DoubleArray output2, DoubleArray output3): """ Evaluate the store the results in the output arrays """ - cdef double result, g, log_g=0.0 + + cdef LongArray tag_arr = self.dest.get_carray('tag') + cdef DoubleArray rhop = self.d_rhop + + cdef double result, g, log_g 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 size_t np = self.dest.get_number_of_particles() + cdef size_t nr = self.dest.num_real_particles + log_g = 0.0 for i in range(np): - self.eval_single(i, kernel, &result) - output1.data[i] = result - log_g += log(result) + if tag_arr.data[i] == LocalReal: + log_g += log( rhop.data[i] ) - log_g /= np + log_g /= nr g = exp(log_g) for i in range(np): - output1.data[i] = self.h0 * self.k * (g/output1.data[i])**self.eps + if tag_arr.data[i] == LocalReal: + output1.data[i] = self.h0 * self.k * (g/rhop.data[i])**self.eps # set the destination's dirty bit since new neighbors are needed - self.dest.set_dirty(True) - - 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 - - cdef double dnr = 0.0 - if self.exclude_self: - if self.src is self.dest: - # this works because nbrs has self particle in last position - nnbrs -= 1 - - result[0] = 0.0 - for j in range(nnbrs): - self.eval_nbr_csph(nbrs.data[j], dest_pid, kernel, result, &dnr) - - if dnr != 0.0: - result[0] /= dnr - + self.dest.set_dirty(True) ############################################################################### # `SPHDivergence` class. @@ -223,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 @@ -239,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): @@ -249,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): @@ -261,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 + if tag_arr.data[i] == LocalReal: + ca = self.d_cs.data[i] + ha = self.d_h.data[i] - 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 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 diff --git a/source/pysph/sph/funcs/energy_funcs.pyx b/source/pysph/sph/funcs/energy_funcs.pyx index f89d695..845f912 100644 --- a/source/pysph/sph/funcs/energy_funcs.pyx +++ b/source/pysph/sph/funcs/energy_funcs.pyx @@ -401,6 +401,10 @@ cdef class ArtificialHeat(SPHFunctionParticle): self.cl_kernel_src_file = "energy_funcs.cl" self.cl_kernel_function_name = "ArtificialHeat" + def set_src_dst_reads(self): + self.src_reads = ['x','y','z','h','m','rho','e,','u','v','w','cs','q'] + self.dst_reads = ['x','y','z','h','rho','e','u','v','w','cs','q'] + cpdef setup_arrays(self): """ Setup the arrays required to read data from source and dest. """ diff --git a/source/pysph/sph/sph_calc.pyx b/source/pysph/sph/sph_calc.pyx index 05f28d3..ff75dc9 100644 --- a/source/pysph/sph/sph_calc.pyx +++ b/source/pysph/sph/sph_calc.pyx @@ -232,6 +232,11 @@ cdef class SPHCalc: self.reset_output_arrays(output1, output2, output3) self.sph_array(output1, output2, output3, exclude_self) + # call an update on the particles if the destination pa is dirty + + if self.dest.is_dirty: + self.particles.update() + cpdef sph_array(self, DoubleArray output1, DoubleArray output2, DoubleArray output3, bint exclude_self=False): """