clawpack/applications/euler_chem/2d/ConvWedge/src/Problem.h

// -*- C++ -*-

// Copyright (C) 2003-2007 California Institute of Technology
// David Hill, Ralf Deiterding

#ifndef AMROC_PROBLEM_H
#define AMROC_PROBLEM_H

#include "eulerrhok2.h"

#define NEQUATIONS  7
#define NEQUSED     7
#define NFIXUP      6
#define NAUX        1

#include "ClpProblem.h"

#define f_track FORTRAN_NAME(track, TRACK)


extern "C" {
  void f_track(const INTEGER &, const DOUBLE *, const DOUBLE *, DOUBLE *, INTEGER&,
               const DOUBLE&, const DOUBLE&, const DOUBLE&,
                DOUBLE&, DOUBLE&, DOUBLE&, DOUBLE&, DOUBLE&, DOUBLE&); 
}

#define f_getangles FORTRAN_NAME(getangles, GETANGLES)

extern "C" {
  void f_getangles(DOUBLE&, DOUBLE&, DOUBLE&); 
}

#define OWN_FLAGGING
#define OWN_GFMAMRSOLVER
#include "ClpStdGFMProblem.h"
#include "AMRGFMInterpolation.h"
#include "F77Interfaces/F77UpdateLevelTransfer.h"
#include "SpecialOutput/AMRGridAccumulation.h"


class FlaggingSpecific : 
  public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
  typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;
public:
  FlaggingSpecific(solver_type& solver) : base(solver) {
      base::AddCriterion(new ScaledGradient<VectorType,FlagType,DIM>());
      base::AddCriterion(new LimiterType<VectorType,FlagType,DIM>());
      base::AddCriterion(new AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver));
      base::AddCriterion(new RelativeError<VectorType,FixupType,FlagType,DIM>(solver));
      base::AddCriterion(new F77ScaledGradient<VectorType,FlagType,DIM>(f_flg));
      base::AddCriterion(new F77LimiterType<VectorType,FlagType,DIM>(f_flg));
      base::AddCriterion(new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,f_flg));
      base::AddCriterion(new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,f_flg));
    }
  ~FlaggingSpecific() { DeleteAllCriterions(); }
};

class SolverSpecific : 
  public AMRGFMSolver<VectorType,FixupType,FlagType,DIM> {
  typedef VectorType::InternalDataType DataType;
  typedef AMRGFMSolver<VectorType,FixupType,FlagType,DIM> base;
  typedef AMRGFMInterpolation<VectorType,FixupType,FlagType,DIM> interpolation_type;
  typedef interpolation_type::point_type point_type;
  typedef F77GFMFileOutput<VectorType,FixupType,FlagType,DIM> output_type;
public:
  SolverSpecific(IntegratorSpecific& integ, 
                 base::initial_condition_type& init,
                 base::boundary_conditions_type& bc) : base(integ, init, bc) {
    SetLevelTransfer(new F77UpdateLevelTransfer<VectorType,DIM>(f_prolong, f_restrict, f_tupdate));
    SetFileOutput(new F77GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this,f_out)); 
    SetFixup(new FixupSpecific(integ));
    SetFlagging(new FlaggingSpecific(*this)); 
    AddGFM(new GhostFluidMethod<VectorType,DIM>(
              new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans),
              new F77GFMLevelSet<DataType,DIM>(f_lset)));
    _Interpolation = new interpolation_type(*this);
  }  
 
  ~SolverSpecific() {
    DeleteGFM(_GFM[0]);
    delete _LevelTransfer;
    delete _Flagging;
    delete _Fixup;
    delete _FileOutput;
    delete _Interpolation;
  }

  virtual void SetupData() {
    base::SetupData();
    _Interpolation->SetupData(base::PGH(), base::NGhosts());
  }

  //  use the Tick function to take pressure readings at
  //  locations that corespond to the transducers in the experiment
  virtual double Tick(int VariableTimeStepping, const double dtv[], const double cflv[],
                      int& Rejections) {

    double cfl = base::Tick(VariableTimeStepping, dtv, cflv, Rejections);

    int me = MY_PROC;
    
    double pi = 4.0*std::atan(1.0);
    // the wegde angles in radians
    // upper plate
    double alpha1; // = 13.9576*pi/180.;
    // lower plate
    double alpha2; // = 9.45*pi/180.;
    double radi; // hinge radius in meters = .25*0.254
    double ra; // hinge radius in inches
    // use fortran function that sees the cuser.i
    // to pull out the wedge angles and hinge radius
    f_getangles(alpha1,alpha2,radi);
    // convert to inches
    ra = radi/0.0254;
    // save trig results
    double c1 = std::cos(alpha1);
    double s1 = std::sin(alpha1);
    double c2 = std::cos(alpha2);
    double s2 = std::sin(alpha2);
    // get the finest grid spacing
    int finelevel = FineLevel(base::GH());
    double dx = DeltaX(base::GH(),0,finelevel);
    double dy = DeltaX(base::GH(),1,finelevel);
    double dl;
    dl = std::sqrt(dx*dx+dy*dy);

    // seven transducers bottom plate
    // and seven on top on wedge plate
    int Npoints = 14;
    // an array for the interpolated pressure
    DataType* p = new DataType[Npoints];
    // an array for the radial distance to transducers
    DataType* r = new DataType[Npoints];
    // an array for the locations of the transducers
    point_type* xc = new point_type[Npoints];
    
    // radial transducer locations in inches
    // measures from hinge outward.  for distance to
    // apex add 1/4 inch.
    //

//  old distances - wedge was altered
//     // load the bottom plate locations first
//     r[0] = 9.5;      r[1] = 7.0;   r[2] = 4.25;
//     r[3] = 2.75;     r[4] = 1.5;   r[5] = 0.6;
//     r[6] = 0.25;
//     //load the top plate locations from tip to apex
//     r[7] = 10.72;    r[8] = 6.0;   r[9] = 3.5;
//     r[10] = 2.0;     r[11] = 0.75; r[12] = 0.4 ;
 
// new distances
    // load the bottom plate locations first
    r[0] = 9.3;      r[1] = 7.0;   r[2] = 4.25;
    r[3] = 2.75;     r[4] = 1.5;   r[5] = 0.6;
    r[6] = 0.25;
    //load the top plate locations from tip to apex
    r[7] = 9.3; 
    r[8] = 8.0;    r[9] = 6.0;   r[10] = 3.5;
    r[11] = 2.0;     r[12] = 0.75; r[13] = 0.5 ;

    // create the x/y locations of the transducers 
    // relative to apexin meters. 
    // offset locations by 'dl' into fluid so we don't 
    // use ghost values when interpolating for the pressure
    //
    //bottom plate
    register int n ;
    for(n=0; n<=6; n++) {
      xc[n](0) =  -dl*s2-0.0254*c2*(r[n]+ra);
      xc[n](1) = dl*c2-0.0254*s2*(r[n]+ra);
    }
    //top plate
     for(n=7; n<=13; n++) {
      xc[n](0) =   -dl*s1 -0.0254*c1*(r[n]+ra);
      xc[n](1) = -dl*c1+0.0254*s1*(r[n]+ra);
    }

    
    // get the pressure from the conserved vector of state
    for (register int l=0; l<=FineLevel(base::GH()); l++) {
      int Time = CurrentTime(base::GH(),l);
      int press_cnt = base::Dim()+4;
      ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l, 
                                              press_cnt, base::t[l]);
    }
    // interpolate the pressure values
    _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
    _Interpolation->ArrayCombine(VizServer,Npoints,p);
    
    // output the pressures to one file (transducer)
    if (me == VizServer) {
      std::ofstream outfile;
      std::ostream* out;
      std::string fname = "transducer.dat";
      outfile.open(fname.c_str(), std::ios::out | std::ios::app);
      out = new std::ostream(outfile.rdbuf());
      // output the time,location,pressure
      *out << base::t[0];
      for( n=0; n<Npoints; n++) {
        *out << " " <<p[n];
      }
      *out << std::endl;
      outfile.close();
      delete out;      
    }
    
    delete [] p;
    delete [] r; 
    delete [] xc;
 
    // now do the shock speed along the centerline.
    

    // the probe diameter is .0016meters
    double y0rdi = 0.0016;
    
    double pmax = 0.0;
    // double rpmaxloc = -1.0;
    double shk_pos = -0.12;
    double reflection_time = 0.0010;
    double aMa,aMb,aMc,aMd,aMe,aMm;
    double fmin = 20000, fmax = 200000.;
    double dfmin = 0.1;

    // get the number of points to be used
    Npoints = base::shape[0]*RefinedBy(base::GH(),finelevel);
    
    // set up the new arrays
    DataType* pc = new DataType[Npoints];
    DataType* rc = new DataType[Npoints];
    DataType* xrc = new DataType[Npoints];
    point_type* xcc = new point_type[Npoints];
  
    for (n=0; n<Npoints; n++) {
      rc[n] =  base::geom[0]+Npoints*dx -(n+0.5)*dx;
      //r[n] = lngth/Npoints*n;
      if((rc[n])>radi) {
        xcc[n](0) = radi-dl/2.;
      } else {
        xcc[n](0) = rc[n];
      }
      xcc[n](1) = y0rdi;
    }
    
        // get the pressure from the conserved vector of state
    for (register int l=0; l<=FineLevel(base::GH()); l++) {
      int Time = CurrentTime(base::GH(),l);
      int press_cnt = base::Dim()+4;
      ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l, 
                                                press_cnt, base::t[l]);
    }

    // interpolate the pressure values
    _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xcc,pc);
    _Interpolation->ArrayCombine(VizServer,Npoints,pc);
    
    if (me == VizServer) {
      int nc;
      for (n=0; n<Npoints; n++) 
        if (pc[n]>pmax) {
          pmax = pc[n];
          // rpmaxloc = x0rdi-r[n];
          // rpmaxloc = xcc[n](0);
        }
      if (base::t[0]<reflection_time) {
        dfmin = 0.0000000001;
        f_track(Npoints,pc,rc,xrc,nc,fmin,fmax,dfmin,aMa,aMb,aMc,aMd,aMe,aMm);
        // 0 is the point with the largest x-coordinate (see indexing above)
        // shk_pos = x0rdi-xr[0];
        shk_pos = xrc[0];
      }
      else {
        dfmin = 0.1;
        f_track(Npoints,pc,rc,xrc,nc,fmin,fmax,dfmin,aMa,aMb,aMc,aMd,aMe,aMm);
        // nc-1 is the point with the smallest x-coordinate (see indexing above)
        //shk_pos = x0rdi-xr[nc-1];
        shk_pos = xrc[nc-1];
      }
    }
  
      
    delete [] pc;
    delete [] rc;
    delete [] xrc;      
    delete [] xcc;
    
    
    // done extracting the shock locations along the ray
    // now output the locations to the file shk.dat

    if (me == VizServer) {
      std::ofstream outfile;
      std::ostream* out;
      std::string fname = "shk.dat";
      outfile.open(fname.c_str(), std::ios::out | std::ios::app);
      out = new std::ostream(outfile.rdbuf());
      // output the time and the max pressure
      *out << base::t[0] << " " << shk_pos << " " << pmax;
      // also output the shock location
      *out << " " << aMa << " " << aMb << " " << aMc;
      *out << " " << aMd << " " << aMe << " " << aMm;
      *out << std::endl;
      outfile.close();
      delete out;      
    }



    return cfl ;
         
  }
protected:
  interpolation_type* _Interpolation;