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_read FORTRAN_NAME(cblread, CBLREAD)
00014 #define f_combl_write FORTRAN_NAME(cblwrite, CBLWRITE)
00015 #define f_ipvel FORTRAN_NAME(ipvel, IPVEL)
00016
00017 extern "C" {
00018 void f_combl_read(INTEGER& NSph, DOUBLE xm[], DOUBLE vm[], DOUBLE cm[], DOUBLE rd[], INTEGER Np[]);
00019 void f_combl_write(INTEGER& NSph, DOUBLE xm[], DOUBLE vm[], DOUBLE cm[], DOUBLE rd[], INTEGER Np[]);
00020 void f_ipvel();
00021 }
00022
00023 #define OWN_GFMAMRSOLVER
00024 #include "ClpStdGFMProblem.h"
00025 #include "AMRGFMInterpolation.h"
00026
00027 class SolverSpecific :
00028 public AMRGFMSolver<VectorType,FixupType,FlagType,DIM> {
00029 typedef VectorType::InternalDataType DataType;
00030 typedef AMRGFMSolver<VectorType,FixupType,FlagType,DIM> base;
00031 typedef AMRGFMInterpolation<VectorType,FixupType,FlagType,DIM> interpolation_type;
00032 typedef interpolation_type::point_type point_type;
00033 typedef F77GFMFileOutput<VectorType,FixupType,FlagType,DIM> output_type;
00034 public:
00035 SolverSpecific(IntegratorSpecific& integ,
00036 base::initial_condition_type& init,
00037 base::boundary_conditions_type& bc) : base(integ, init, bc) {
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 }
00051
00052 ~SolverSpecific() {
00053 DeleteGFM(_GFM[0]);
00054 delete _Flagging;
00055 delete _Fixup;
00056 delete _Interpolation;
00057 }
00058
00059 virtual void SetupData() {
00060 base::SetupData();
00061 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00062 }
00063
00064
00065
00066
00067 virtual double Tick(int VariableTimeStepping, const double dtv[],
00068 const double cflv[], int& Rejections) {
00069
00070 int me = MY_PROC;
00071
00072 point_type force_old[5];
00073 double t_old = base::t[0];
00074 CalculateForce(force_old);
00075
00076 double cfl = base::Tick(VariableTimeStepping, dtv, cflv, Rejections);
00077
00078 point_type a, force_new[5];
00079 double dt = base::t[0]-t_old;
00080 CalculateForce(force_new);
00081
00082 int NSph, Np[5];
00083 double xm[15], vm[15], cm[5], rd[5];
00084 f_combl_read(NSph,xm,vm,cm,rd,Np);
00085
00086 for (int ns=0; ns<NSph; ns++) {
00087 a = 0.5/cm[ns]*(force_old[ns]+force_new[ns]);
00088 vm[base::Dim()*ns ] += a(0)*dt;
00089 vm[base::Dim()*ns+1] += a(1)*dt;
00090 vm[base::Dim()*ns+2] += a(2)*dt;
00091 xm[base::Dim()*ns ] += vm[base::Dim()*ns ]*dt;
00092 xm[base::Dim()*ns+1] += vm[base::Dim()*ns+1]*dt;
00093 xm[base::Dim()*ns+2] += vm[base::Dim()*ns+2]*dt;
00094 }
00095
00096 f_combl_write(NSph,xm,vm,cm,rd,Np);
00097
00098
00099 if (me == VizServer) {
00100 std::ofstream outfile;
00101 std::ostream* out;
00102 std::string fname = "drag.dat";
00103 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00104 out = new std::ostream(outfile.rdbuf());
00105 *out << base::t[0];
00106 for (int ns=0; ns<NSph; ns++) {
00107 int d;
00108 for (d=0; d<base::Dim(); d++)
00109 *out << " " << force_new[ns](d);
00110 for (d=0; d<base::Dim(); d++)
00111 *out << " " << xm[base::Dim()*ns+d];
00112 for (d=0; d<base::Dim(); d++)
00113 *out << " " << vm[base::Dim()*ns+d];
00114 }
00115 *out << std::endl;
00116 outfile.close();
00117 delete out;
00118 }
00119
00120 return cfl;
00121 }
00122
00123
00124 void CalculateForce(point_type* sum) {
00125 int me = MY_PROC;
00126
00127 double pi = 4.0*std::atan(1.0);
00128
00129
00130 int NSph, Np[5];
00131 double xm[15], vm[15], cm[5], rd[5];
00132 f_combl_read(NSph,xm,vm,cm,rd,Np);
00133
00134
00135 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00136 int Time = CurrentTime(base::GH(),l);
00137 int press_cnt = base::Dim()+4;
00138 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00139 press_cnt, base::t[l]);
00140 }
00141
00142 for (int ns=0; ns<NSph; ns++) {
00143 int Npoints = 8*Np[ns]*Np[ns];
00144 DataType* p = new DataType[Npoints];
00145 point_type* xc = new point_type[Npoints];
00146 point_type* normals = new point_type[Npoints];
00147 DataType* ds = new DataType[Npoints];
00148 DataType* dist = (DataType*) 0;
00149 double dalpha = 0.5*pi/Np[ns];
00150 double dbeta = 0.5*pi/Np[ns];
00151 register int n = 0;
00152
00153
00154 for (register int l=0; l<4*Np[ns]; l++) {
00155 double alpha = dalpha*(l+0.5);
00156 for (register int m=0; m<2*Np[ns]; m++) {
00157 double beta = dalpha*(m+0.5);
00158 xc[n](0) = xm[base::Dim()*ns ]+rd[ns]*std::cos(alpha)*std::sin(beta);
00159 xc[n](1) = xm[base::Dim()*ns+1]+rd[ns]*std::sin(alpha)*std::sin(beta);
00160 xc[n](2) = xm[base::Dim()*ns+2]+rd[ns]*std::cos(beta);
00161 ds[n] = pi*rd[ns]*rd[ns]*dbeta*std::sin(beta)/(2.*Np[ns]);
00162 n++;
00163 }
00164 }
00165
00166
00167 _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
00168 _Interpolation->ArrayCombine(VizServer,Npoints,p);
00169
00170
00171 _Interpolation->PointsGeom(base::GFM(0).Phi(),base::t[0],Npoints,xc,normals,dist);
00172 _Interpolation->ArrayCombine(VizServer,base::Dim()*Npoints,normals[0].data());
00173
00174
00175 if (me == VizServer) {
00176 sum[ns] = 0.0;
00177 for (n=0; n<Npoints; n++)
00178 if (p[n]>-1.e37)
00179 sum[ns] += p[n]*normals[n]*ds[n];
00180 }
00181
00182 delete [] p;
00183 delete [] xc;
00184 delete [] normals;
00185 delete [] ds;
00186 }
00187
00188 #ifdef DAGH_NO_MPI
00189 #else
00190 if (comm_service::dce() && comm_service::proc_num() > 1) {
00191 int num = comm_service::proc_num();
00192 MPI_Status status;
00193 int R;
00194 int size = sizeof(DataType)*5;
00195 int tag = 10001;
00196 if (me == VizServer)
00197 for (int proc=0; proc<num; proc++) {
00198 if (proc==VizServer) continue;
00199 R = MPI_Send((void *)sum, size, MPI_BYTE, proc, tag, comm_service::comm());
00200 if ( MPI_SUCCESS != R )
00201 comm_service::error_die( "SumCombine", "MPI_Send", R );
00202 }
00203 else {
00204 R = MPI_Recv((void *)sum, size, MPI_BYTE, VizServer, tag, comm_service::comm(), &status);
00205 if ( MPI_SUCCESS != R )
00206 comm_service::error_die( "ArrayCombine", "MPI_Recv", R );
00207 }
00208 }
00209 #endif
00210 }
00211
00212 protected:
00213 interpolation_type* _Interpolation;