OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
Sph.cpp
1// ------------------------------------------------------------------------------------------
2// Copyright (C) 2021 ---- absingh
3//
4// This file is part of the Surface DCPSE Paper.
5//
6// This program is free software: you can redistribute it and/or modify
7// it under the terms of the GNU General Public License as published by
8// the Free Software Foundation, either version 3 of the License, or
9// (at your option) any later version.
10//
11// This program is distributed in the hope that it will be useful,
12// but WITHOUT ANY WARRANTY; without even the implied warranty of
13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14// GNU General Public License for more details.
15//
16// You should have received a copy of the GNU General Public License
17// along with this program. If not, see <http://www.gnu.org/licenses/>.
18// ------------------------------------------------------------------------------------------
19
20#include "Vector/vector_dist_subset.hpp"
21#include "DCPSE/DCPSE_op/DCPSE_surface_op.hpp"
22#include "DCPSE/DCPSE_op/DCPSE_Solver.hpp"
23#include "util/SphericalHarmonics.hpp"
24#include <util/PathsAndFiles.hpp>
25
26std::string p_output;
27
28constexpr int DF=0;
29constexpr int F=1;
30constexpr int ANADF=2;
31constexpr int ERR=3;
32constexpr int NORMAL=4;
33constexpr int PCOORD=5;
34
35typedef aggregate<double,double,double,double,Point<3,double>,Point<3,double>> prop;
36openfpm::vector<std::string> propNames{{"Df","f","AnaDf","error","Normal","PCOORD"}};
37typedef vector_dist_ws<3,double,prop> vector_type;
38
39int main(int argc, char * argv[]) {
40 openfpm_init(&argc,&argv);
41 auto & v_cl = create_vcluster();
42 timer tt;
43 tt.start();
44 size_t n;
45 double spL;
46 double grid_spacing_surf;
47 double rCut;
48 constexpr int K = 4;
49 // Get command line arguments
50 if (argc > 3) {
51 std::cout << "Maybe Error: The executable requires the following arguments: number_grid_points order_of_operator" << std::endl;
52 }
53 n = std::stoi(argv[1]);
54 unsigned int ord=std::stoi(argv[2]);
55 //spL=std::stof(argv[2]);
56 //grid_spacing_surf=std::stof(argv[3]);
57 //rCut=std::stof(argv[4]);
58
59 // Files and directories
60 std::string cwd{boost::filesystem::current_path().string()};
61 std::string output_dir{cwd + "/output_Sphere"};
62 create_directory_if_not_exist(output_dir);
63 p_output = output_dir + "/particles";
64 size_t n_sp=n;
65 // Domain
66 double boxP1{-1.5}, boxP2{1.5};
67 double boxSize{boxP2 - boxP1};
68 size_t sz[3] = {n,n,n};
69 double grid_spacing{0.8/((std::pow(sz[0],1/3.0)-1))};;//grid_spacing{boxSize/(sz[0]-1)};
70 grid_spacing_surf=grid_spacing;//grid_spacing*50;
71 rCut=2.9 * grid_spacing_surf; //3.1
72
73 Box<3,double> domain{{boxP1,boxP1,boxP1},{boxP2,boxP2,boxP2}};
74 size_t bc[3] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
75 Ghost<3,double> ghost{rCut + grid_spacing/8.0};
76
77 // particles
78 vector_type particles{0,domain,bc,ghost};
79 particles.setPropNames(propNames);
80
81 Point<3,double> coord;
82 // 1. particles on the Circular surface
83 double theta{0.0};
84 double Golden_angle=M_PI * (3.0 - sqrt(5.0));
85 if (v_cl.rank() == 0) {
86 for(int i=1;i<n_sp-1;i++)
87 {
88 double y = 1.0 - (i /double(n_sp - 1.0)) * 2.0;
89 double radius = sqrt(1 - y * y);
90 double Golden_theta = Golden_angle * i;
91 double x = cos(Golden_theta) * radius;
92 double z = sin(Golden_theta) * radius;
93 particles.add();
94 particles.getLastPos()[0] = x;
95 particles.getLastPos()[1] = y;
96 particles.getLastPos()[2] = z;
97 double rm=sqrt(x*x+y*y+z*z);
98 particles.getLastProp<NORMAL>()[0] = x/rm;
99 particles.getLastProp<NORMAL>()[1] = y/rm;
100 particles.getLastProp<NORMAL>()[2] = z/rm;
101 particles.getLastProp<PCOORD>()[0] = 1.0 ;
102 particles.getLastProp<PCOORD>()[1] = std::atan2(sqrt(x*x+y*y),z);
103 particles.getLastProp<PCOORD>()[2] = std::atan2(y,x);
104 //if((i>n_sp/4.0 && i<=n_sp/4.0+K) || (i>3*n_sp/4.0 && i<=3*n_sp/4.0+K+1))
105 //if(i<20*K+1)
106 double tol=1e-2;
107 if(fabs(particles.getLastProp<PCOORD>()[2])<M_PI+tol && fabs(particles.getLastProp<PCOORD>()[2])>M_PI-tol)
108 {particles.getLastSubset(0);}
109 else if(fabs(particles.getLastProp<PCOORD>()[2])<M_PI/2.+tol && fabs(particles.getLastProp<PCOORD>()[2])>M_PI/2.-tol)
110 {particles.getLastSubset(1);}
111 else if(fabs(particles.getLastProp<PCOORD>()[2])<tol && fabs(particles.getLastProp<PCOORD>()[2])>-tol)
112 {particles.getLastSubset(0);}
113 else
114 {particles.getLastSubset(0);}
115
116
117
118
119 }
120
121 double x,y,z;
122 x=0;
123 y=0;
124 z=-1;
125 particles.add();
126 particles.getLastPos()[0] = x;
127 particles.getLastPos()[1] = y;
128 particles.getLastPos()[2] = z;
129 double rm=sqrt(x*x+y*y+z*z);
130 particles.getLastProp<NORMAL>()[0] = x/rm;
131 particles.getLastProp<NORMAL>()[1] = y/rm;
132 particles.getLastProp<NORMAL>()[2] = z/rm;
133 particles.getLastProp<PCOORD>()[0] = 1.0 ;
134 particles.getLastProp<PCOORD>()[1] = std::atan2(sqrt(x*x+y*y),z);
135 particles.getLastProp<PCOORD>()[2] = std::atan2(y,x);
136 particles.getLastSubset(1);
137 x=0;
138 y=0;
139 z=1;
140 particles.add();
141 particles.getLastPos()[0] = x;
142 particles.getLastPos()[1] = y;
143 particles.getLastPos()[2] = z;
144 rm=sqrt(x*x+y*y+z*z);
145 particles.getLastProp<NORMAL>()[0] = x/rm;
146 particles.getLastProp<NORMAL>()[1] = y/rm;
147 particles.getLastProp<NORMAL>()[2] = z/rm;
148 particles.getLastProp<PCOORD>()[0] = 1.0 ;
149 particles.getLastProp<PCOORD>()[1] = std::atan2(sqrt(x*x+y*y),z);
150 particles.getLastProp<PCOORD>()[2] = std::atan2(y,x);
151 particles.getLastSubset(1);
152 std::cout << "n: " << n << " - grid spacing: " << grid_spacing << " - rCut: " << rCut << "YDtheta" << grid_spacing_surf<<std::endl;
153 }
154
155 particles.map();
156 particles.ghost_get<F>();
157
158 vector_dist_subset<3,double,prop> particles_bulk(particles,0);
159 vector_dist_subset<3,double,prop> particles_boundary(particles,1);
160 auto &bulkIds=particles_bulk.getIds();
161 auto &bdrIds=particles_boundary.getIds();
162 std::unordered_map<const lm,double,key_hash,key_equal> Alm;
163 //Setting max mode l_max
164 //Setting amplitudes to 1
165 for(int l=0;l<=K;l++){
166 for(int m=-l;m<=l;m++){
167 Alm[std::make_tuple(l,m)]=0;
168 }
169 }
170 //2Alm[std::make_tuple(1,0)]=1;
171 Alm[std::make_tuple(K,0)]=1;
172 //Alm[std::make_tuple(3,2)]=1;
173
174 spL=K;
175 auto it2 = particles.getDomainIterator();
176 while (it2.isNext()) {
177 auto p = it2.get();
178 Point<3, double> xP = particles.getProp<PCOORD>(p);
179 /*double Sum=0;
180 for(int m=-spL;m<=spL;++m)
181 {
182 Sum+=openfpm::math::Y(spL,m,xP[1],xP[2]);
183 }*/
184 //particles.getProp<ANADF>(p) = Sum;//openfpm::math::Y(K,K,xP[1],xP[2]);openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);;
185 particles.getProp<ANADF>(p)=openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);
186 particles.getProp<DF>(p)=-(K)*(K+1)*openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Alm);
187 ++it2;
188 }
189 auto DErr=getV<ERR>(particles);
190 auto Af=getV<ANADF>(particles);
191 auto Df=getV<DF>(particles);
192 auto f=getV<F>(particles);
193
194 SurfaceDerivative_xx<NORMAL> Sdxx{particles,ord,rCut,grid_spacing_surf};
195 SurfaceDerivative_yy<NORMAL> Sdyy{particles,ord,rCut,grid_spacing_surf};
196 SurfaceDerivative_zz<NORMAL> Sdzz{particles,ord,rCut,grid_spacing_surf};
197
198 DErr=(Sdxx(Af)+Sdyy(Af)+Sdzz(Af))-f;
199 particles.deleteGhost();
200 particles.write_frame(p_output,0,BINARY);
201
202
203 DCPSE_scheme<equations3d1, decltype(particles)> Solver(particles);
204 auto Eig = Sdxx(f)+Sdyy(f)+Sdzz(f);//+K*(K+1)*f;
205 Solver.impose(Eig, bulkIds, Df);
206 Solver.impose(f, bdrIds, Af);
207 Solver.solve(f);
208
209 // Error
210 double MaxError=0,L2=0;
211 for (int j = 0; j < bulkIds.size(); j++) {
212 auto p = bulkIds.get<0>(j);
213 particles.getProp<ERR>(p) = fabs(particles.getProp<F>(p)-particles.getProp<ANADF>(p));
214 if (particles.getProp<ERR>(p)>MaxError){
215 MaxError=particles.getProp<ERR>(p);
216 }
217 L2+=particles.getProp<ERR>(p)*particles.getProp<ERR>(p);
218 }
219 v_cl.sum(L2);
220 v_cl.max(MaxError);
221 v_cl.execute();
222 L2=sqrt(L2);
223 std::cout.precision(16);
224 if (v_cl.rank()==0){
225 std::cout<<"L2:"<<L2<<std::endl;
226 std::cout<<"L_inf:"<<MaxError<<std::endl;
227 }
228 particles.deleteGhost();
229 particles.write(p_output,BINARY);
230 Sdxx.deallocate(particles);
231 Sdyy.deallocate(particles);
232
233 tt.stop();
234 if (v_cl.rank() == 0)
235 std::cout << "Simulation took: " << tt.getcputime() << " s (CPU) - " << tt.getwct() << " s (wall)\n";
236
237 openfpm_finalize();
238 return 0;
239
240}
PathsAndFiles.hpp
Header file containing functions for creating files and folders.
create_directory_if_not_exist
static void create_directory_if_not_exist(std::string path, bool silent=0)
Creates a directory if not already existent.
Definition PathsAndFiles.hpp:103
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
Point
This class implement the point shape in an N-dimensional space.
Definition Point.hpp:28
openfpm::vector
Implementation of 1-D std::vector like structure.
Definition map_vector.hpp:203
timer
Class for cpu time benchmarking.
Definition timer.hpp:28
timer::stop
void stop()
Stop the timer.
Definition timer.hpp:119
timer::getcputime
double getcputime()
Return the cpu time.
Definition timer.hpp:142
timer::start
void start()
Start the timer.
Definition timer.hpp:90
timer::getwct
double getwct()
Return the elapsed real time.
Definition timer.hpp:130
vector_dist_iterator::get
vect_dist_key_dx get()
Get the actual key.
Definition vector_dist_iterator.hpp:75
vector_dist_subset
Definition vector_dist_subset.hpp:84
vector_dist_ws
Definition vector_dist_subset.hpp:17
vector_dist
Distributed vector.
Definition vector_dist.hpp:305
vector_dist::write_frame
bool write_frame(std::string out, size_t iteration, int opt=VTK_WRITER)
Output particle position and properties.
Definition vector_dist.hpp:2874
vector_dist::getProp
auto getProp(vect_dist_key_dx vec_key) -> decltype(v_prp.template get< id >(vec_key.getKey()))
Get the property of an element.
Definition vector_dist.hpp:826
vector_dist::getLastProp
auto getLastProp() -> decltype(v_prp.template get< id >(0))
Get the property of the last element.
Definition vector_dist.hpp:1153
vector_dist::setPropNames
void setPropNames(const openfpm::vector< std::string > &names)
Set the properties names.
Definition vector_dist.hpp:3191
vector_dist::getDomainIterator
vector_dist_iterator getDomainIterator() const
Get an iterator that traverse the particles in the domain.
Definition vector_dist.hpp:2165
vector_dist::ghost_get
void ghost_get(size_t opt=WITH_POSITION)
It synchronize the properties and position of the ghost particles.
Definition vector_dist.hpp:2481
vector_dist::map
void map(size_t opt=NONE)
It move all the particles that does not belong to the local processor to the respective processor.
Definition vector_dist.hpp:2436
vector_dist::write
bool write(std::string out, int opt=VTK_WRITER)
Output particle position and properties.
Definition vector_dist.hpp:2759
vector_dist::getLastPos
auto getLastPos() -> decltype(v_pos.template get< 0 >(0))
Get the position of the last element.
Definition vector_dist.hpp:1130
vector_dist::add
void add()
Add local particle.
Definition vector_dist.hpp:1094
vector_dist::deleteGhost
void deleteGhost()
Delete the particles on the ghost.
Definition vector_dist.hpp:2820
F
[v_transform metafunction]
Definition variadic_to_vmpl_unit_test.hpp:18
aggregate
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition aggregate.hpp:215