1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35 Andreas Aigner (JKU Linz))
36
37 Copyright 2009-2012 JKU Linz
38 ------------------------------------------------------------------------- */
39
40 #include <cmath>
41 #include <mpi.h>
42 #include <string.h>
43 #include <stdlib.h>
44 #include "fix_sph.h"
45 #include "update.h"
46 #include "respa.h"
47 #include "atom.h"
48 #include "atom_vec.h"
49 #include "force.h"
50 #include "pair_sph.h"
51 #include "comm.h"
52 #include "neighbor.h"
53 #include "neigh_list.h"
54 #include "neigh_request.h"
55 #include "memory.h"
56 #include "error.h"
57 #include "sph_kernels.h"
58 #include "modify.h"
59 #include "fix_property_atom.h"
60 #include "fix_property_global.h"
61
62 using namespace LAMMPS_NS;
63 using namespace FixConst;
64
65 /* ---------------------------------------------------------------------- */
66
FixSph(LAMMPS * lmp,int narg,char ** arg)67 FixSph::FixSph(LAMMPS *lmp, int narg, char **arg) :
68 Fix(lmp, narg, arg)
69 {
70 kernel_flag = 1; // default: kernel is used
71 kernel_id = -1; // default value
72 kernel_cut = -1;
73 kernel_style = NULL;
74
75 fppaSl = NULL;
76 fppaSlType = NULL;
77 sl = NULL;
78 slComType = NULL;
79 }
80
81 /* ---------------------------------------------------------------------- */
82
~FixSph()83 FixSph::~FixSph()
84 {
85 if(kernel_style) delete []kernel_style;
86 if(mass_type) memory->destroy(slComType);
87 }
88
89 /* ---------------------------------------------------------------------- */
90
setmask()91 int FixSph::setmask()
92 {
93 int mask = 0;
94 mask |= POST_INTEGRATE;
95 mask |= POST_INTEGRATE_RESPA;
96 return mask;
97 }
98
99 /* ---------------------------------------------------------------------- */
100
init()101 void FixSph::init()
102 {
103 mass_type = atom->avec->mass_type;
104 int ntypes = atom->ntypes;
105 // need a half neighbor list, built when ever re-neighboring occurs
106
107 int irequest = neighbor->request((void *) this);
108 neighbor->requests[irequest]->pair = 0;
109 neighbor->requests[irequest]->fix = 1;
110
111 if (strcmp(update->integrate_style,"respa") == 0)
112 nlevels_respa = ((Respa *) update->integrate)->nlevels;
113
114 // check if kernel id is set
115 if (kernel_flag && kernel_id < 0) error->all(FLERR,"No sph kernel for fixes is set.");
116 // set kernel_cut
117 kernel_cut = SPH_KERNEL_NS::sph_kernel_cut(kernel_id);
118
119 // get the fix_property containing the smoothing length
120 if (mass_type) {
121 if (fppaSlType == NULL) {
122 fppaSlType=static_cast<FixPropertyGlobal*>(modify->find_fix_property("sl","property/global","peratomtype",ntypes,0,force->pair_style));
123 }
124 if (!fppaSlType) error->all(FLERR,"Fix sph only works with a fix property/global that defines sl");
125
126 // allocate memory for per atom-type property
127 //TODO: copy slComType from pair?
128 if (!slComType) memory->create(slComType,ntypes+1,ntypes+1,"fix:slComType");
129
130 for (int i = 1; i <= ntypes; i++)
131 for (int j = i; j <= ntypes; j++) {
132 double sli = fppaSlType->compute_vector(i-1);
133 double slj = fppaSlType->compute_vector(j-1);
134
135 slComType[i][j] = slComType[j][i] = interpDist(sli,slj);;
136 }
137
138 } else {
139 if (fppaSl == NULL) {
140 fppaSl=static_cast<FixPropertyAtom*>(modify->find_fix_property("sl","property/atom","scalar",0,0,"FixSph",false));
141 }
142 if(!fppaSl) error->all(FLERR,"Fix sph only works with a fix property/atom that defines sl. Internal error!");
143 }
144 }
145
146 /* ---------------------------------------------------------------------- */
147
init_list(int id,NeighList * ptr)148 void FixSph::init_list(int id, NeighList *ptr)
149 {
150 list = ptr;
151 }
152
153 /* ---------------------------------------------------------------------- */
154
post_integrate_respa(int ilevel,int iloop)155 void FixSph::post_integrate_respa(int ilevel, int iloop)
156 {
157 if (ilevel == nlevels_respa-1) post_integrate();
158 }
159
160 /* ---------------------------------------------------------------------- */
161
updatePtrs()162 void FixSph::updatePtrs()
163 {
164 if (fppaSl) sl = fppaSl->vector_atom;
165 if (fppaSlType) sl = fppaSlType->values;
166 }
167
168 /* ----------------------------------------------------------------------
169 return common radius for types with different smoothing lengths
170 atm: simple arithmetic mean (compare Morris)
171 ------------------------------------------------------------------------- */
172 /*
173 inline double FixSph::interpDist(double disti, double distj)
174 {
175 return 0.5*(disti+distj);
176 }
177 */
178