00001
00002
00003
00004
00005
00006 #ifndef AMROC_PROBLEM_H
00007 #define AMROC_PROBLEM_H
00008
00009 #include <cmath>
00010 #include "euler3.h"
00011 #include "ClpProblem.h"
00012
00013 #define f_combl_write FORTRAN_NAME(cblwrite, CBLWRITE)
00014 #define f_ipvel FORTRAN_NAME(ipvel, IPVEL)
00015
00016 extern "C" {
00017 void f_combl_write(INTEGER& NSph, DOUBLE* xm, DOUBLE* vm, DOUBLE rd[]);
00018 void f_ipvel();
00019 }
00020
00021 #define OWN_GFMAMRSOLVER
00022 #include "ClpStdGFMProblem.h"
00023 #include "AMRGFMInterpolation.h"
00024
00025 class SolverSpecific :
00026 public AMRGFMSolver<VectorType,FixupType,FlagType,DIM> {
00027 typedef VectorType::InternalDataType DataType;
00028 typedef AMRGFMSolver<VectorType,FixupType,FlagType,DIM> base;
00029 typedef AMRGFMInterpolation<VectorType,FixupType,FlagType,DIM> interpolation_type;
00030 typedef interpolation_type::point_type point_type;
00031 typedef F77GFMFileOutput<VectorType,FixupType,FlagType,DIM> output_type;
00032 public:
00033 SolverSpecific(IntegratorSpecific& integ,
00034 base::initial_condition_type& init,
00035 base::boundary_conditions_type& bc) : base(integ, init, bc), CouplingLevel(0),
00036 UpdateForceEvery(1), NSph(0), Np(0), xm(0), vm(0), fm(0),
00037 xm_old(0), vm_old(0), fm_old(0), cm(0), rd(0), tm(0) {
00038 SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
00039 #ifdef f_flgout
00040 SetFileOutput(new F77GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this,f_flgout));
00041 #else
00042 SetFileOutput(new GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this));
00043 #endif
00044 SetFixup(new FixupSpecific(integ));
00045 SetFlagging(new FlaggingSpecific(*this));
00046 AddGFM(new GhostFluidMethod<VectorType,DIM>(
00047 new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans,f_ipvel),
00048 new F77GFMLevelSet<DataType,DIM>(f_lset)));
00049 _Interpolation = new interpolation_type(*this);
00050 pi = 4.0*std::atan(1.0);
00051 }
00052
00053 ~SolverSpecific() {
00054 DeleteGFM(_GFM[0]);
00055 delete _Flagging;
00056 delete _Fixup;
00057 delete _Interpolation;
00058 if (Np) delete [] Np;
00059 if (xm) delete [] xm;
00060 if (vm) delete [] vm;
00061 if (fm) delete [] fm;
00062 if (xm_old) delete [] xm_old;
00063 if (vm_old) delete [] vm_old;
00064 if (fm_old) delete [] fm_old;
00065 if (cm) delete [] cm;
00066 if (rd) delete [] rd;
00067 if (tm) delete [] tm;
00068 }
00069
00070 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00071 base::register_at(Ctrl,prefix);
00072 RegisterAt(base::LocCtrl,"CouplingLevel",CouplingLevel);
00073 RegisterAt(base::LocCtrl,"UpdateForceEvery",UpdateForceEvery);
00074 }
00075 virtual void register_at(ControlDevice& Ctrl) {
00076 base::register_at(Ctrl);
00077 }
00078
00079 virtual void SetupData() {
00080 base::SetupData();
00081 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00082 if (CouplingLevel>base::MaxLev-1)
00083 CouplingLevel=base::MaxLev-1;
00084 }
00085
00086 virtual void Initialize(double& t, double &dt) {
00087 std::ifstream ifs("spheres.dat", std::ios::in);
00088 ifs >> NSph;
00089 Np = new int[NSph];
00090 xm = new point_type[NSph];
00091 vm = new point_type[NSph];
00092 fm = new point_type[NSph];
00093 cm = new DataType[NSph];
00094 rd = new DataType[NSph];
00095 tm = new DataType[NSph];
00096 for (int n=0; n<NSph; n++) {
00097 ifs >> xm[n](0) >> xm[n](1) >> xm[n](2) >> rd[n]
00098 >> vm[n](0) >> vm[n](1) >> vm[n](2) >> cm[n]
00099 >> tm[n] >> Np[n];
00100 cm[n] = cm[n]*pi*4./3.*rd[n]*rd[n]*rd[n];
00101 fm[n] = 0.;
00102 }
00103 f_combl_write(NSph,&xm[0](0),&vm[0](0),rd);
00104 ifs.close();
00105 base::Initialize(t,dt);
00106 OutputFiles(t);
00107 }
00108
00109 virtual void Restart(double& t, double& dt) {
00110 std::ifstream ifs("spheres.cp", std::ios::in);
00111 ifs.read((char*)&NSph,sizeof(int));
00112 Np = new int[NSph];
00113 xm = new point_type[NSph];
00114 vm = new point_type[NSph];
00115 fm = new point_type[NSph];
00116 cm = new DataType[NSph];
00117 rd = new DataType[NSph];
00118 tm = new DataType[NSph];
00119 ifs.read((char*)Np,NSph*sizeof(int));
00120 ifs.read((char*)xm,NSph*sizeof(point_type));
00121 ifs.read((char*)vm,NSph*sizeof(point_type));
00122 ifs.read((char*)fm,NSph*sizeof(point_type));
00123 ifs.read((char*)cm,NSph*sizeof(DataType));
00124 ifs.read((char*)rd,NSph*sizeof(DataType));
00125 ifs.read((char*)tm,NSph*sizeof(DataType));
00126 ifs.close();
00127 f_combl_write(NSph,&xm[0](0),&vm[0](0),rd);
00128 base::Restart(t,dt);
00129 }
00130 virtual void Checkpointing() {
00131 int me = MY_PROC;
00132 if (me == VizServer) {
00133 std::ofstream ofs("spheres.cp", std::ios::out);
00134 ofs.write((char*)&NSph,sizeof(int));
00135 ofs.write((char*)Np,NSph*sizeof(int));
00136 ofs.write((char*)xm,NSph*sizeof(point_type));
00137 ofs.write((char*)vm,NSph*sizeof(point_type));
00138 ofs.write((char*)fm,NSph*sizeof(point_type));
00139 ofs.write((char*)cm,NSph*sizeof(DataType));
00140 ofs.write((char*)rd,NSph*sizeof(DataType));
00141 ofs.write((char*)tm,NSph*sizeof(DataType));
00142 ofs.close();
00143 }
00144 base::Checkpointing();
00145 }
00146
00147 virtual void Restart_(std::stringstream& CheckpointStr) {
00148 base::Restart_(CheckpointStr);
00149 for (register int n=0; n<NSph; n++) {
00150 xm[n] = xm_old[n]; vm[n] = vm_old[n]; fm[n] = fm_old[n];
00151 }
00152 }
00153 virtual void Checkpointing_(std::stringstream& CheckpointStr) {
00154 base::Checkpointing_(CheckpointStr);
00155 if (!xm_old) xm_old = new point_type[NSph];
00156 if (!vm_old) vm_old = new point_type[NSph];
00157 if (!fm_old) fm_old = new point_type[NSph];
00158 for (register int n=0; n<NSph; n++) {
00159 xm_old[n] = xm[n]; vm_old[n] = vm[n]; fm_old[n] = fm[n];
00160 }
00161 }
00162
00163 virtual double IntegrateLevel(vec_grid_fct_type& u, const int Time, const int Level,
00164 double t, double dt, bool DoFixup, double tc,
00165 const int which) {
00166
00167 double cfl = base::IntegrateLevel(u,Time,Level,t,dt,DoFixup,tc,which);
00168 #ifdef DEBUG_PRINT
00169 ( comm_service::log() << "*** base::IntegrateLevel() " << Level << " finished.\n" ).flush();
00170 #endif
00171
00172 if (which==0 && Level==CouplingLevel) {
00173 point_type* fmold = new point_type[NSph];
00174 for (int ns=0; ns<NSph; ns++)
00175 fmold[ns] = fm[ns];
00176 point_type a;
00177 if (Time % UpdateForceEvery == 0)
00178 CalculateForce();
00179 for (int ns=0; ns<NSph; ns++) {
00180 if (t+dt<=tm[ns]) continue;
00181 a = 0.5/cm[ns]*(fm[ns]+fmold[ns]);
00182 vm[ns] += a*dt;
00183 xm[ns] += vm[ns]*dt;
00184 }
00185 delete [] fmold;
00186 f_combl_write(NSph,&xm[0](0),&vm[0](0),rd);
00187 OutputFiles(t+dt);
00188 }
00189
00190 return cfl;
00191 }
00192
00193 void CalculateForce() {
00194
00195 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00196 int Time = CurrentTime(base::GH(),l);
00197 int press_cnt = base::Dim()+4;
00198 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00199 press_cnt, base::t[l]);
00200 }
00201
00202 double *force = new double[3*NSph];
00203 for (int ns=0; ns<NSph; ns++) {
00204 int Npoints = 8*Np[ns]*Np[ns];
00205 DataType* p = new DataType[Npoints];
00206 point_type* xc = new point_type[Npoints];
00207 point_type* normals = new point_type[Npoints];
00208 DataType* ds = new DataType[Npoints];
00209 DataType* dist = (DataType*) 0;
00210 double dalpha = 0.5*pi/Np[ns];
00211 double dbeta = 0.5*pi/Np[ns];
00212 register int n = 0;
00213
00214
00215 for (register int l=0; l<4*Np[ns]; l++) {
00216 double alpha = dalpha*(l+0.5);
00217 for (register int m=0; m<2*Np[ns]; m++) {
00218 double beta = dalpha*(m+0.5);
00219 xc[n](0) = xm[ns](0)+rd[ns]*std::cos(alpha)*std::sin(beta);
00220 xc[n](1) = xm[ns](1)+rd[ns]*std::sin(alpha)*std::sin(beta);
00221 xc[n](2) = xm[ns](2)+rd[ns]*std::cos(beta);
00222 ds[n] = pi*rd[ns]*rd[ns]*dbeta*std::sin(beta)/(2.*Np[ns]);
00223 n++;
00224 }
00225 }
00226
00227
00228 _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
00229 _Interpolation->PointsGeom(base::GFM(0).Phi(),base::t[0],Npoints,xc,normals,dist);
00230
00231 force[3*ns ] = 0.0;
00232 force[3*ns+1] = 0.0;
00233 force[3*ns+2] = 0.0;
00234 for (n=0; n<Npoints; n++)
00235 if (p[n]>-1.e37) {
00236 force[3*ns ] += p[n]*normals[n](0)*ds[n];
00237 force[3*ns+1] += p[n]*normals[n](1)*ds[n];
00238 force[3*ns+2] += p[n]*normals[n](2)*ds[n];
00239 }
00240
00241 delete [] p;
00242 delete [] xc;
00243 delete [] normals;
00244 delete [] ds;
00245 }
00246
00247 #ifdef DAGH_NO_MPI
00248 for (int ns=0; ns<NSph; ns++) {
00249 fm[ns](0) = force[3*ns ];
00250 fm[ns](1) = force[3*ns+1];
00251 fm[ns](2) = force[3*ns+2];
00252 }
00253 #else
00254 int Nval = 3*NSph;
00255 double *global_force = new double[3*NSph];
00256 MPI_Allreduce(force, global_force, Nval, MPI_DOUBLE, MPI_SUM, base::Comm());
00257 for (int ns=0; ns<NSph; ns++) {
00258 fm[ns](0) = global_force[3*ns ];
00259 fm[ns](1) = global_force[3*ns+1];
00260 fm[ns](2) = global_force[3*ns+2];
00261 }
00262 delete [] global_force;
00263 #endif
00264
00265 delete [] force;
00266 }
00267
00268 void OutputFiles(double t) {
00269 int me = MY_PROC;
00270 if (me == VizServer) {
00271 std::ofstream ofs1("positions.txt", std::ios::out | std::ios::app );
00272 ofs1 << t;
00273 for (int ns=0; ns<NSph; ns++)
00274 ofs1 << " " << xm[ns](0) << " " << xm[ns](1) << " " << xm[ns](2);
00275 ofs1 << std::endl;
00276 ofs1.close();
00277
00278 std::ofstream ofs2("forces.txt", std::ios::out | std::ios::app );
00279 ofs2 << t;
00280 for (int ns=0; ns<NSph; ns++)
00281 ofs2 << " " << fm[ns](0) << " " << fm[ns](1) << " " << fm[ns](2);
00282 ofs2 << std::endl;
00283 ofs2.close();
00284
00285 std::ofstream ofs3("velocities.txt", std::ios::out | std::ios::app );
00286 ofs3 << t;
00287 for (int ns=0; ns<NSph; ns++)
00288 ofs3 << " " << vm[ns](0) << " " << vm[ns](1) << " " << vm[ns](2);
00289 ofs3 << std::endl;
00290 ofs3.close();
00291 }
00292 }
00293
00294 protected:
00295 interpolation_type* _Interpolation;
00296 int CouplingLevel, UpdateForceEvery, NSph, *Np;
00297 point_type *xm, *vm, *fm, *xm_old, *vm_old, *fm_old;
00298 DataType *cm, *rd, *tm, pi;
00299 };