! $Id: particles_potential dhruba.mitra@gmail.com$ ! ! This module takes care of everything related to pairwise interaction ! of particles. It is experimental now (April 2016) ! !** AUTOMATIC CPARAM.INC GENERATION **************************** ! ! CPARAM logical, parameter :: lparticles_potential=.true. ! !*************************************************************** module Particles_potential ! use Cdata use General, only: keep_compiler_quiet use Messages use Particles_cdata use Particles_sub ! use Particles_radius ! implicit none ! include 'particles_potential.h' ! character (len=labellen) :: ppotential='nothing' integer :: sigma_in_grid = 1, cell_in_grid=1 real :: psigma_by_dx=0.1,ppower=19,skin_factor=2.,fampl=1. real :: Rcutoff=0.,dRbin=1.,cell_length=0. real :: pmom_max=6, pmom_min=-2, mom_step=0.25 integer :: Rcutoff_in_grid=1,Nbin_in_Rcutoff=100 real :: rescale_diameter=1. logical :: lpotential integer :: mom_max,mpface,mpedge,mpcorner integer :: npbufl,npbufu ! ! ! ----------- About Potential ---------------------- ! This module calculates all quantities related to particle-particle ! interaction. There may or may not actually be a potential, but ! we shall use this module to calculate diagnostic of particle ! pairs. !--------But if we do have a potential then ..---------- ! Note: psigma below is psigma = psigma_by_dx * dx ; because we always like to think of the ! range of the potential in units of dx. While calculating neighbours we look around for ! sigma_in_grid number of grid points. This should be >= 1 . Also, if ! psigma is larger than dx then sigma_in_grid must be larger too. ! Note : psigma should be much smaller than the dissipation range, maybe even less than a grid-spacing ! as we are not including many terms in the Maxey-Riley equations. As far as the interaction ! with the fluid is concerned our particles are point particles with inertia. But the ! interaction potential gives them an effective radius. The interaction potential is ! typically of the form ! V(r) = function of (r/psigma) , \xi = r/psigma ! The default potential is repulsive ! V(r) = fampl*(1/xi)^(beta) ! with beta = 2*ppowerby2 ! This potential is quite steep (almost hard-sphere) hence the effective force on a particle ! due to other particles which are within a distance of skin_factor*psigma. Particles ! within this distance are included in the neighbourlist. ! integer :: ncell=0 integer :: mcellx=0,mcelly=0,mcellz=0 integer :: arb_factor=10 logical :: lhead_allocated=.false. integer, allocatable, dimension(:,:,:) :: head integer, allocatable,dimension(:) :: link_list real, allocatable,dimension(:) :: mom_array real, allocatable,dimension(:,:,:) :: MomJntPDF,MomColJntPDF integer :: ysteps_int,zsteps_int namelist /particles_potential_init_pars/ & arb_factor,ppotential, cell_in_grid, psigma_by_dx, skin_factor, & sigma_in_grid,fampl,Rcutoff_in_grid,Nbin_in_Rcutoff,rescale_diameter,lpotential ! namelist /particles_potential_run_pars/ & ppotential, cell_in_grid, psigma_by_dx, skin_factor, & sigma_in_grid,fampl,Rcutoff_in_grid,Nbin_in_Rcutoff,rescale_diameter,lpotential ! logical :: idiag_particles_vijm=.false.,idiag_particles_vijrms=.false.,idiag_abs_mom=.false. logical :: idiag_colvel_mom=.false.,idiag_gr=.false. ! contains !*********************************************************************** subroutine register_particles_potential() ! ! Set up indices for access to the fp and dfp arrays ! if (lroot) call svn_id( & "$Id$") endsubroutine register_particles_potential !*********************************************************************** subroutine initialize_particles_potential(fp) ! ! Perform any post-parameter-read initialization i.e. calculate derived ! parameters. ! use particles_radius, only: get_maxrad real, dimension (mpar_loc,mparray), intent (in) :: fp integer :: mom_tmp,imom real :: rmax ! ! assume isotropy ! Rcutoff=Rcutoff_in_grid*dx cell_length=cell_in_grid*dx dRbin=Rcutoff/real(Nbin_in_Rcutoff) ! ! The size of a cell is twice the radius of the biggest particle. This assumes that ! the size of the particles are NOT going to change over time. Otherwise the input ! parameter cell_length, if not equal to zero, sets the size of the cell. ! call get_maxrad(rmax) if (cell_length.eq.0.) then cell_length=4.*rmax endif if (lroot) print*, 'rmax, cell_length',rmax,cell_length ! ! the following line assumes that the domain is roughly size in all three ! directions. If not, we need to code some more ! ncell=int((x(l2)-x(l1))/cell_length)+1 cell_length=(x(l2)-x(l1))/ncell ! ! Assuming uniform distribution we can estimate the number of particles ! in a slab. These number are then multiplied ! an arbitrary factor (arb_factor) for which the default value is 10 ! nslab=arb_factor*(npar/ncpus)/ncell ! ! If we are using many processors then our domain effectively includes ! three (two ?) neighbouring processors in each directions. ! mcellx=ncell;mcelly=ncell;mcellz=ncell ! ! Allocate the arrays head and link_list (only if they have ! not been allocated before) ! if(.not.lhead_allocated) then if (lmpicomm) then lpar_max=max(arb_factor*(npar/ncpus)+6*nslab,mpar_loc) allocate(fp_buffer_in(nslab,mparray)) allocate(fp_buffer_out(nslab,mparray)) else lpar_max=mpar_loc endif allocate(head(-1:mcellx,-1:mcelly,-1:mcellz)) allocate(link_list(lpar_max)) ! ! We also need to allocate a larger array in case of parallel communications ! allocate(fpwn(lpar_max,mparray)) lhead_allocated=.true. fpwn=0. head=0 link_list=0 endif ! ! The following are necessary to calculate the moments ! of the joint PDF on the fly. ! ! mom_max=int((pmom_max-pmom_min)/mom_step)+1 ! allocate(mom_array(mom_max)) ! mom_array=0. ! imom=2 ! mom_tmp=pmom_min ! do while (mom_tmp .le. pmom_max) ! if (mom_tmp .ne. 0) then ! mom_array(imom) = mom_tmp ! imom=imom+1 ! endif ! mom_tmp = mom_tmp+mom_step ! enddo ! write(*,*) mom_array !! ! allocate(MomJntPDF(mom_max,Nbin_in_Rcutoff,ndustrad*(ndustrad+1)/2)) ! MomJntPDF=0. ! allocate(MomColJntPDF(mom_max,Nbin_in_Rcutoff,ndustrad*(ndustrad+1)/2)) ! MomColJntPDF=0. ! endsubroutine initialize_particles_potential !*********************************************************************** subroutine particles_potential_clean_up() ! ! cleanup after the particles_potential module ! if(lhead_allocated) then deallocate(head) deallocate(link_list) deallocate(fpwn) if (lmpicomm) then deallocate(fp_buffer_in) deallocate(fp_buffer_out) endif lhead_allocated=.false. endif deallocate(MomJntPDF) deallocate(MomColJntPDF) ! endsubroutine particles_potential_clean_up !*********************************************************************** subroutine init_particles_potential(f,fp,ineargrid) real, dimension (mx,my,mz,mfarray),intent(in) :: f real, dimension (mpar_loc,mparray),intent(in) :: fp integer, dimension(mpar_loc,3),intent(in) :: ineargrid ! ! initial setting for potential ! call keep_compiler_quiet(f) call keep_compiler_quiet(fp) call keep_compiler_quiet(ineargrid) ! endsubroutine init_particles_potential !*********************************************************************** subroutine construct_link_list(plist_min,plist_max) integer :: ip,imom integer, intent(in) :: plist_min,plist_max integer,dimension(3) :: cell_vec real,dimension(3) :: xxip,my_origin ! if (lroot) print*,'plist_min,plist_max',plist_min,plist_max do ip=plist_min,plist_max xxip=fpwn(ip,ixp:izp) my_origin(1)=x(l1); my_origin(2)=y(m1); my_origin(3)=z(n1) cell_vec=floor((xxip-my_origin)/cell_length) if ( (cell_vec(1) > mcellx) .or. (cell_vec(2) > mcelly) .or. (cell_vec(3) > mcellz)) then print*,'1 mcellx,mcelly,mcellz,cell_vec',mcellx,mcelly,mcellz,cell_vec print*,'my_origin,xxip',my_origin,xxip print*,'cell_length',cell_length endif if ( (cell_vec(1) < -1) .or. (cell_vec(2) < -1 ) .or. (cell_vec(3) < -1) ) then print*,'2 mcellx,mcelly,mcellz,cell_vec',mcellx,mcelly,mcellz,cell_vec print*,'my_origin,xxip',my_origin,xxip print*,'cell_length',cell_length endif link_list(ip)=head(cell_vec(1),cell_vec(2),cell_vec(3)) head(cell_vec(1),cell_vec(2),cell_vec(3))=ip enddo if (lroot) then open(unit=1,file='head.tst',status='unknown') write(1,*) 'head inside construct_link_list',head close(1) endif ! endsubroutine construct_link_list !*********************************************************************** subroutine dxxp_dt_potential(f,df,fp,dfp,ineargrid) ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real, dimension (mpar_loc,mparray) :: fp real, dimension (mpar_loc,mpvar) :: dfp integer, dimension (mpar_loc,3) :: ineargrid ! logical :: lheader, lfirstcall=.true. ! intent (in) :: f, fp, ineargrid intent (inout) :: df, dfp call keep_compiler_quiet(f) call keep_compiler_quiet(df) call keep_compiler_quiet(ineargrid) ! endsubroutine dxxp_dt_potential !*********************************************************************** subroutine dvvp_dt_potential_pencil(f,df,fp,dfp,ineargrid) ! ! ! Jan-2017/dhruba: dummy ! real, dimension (mx,my,mz,mfarray), intent (inout) :: f real, dimension (mx,my,mz,mvar), intent (inout) :: df real, dimension (mpar_loc,mparray), intent (in) :: fp real, dimension (mpar_loc,mpvar), intent (inout) :: dfp integer, dimension (mpar_loc,3), intent (in) :: ineargrid ! call keep_compiler_quiet(f) call keep_compiler_quiet(df) call keep_compiler_quiet(fp) call keep_compiler_quiet(dfp) call keep_compiler_quiet(ineargrid) ! endsubroutine dvvp_dt_potential_pencil !*********************************************************************** subroutine dvvp_dt_potential(f,df,fp,dfp,ineargrid) ! ! Evolution of dust particle velocity. ! ! Jan-2017/dhruba: coded ! use Diagnostics use EquationOfState, only: cs20 ! real, dimension (mx,my,mz,mfarray), intent (inout) :: f real, dimension (mx,my,mz,mvar), intent (inout) :: df real, dimension (mpar_loc,mparray), intent (in) :: fp real, dimension (mpar_loc,mpvar), intent (inout) :: dfp integer, dimension (mpar_loc,3), intent (in) :: ineargrid ! logical :: lheader, lfirstcall=.true. ! ! We need to calculate pairwise quantities only if we are calculating ! pairwise diagnostic or if we actually have a potential of interaction ! if (ldiagnos.or.lparticles_potential) then ! ! first fill up fpwn with particles in local fp (do this even ! in the serial case ) ! lpar_loc=npar_loc fpwn(1:npar_loc,:) = fp(1:npar_loc,:) ! ! Now make the linked list for the local particles ! but set head to zero before you begin head = 0 if (lroot) print*,'going into construct_link_list' call construct_link_list(1,npar_loc) if (lroot) print*,'coming out of construct_link_list' ! ! Now load the data from boundaries to buffers (this is ! different in the serial and parallel case. ! if (lroot) print*,'going into particles_neighbour_proc' call particles_neighbour_proc() if (lroot) print*,'coming out of particles_neighbour_proc' ! if (lroot) print*,'going into calculate_potential' call calculate_potential (dfp,ineargrid) if (lroot) print*,'coming out of calculate_potential' ! endif endsubroutine dvvp_dt_potential !*********************************************************************** subroutine assimilate_incoming(npbuf) integer,intent(in) :: npbuf fpwn(lpar_loc+1:lpar_loc+npbuf,:)=fp_buffer_in(1:npbuf,:) call construct_link_list(lpar_loc+1,lpar_loc+npbuf) lpar_loc=lpar_loc+npbuf ! endsubroutine assimilate_incoming !*********************************************************************** subroutine get_boundary_particles(idirn,porm,npbuf) ! ! fp_buffer in known globally ! integer, intent(in) :: idirn,porm integer, intent(out) :: npbuf ! select case(idirn) case(3) if (porm.eq.1) then call make_fpbuffer(mcellz-1,mcellz-1,-1,mcelly,-1,mcellx,npbuf) else call make_fpbuffer(0,0,-1,mcelly,-1,mcellx,npbuf) endif case(2) if (porm.eq.1) then call make_fpbuffer(-1,mcellz,mcelly-1,mcelly-1,-1,mcellx,npbuf) else call make_fpbuffer(-1,mcellz,0,0,-1,mcellx,npbuf) endif case(1) if (porm.eq.1) then call make_fpbuffer(-1,mcellz,-1,mcelly,mcellx-1,mcellx-1,npbuf) else call make_fpbuffer(-1,mcellz,-1,mcelly,0,0,npbuf) endif case default ! ! Catch unknown values ! call fatal_error("particles_potential", & "get_boundary_particles is called with wrong idirn") ! endselect ! endsubroutine get_boundary_particles !*********************************************************************** subroutine make_fpbuffer(izmin,izmax,iymin,iymax,ixmin,ixmax,npbuf) ! ! fp_buffer is known globally ! integer, intent(in) :: izmin,izmax,iymin,iymax,ixmin,ixmax integer, intent(out) :: npbuf integer :: ipbuf,ix,iy,iz fp_buffer_out=0. npbuf=0 ipbuf=0 do iz=izmin,izmax do iy=iymin,iymax do ix=ixmin,ixmax ip = head(ix,iy,iz) do while (ip.ne.0) ipbuf=ipbuf+1 fp_buffer_out(ipbuf,:)=fpwn(ip,:) ip = link_list(ip) enddo ! loop all the particles in a cell ends enddo enddo enddo npbuf=ipbuf endsubroutine make_fpbuffer !*********************************************************************** subroutine particles_neighbour_proc() ! ! calls the same subroutine 3 times for each direction. ! The sequence of calls is crucial, the inner coding assumes ! this sequence and the results will go wrong if the order ! is changed. ! if (nprocz > 1) then call particles_neighbour_proc_dirn(3) else call apply_pbc_head(3) endif ! if (nprocy > 1) then call particles_neighbour_proc_dirn(2) else call apply_pbc_head(2) endif ! if (nprocx > 1) then call particles_neighbour_proc_dirn(1) else call apply_pbc_head(1) endif endsubroutine particles_neighbour_proc !*********************************************************************** subroutine apply_pbc_head(idirn) integer, intent(in) :: idirn ! select case (idirn) ! copying along first direction ! first the faces case(1) head(-1,-1:mcelly,-1:mcellz) = head(mcellx-1,-1:mcelly,-1:mcellz) head(mcellx,-1:mcelly,-1:mcellz) = head(0,-1:mcelly,-1:mcellz) ! case(2) head(-1:mcellx,-1,-1:mcellz) = head(-1:mcellx,mcelly-1,-1:mcellz) head(-1:mcellx,mcelly,-1:mcellz) = head(-1:mcellx,0,-1:mcellz) ! case(3) head(-1:mcellx,-1:mcelly,-1) = head(-1:mcellx,-1:mcelly,mcellz-1) head(-1:mcellx,-1:mcelly,mcellz) = head(-1:mcellx,-1:mcelly,0) ! endselect end subroutine apply_pbc_head !*********************************************************************** subroutine particles_neighbour_proc_dirn(idirn) use Particles_mpicomm ! use Particles integer, intent(in) :: idirn integer :: uneigh,lneigh integer :: npbuf,my_npbuf,her_npbuf !--------- select case(idirn) case(3) uneigh=zuneigh;lneigh=zlneigh case(2) uneigh=yuneigh;lneigh=ylneigh case(1) uneigh=xuneigh;lneigh=xlneigh case default ! ! Catch unknown values ! call fatal_error("particles_mpicomm", & "wrong value of idirn") ! endselect ! ! buffer for UPPER boundary ! call get_boundary_particles(idirn,1,npbuf) my_npbuf=npbuf ! ! send and receive buffers ! call communicate_fpbuf(uneigh,lneigh,her_npbuf,my_npbuf) ! ! Now my buffer size is changed ! npbuf=her_npbuf call assimilate_incoming(npbuf) ! ! buffer for LOWER boundary ! call get_boundary_particles(idirn,-1,npbuf) my_npbuf=npbuf ! ! send and receive buffers ! call communicate_fpbuf(lneigh,uneigh,her_npbuf,my_npbuf) ! ! Now my buffer size is changed ! npbuf=her_npbuf call assimilate_incoming(npbuf) ! endsubroutine particles_neighbour_proc_dirn !*********************************************************************** subroutine calculate_potential(dfp,ineargrid) ! ! contribution particle acceleration due to particle-particle interaction ! ! use Diagnostics use Sub, only: periodic_fold_back ! real, dimension (mpar_loc,mpvar) :: dfp integer, dimension (mpar_loc,3) :: ineargrid ! logical :: lheader, lfirstcall=.true. ! intent (in) :: ineargrid intent (inout) :: dfp !---------------------------------- integer :: ix,iy,iz,ip,jp,kp integer,parameter :: Nnab=13 integer,dimension(Nnab,3) :: neighbours integer :: inab,ix2 real,dimension(3) :: xxij,vvij ! ! Print out header information in first time step. ! lheader=lfirstcall .and. lroot if (lheader) then print*,'dvvp_dt: Calculate dvvp_dt_potential' endif ! ! At this stage we know the linked list so we access by particles ! by it. ! ! check if the link list has been allocated, otherwise abort if (.not.lhead_allocated) & call fatal_error('dvvp_dt_potential', 'The linked list is not allocated; ABORTING') ! if (lroot) print*,'mcellx,mcelly,mcellz',mcellx,mcelly,mcellz if (lroot) then open(unit=1,file='head2.tst',status='unknown') write(1,*) 'head before mcell loop',head close(1) endif if (lroot) print*,'calculate_potential 1' do iz=0,mcellz-1;do iy=0,mcelly-1;do ix=0,mcellx-1 ip = head(ix,iy,iz) ! ! within the same cell ! do while (ip.ne.0) jp = link_list(ip) do while (jp.ne.0) xxij= fpwn(jp,ixp:izp)-fpwn(ip,ixp:izp) vvij=fpwn(jp,ivpx:ivpz)-fpwn(ip,ivpx:ivpz) call two_particle_int(dfp,ip,jp,xxij,vvij) jp = link_list(jp) enddo ip = link_list(ip) enddo ! loop all the particles in a cell ends ! ! Now for neighbouring cells ! ip = head(ix,iy,iz) call get_cell_neighbours(ix,iy,iz,neighbours) do while (ip.ne.0) do inab = 1,Nnab ix2 = neighbours(inab,1) iy2 = neighbours(inab,2) iz2 = neighbours(inab,3) jp = head(ix2,iy2,iz2) do while (jp.ne.0) xxij= fpwn(jp,ixp:izp)-fpwn(ip,ixp:izp) call periodic_fold_back(xxij, Lxyz) vvij=fpwn(jp,ivpx:ivpz)-fpwn(ip,ivpx:ivpz) call two_particle_int(dfp,ip,jp,xxij,vvij) jp = link_list(jp) enddo enddo enddo ! enddo; enddo; enddo if (lroot) print*,'calculate_potential 2' ! ! loop over cells done ! ! if (lfirstcall) lfirstcall=.false. ! call keep_compiler_quiet(ineargrid) ! endsubroutine calculate_potential !*********************************************************************** subroutine two_particle_int(dfp,ip,jp,xxij,vvij) ! use particles_radius, only: get_stbin use Diagnostics use Sub, only: lower_triangular_index ! real, dimension (mpar_loc,mpvar) :: dfp ! integer,intent(in) :: ip,jp real, dimension(3),intent(in) :: xxij,vvij real,dimension(3) :: accni,accnj real :: Rsqr,Vsqr,Vparallel,pmom,RR,api,apj integer :: iStbin,jStbin,ijSt,iRbin,jmom !!--------------------------------- Rsqr=xxij(1)*xxij(1)+xxij(2)*xxij(2)+xxij(3)*xxij(3) RR = sqrt(Rsqr) if (lpotential) then call get_interparticle_accn(accni,accnj,ip,jp,xxij) dfp(ip,ivpx:ivpz) = dfp(ip,ivpx:ivpz)+accni ! ! The other(jp) particle may be in a different processor. If that is the ! case then we do not add to the dfp ! if (jp .le. npar_loc) & dfp(jp,ivpx:ivpz) = dfp(jp,ivpx:ivpz)+accnj endif ! if (ldiagnos) then ! if (RR .lt. Rcutoff) then ! Vsqr=vvij(1)*vvij(1)+vvij(2)*vvij(2)+vvij(3)*vvij(3) ! Vparallel=dot_product(xxij,vvij)/RR ! iRbin=floor(RR/dRbin) ! api = fpwn(ip,iap) ! call get_stbin(api,iStbin) ! apj = fpwn(jp,iap) ! call get_stbin(apj,jStbin) ! call lower_triangular_index(ijSt,iStbin,jStbin) ! if (idiag_gr) & ! MomJntPDF(1,iRbin,ijSt) = MomJntPDF(1,iRbin,ijSt)+1. ! if (idiag_colvel_mom) then ! if (Vparallel .lt. 0) & ! MomColJntPDF(1,iRbin,ijSt) = MomColJntPDF(1,iRbin,ijSt)+1. ! endif ! do jmom=2,mom_max ! pmom=mom_array(jmom) ! if (idiag_abs_mom) & ! MomJntPDF(jmom,iRbin,ijSt) = MomJntPDF(jmom,iRbin,ijSt)+(abs(Vparallel))**pmom ! if (Vparallel .lt. 0) then ! MomColJntPDF(jmom,iRbin,ijSt) = MomColJntPDF(jmom,iRbin,ijSt)+(-Vparallel)**pmom ! endif ! enddo ! endif ! endif ! endsubroutine two_particle_int !*********************************************************************** subroutine get_cell_neighbours(ix,iy,iz,neighbours) integer,intent(in) :: ix,iy,iz integer,parameter :: Nnab=13 integer,dimension(Nnab,3) :: neighbours ! neighbours(1,1) = ix + 1 neighbours(1,2) = iy neighbours(1,3) = iz ! neighbours(2,1) = ix - 1 neighbours(2,2) = iy - 1 neighbours(2,3) = iz + 1 ! neighbours(3,1) = ix neighbours(3,2) = iy - 1 neighbours(3,3) = iz + 1 ! neighbours(4,1) = ix + 1 neighbours(4,2) = iy - 1 neighbours(4,3) = iz + 1 ! neighbours(5,1) = ix - 1 neighbours(5,2) = iy neighbours(5,3) = iz + 1 ! neighbours(6,1) = ix neighbours(6,2) = iy neighbours(6,3) = iz + 1 ! neighbours(7,1) = ix + 1 neighbours(7,2) = iy neighbours(7,3) = iz + 1 ! neighbours(8,1) = ix - 1 neighbours(8,2) = iy + 1 neighbours(8,3) = iz + 1 ! neighbours(9,1) = ix neighbours(9,2) = iy + 1 neighbours(9,3) = iz + 1 ! neighbours(10,1) = ix + 1 neighbours(10,2) = iy + 1 neighbours(10,3) = iz + 1 ! neighbours(11,1) = ix - 1 neighbours(11,2) = iy + 1 neighbours(11,3) = iz ! neighbours(12,1) = ix neighbours(12,2) = iy + 1 neighbours(12,3) = iz ! neighbours(13,1) = ix + 1 neighbours(13,2) = iy + 1 neighbours(13,3) = iz endsubroutine get_cell_neighbours !*********************************************************************** subroutine get_interparticle_accn(accni,accnj,ip,jp,RR) use particles_radius, only: get_mass_from_radius integer,intent(in) :: ip,jp real, dimension(3),intent(in) :: RR real,dimension(3), intent(out) :: accni,accnj real, dimension (mpar_loc,mparray) :: fp real :: Vij real,dimension(3) :: force_ij real :: mp_i,mp_j ! call get_interaction_force(force_ij,RR,ip,jp) call get_mass_from_radius(mp_i,fpwn,ip) accni=force_ij/mp_i call get_mass_from_radius(mp_j,fpwn,jp) accnj=-force_ij/mp_j ! endsubroutine get_interparticle_accn !*********************************************************************** subroutine get_interaction_force(force_ij,RR,ip,jp) integer, intent(in) :: ip,jp real,dimension(3),intent(in) :: RR real,dimension(3),intent(out) :: force_ij real :: RR_mod real,dimension(3) :: Rcap real :: radiusi,radiusj,diameter_ij,force_amps ! select case (ppotential) case ('rep-power-law-cutoff') ! ! repulsive power law ! RR_mod=sqrt(RR(1)*RR(1)+RR(2)*RR(2)+RR(3)*RR(3)) Rcap=RR/RR_mod radiusi=fpwn(ip,iap) radiusj=fpwn(jp,iap) diameter_ij=rescale_diameter*(radiusi+radiusj) if (RR_mod .gt. diameter_ij) then force_ij=0. else force_amps=fampl*RR_mod**(-ppower) force_ij=-force_amps*Rcap endif case default call fatal_error('particles_potential: no potential coded ','get_interaction_force') endselect ! endsubroutine get_interaction_force !*********************************************************************** subroutine read_particles_pot_init_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=particles_potential_init_pars, IOSTAT=iostat) ! endsubroutine read_particles_pot_init_pars !*********************************************************************** subroutine write_particles_pot_init_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=particles_potential_init_pars) ! endsubroutine write_particles_pot_init_pars !*********************************************************************** subroutine read_particles_pot_run_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=particles_potential_run_pars, IOSTAT=iostat) ! endsubroutine read_particles_pot_run_pars !*********************************************************************** subroutine write_particles_pot_run_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=particles_potential_run_pars) ! endsubroutine write_particles_pot_run_pars !*********************************************************************** subroutine rprint_particles_potential(lreset,lwrite) ! ! Read and register print parameters relevant for particles. ! ! 29-dec-04/anders: coded ! use Diagnostics ! logical :: lreset logical, optional :: lwrite ! integer :: iname,inamez,inamey,inamex,inamexy,inamexz,inamer,inamerz logical :: lwr ! ! Write information to index.pro. ! lwr = .false. if (present(lwrite)) lwr=lwrite ! ! ! Reset everything in case of reset. ! if (lreset) then ! idiag_xpm=0; idiag_ypm=0; idiag_zpm=0 endif ! ! Run through all possible names that may be listed in print.in. ! if (lroot .and. ip<14) print*,'rprint_particles: run through parse list' ! do iname=1,nname ! call parse_name(iname,cname(iname),cform(iname),'nparmin',idiag_nparmin) ! enddo ! endsubroutine rprint_particles_potential !*********************************************************************** subroutine periodic_boundcond_on_aux(f) ! ! Impose periodic boundary condition on bb use Boundcond, only: set_periodic_boundcond_on_aux real, dimension(mx,my,mz,mfarray), intent(in) :: f ! dummy call keep_compiler_quiet(f) endsubroutine periodic_boundcond_on_aux !*********************************************************************** endmodule Particles_potential