LSProblemImplem.H

Go to the documentation of this file.

#ifdef CH_LANG_CC
/*
*      _______              __
*     / ___/ /  ___  __ _  / /  ___
*    / /__/ _ \/ _ \/  V \/ _ \/ _ \
*    \___/_//_/\___/_/_/_/_.__/\___/
*    Please refer to Copyright.txt, in Chombo's root directory.
*/
#endif

#ifndef _LSPROBLEMIMPLEM_H_
#define _LSPROBLEMIMPLEM_H_

#if defined(CH_Darwin) && defined(__GNUC__) && ( __GNUC__ == 3 )
// deal with the broken isnan()/isinf() in GCC on MacOS
#include <unistd.h>
#define _GLIBCPP_USE_C99 1
#endif

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <string>
using std::endl;

#include "ConstrainedLS.H"
#include "LSquares.H"
#include "IFData.H"
#include "CutCellMoments.H"

#include "NamespaceHeader.H"

template<int dim>
int LSProblem<dim>::nChooseR(int n, int r)
{
  if (r == 0) return 1;
  int num = 1;
  int den = 1;
  for (int i = 1; i <=r ; ++i) {
    num *= (n+1-i);
    den *= i;
  }
  return num/den;
}




template<int dim> void LSProblem<dim>::computeBounds(const IndexTM<Real,dim> & a_dx,
                                                     const CutCellMoments<dim> & a_ccm)
{
  int nEBBounds = 0;
  if (m_degreeP == 0 && dim == 2) // || m_degreeP == 1)
    {
      nEBBounds = m_numP;       // A lower bound for each EB. Haven't figured out a useful u.b. yet
    }


  int nVar=m_numP + m_numPLess1; // total variables in problem, the first

  m_nInequality= nEBBounds+2*m_numPLess1;
  allocArray(nVar,m_nInequality,m_constraintMat);
  m_constraintVec = new Real[m_nInequality];

  // This constraint is a lower bound on the zero monomial of the EB.
  // It is based on the fact that integral(mono)/integral(1) < Max(mono)
  if (m_degreeP == 0 && dim == 2)
    {
      Real lobnd = 0.; // minimum, but we can do better
      IvDim whichMono; // for tracking
      for(typename CutCellMoments<dim>::PthMoment::const_iterator it = a_ccm.m_EBmoments.begin();
          it != a_ccm.m_EBmoments.end();++it)
        {
          IvDim mono = it->first;
          Real val = it->second;
          if (val < 0) val = -val;
          //pout() << "In zero degree, " << mono <<  " val=" << val <<endl;
          for (int idim=0;idim < dim;++idim){
            val /= pow(a_ccm.m_iFData.m_dx[idim]/2.,mono[idim]);
          }
          if (val > lobnd){
            whichMono = mono;
            lobnd = (val> lobnd) ? val : lobnd;
          } //todo: this wasfor debug
        }
      //pout() << "New EB low bound:: " << lobnd << ", " <<whichMono << endl;

      // The constraints matrix fits CiT x + ci >=0
      // so the rows are variables and the columns are constraints (may be a FORTRAN holdover)
      IvDim mono = IvDim::Zero;

      //pout() << "Monomial: " << mono << a_dx << endl;
      int ndx = 0; //todo: guarantee there should be only one for (0,0)
      int variableNdx = ndx;
      int loBoundNdx = ndx;
      m_constraintMat[variableNdx][loBoundNdx] = 1.0;
      m_constraintVec[loBoundNdx]=-lobnd;
    }

  // Now create constraints for the volume integrals.
  Vector<Real> bounds(2*m_monoLocPLess1.size());
  for(typename PthMomentLoc::const_iterator it = m_monoLocPLess1.begin();
      it != m_monoLocPLess1.end();++it)
    {

      // The constraints matrix fits CiT x + ci >=0
      // so the rows are variables and the columns are constraints (may be a FORTRAN holdover)
      IvDim mono = it->first;

      //pout() << "Monomial: " << mono << a_dx << endl;
      int ndx = it->second;
      int variableNdx = ndx + m_numP;
      int loBoundNdx = nEBBounds + ndx*2;
      int hiBoundNdx = loBoundNdx+1;
      Real lobnd, hibnd;

      // Compute the lower and upper bounds
      // based on quadrature of all the regions that contribute
      // negatively or positively.
      momentBounds(lobnd,hibnd,mono, a_dx);

      // Refine the lower bound of the zero degree monomial 
      // based on the identity
      // integral(mono)/integral(zeromono) < Max(mono)
      bool isZeroDegree = true;
      for (int idim = 0; idim < dim; ++idim){
        if (mono[idim] >0) isZeroDegree =false;
      }
      if (isZeroDegree){
        IvDim whichMono;
        for(typename CutCellMoments<dim>::PthMoment::const_iterator it = a_ccm.m_moments.begin();
              it != a_ccm.m_moments.end();++it)
            {
              IvDim mono = it->first;
              Real val = abs(it->second);
              int degree=0;
              //pout() << "In zero degree, " << mono <<  " val=" << val <<endl;
              for (int idim=0;idim < dim;++idim){
                val /= pow(a_ccm.m_iFData.m_dx[idim]/2.,mono[idim]);
                degree+=mono[idim];
              }
              if (val > lobnd && degree < 2){
                whichMono = mono;
                lobnd = Max(val, lobnd);
                //pout() << "whichMono="<<whichMono<<endl;
              }
            }

        // Refine the maximum volume or area based on the no undercut 
        // assumption, which is appropriate only for "top face" 
        // calculations such as shallow water. 
        // No undercut means that the volume in 2D or 3D is 
        // constrained by the top area times a maximum plausible height
        // Originally this maximum was the reference height, but now it is the 
        // distance from the lowest z intersection to the reference height.
        if ( dim == 3 || (dim == 2 && a_ccm.m_iFData.m_normalDir[0] != 2)){
          // first compute the top face area.
          int zdir = dim - 1;    // todo: verify/guarantee this is z direction
          Iv2 bndMoment;         // for fetching a moment
          bndMoment[0] = zdir;  
          bndMoment[1] = 1;      // this means on the hiSide
          //pout() << "Area boundary, dim=" << dim << endl;
          const CutCellMoments<dim-1> & refSide = a_ccm.getBdCutCellMoments(bndMoment);
          Real areaBnd = refSide.getMoment(IndexTM<int,dim-1>::Zero,VolMoment);
          //pout() <<"dx (dim=" << dim << ")" << a_ccm.m_iFData.m_dx << endl;
          //pout() << "Dim=" << dim <<" Area boundary: "<< areaBnd << endl;
          //pout() << "Ref side: " << refSide << endl;
          // now compute the maximum plausible height
          Real refHeight = a_ccm.m_iFData.m_dx[zdir];
          Real maxHeight = refHeight;
          int nintersect = 0;
          if (dim == 3)
            {
              for (typename IFData<dim>::EdgeIntersections::const_iterator 
                     it = a_ccm.m_iFData.m_intersections.begin();
                   it != a_ccm.m_iFData.m_intersections.end(); ++ it)
                {
                  if (it->first[0] == 2)
                    { 
                      //pout()<< it->first << " " << it->second << endl;
                      Real zIntersectFrac = it->second;
                      Real zIntersectHt = refHeight * (1.0 - zIntersectFrac)/2.0;
                      maxHeight = Min(maxHeight,zIntersectHt);
                      nintersect++;
                      //pout() << "Max Height revised: " << maxHeight << endl;
                    }
                }
            }

          //pout() << "Ref height: " << refHeight <<endl;
          //pout() <<"Area normal to ref height (currently hardwired): " <<areaBnd <<endl;
          areaBnd *=  refHeight;
          if (nintersect > 1)areaBnd *= 1.0;

          if (areaBnd < 0.0){
            MayDay::Abort( "Area bound negative");
          }
          hibnd = Min(areaBnd,hibnd);
          if (hibnd < lobnd)lobnd=hibnd*0.9;
          if (hibnd < lobnd){
            pout() << "Lower Bound: " << lobnd << " Upper Bound: " << hibnd << endl;
            MayDay::Abort( "Lower bound is higher than high bound");
          }
          //pout() << "New high bound for " << mono <<": " << hibnd << endl;
        }
      }
      m_constraintMat[variableNdx][loBoundNdx] = 1.0;
      m_constraintMat[variableNdx][hiBoundNdx] = -1.0;
      m_constraintVec[loBoundNdx]=-lobnd;
      m_constraintVec[hiBoundNdx]=hibnd;

      //pout() << "Low: " << lobnd << ".....Hi" << hibnd << endl;

    }


}

//Calculate upper and lower bounds for a moment
template<int dim>
void LSProblem<dim>::momentBounds(Real & lobnd,
                                  Real & hibnd,
                                  const IvDim & mono,
                                  const IndexTM<Real,dim>& a_dx)
{

  // Calculate the contribution to the monomial bound of the 
  // even degrees. This contribution will modulate the odd part.
  // At the same time, we accumulate a list of odd dimensions
  Real evenbnd = 1.;
  Vector<int> oddDims;
  for (int idim = 0; idim < dim ; ++idim){
    if (mono[idim] % 2 == 0){
      Real pPlus1 = (Real) mono[idim]+1;
      evenbnd *= pow(((Real)0.5)*a_dx[idim],pPlus1)/pPlus1;
    }
    else
      {
        oddDims.push_back(idim);
      }
  }

  int nOddDims = oddDims.size();
  int nEvenDims = dim - nOddDims;
  evenbnd *= pow(2.0, nEvenDims);

  // integrate the odd powers from middle to edge of cell, assuming middle is zero
  Real partialCellMag = 1;
  for (int iodd = 0; iodd < nOddDims ; iodd++){
    Real pPlus1 =(Real)( mono[oddDims[iodd]]+1);
    partialCellMag *= pow(((Real)0.5)*a_dx[oddDims[iodd]],pPlus1)/pPlus1;
  }


  // Now figure out number of sections (quadrants, octants...)
  // that will contribute positively because they 
  // yield an even number of negative numbers
  Real posContrib=0.0;
  Real negContrib=0.0;

  if (nOddDims > 0){
    for (int j=0; j <= nOddDims; j+=2){
      int nContrib = nChooseR(nOddDims, j);
      posContrib+= nContrib* partialCellMag;
    }

    for (int j=1; j <= nOddDims; j+=2){
      int nContrib = nChooseR(nOddDims, j);
      negContrib-= nContrib* partialCellMag;
    }
  }
  hibnd =  nOddDims > 0 ? posContrib*evenbnd : evenbnd;
  lobnd = nOddDims > 0 ? negContrib*evenbnd : 0.0;

}



template<int dim> LSProblem<dim>::~LSProblem()
{
  if (m_matrix != NULL)
    {
      //free m_matrix
      int numRows = dim*m_numP;
      int numCols = m_numP + m_numPLess1;
      freeArray(numRows,numCols,m_matrix);

    }
  if (m_constraintMat != NULL){
      freeArray(m_numP+m_numPLess1, m_nInequality,m_constraintMat);
      delete [] m_constraintVec;
  }

}
//this constructor is used when just a list of monomials is wanted
template<int dim> LSProblem<dim>::LSProblem(const int  & a_degreeP,
                                            const bool & a_useConstraints)
  :m_degreeP(a_degreeP),
   m_useConstraints(a_useConstraints)
{
   fillMap(m_monoLocP,m_locMonoP,m_degreeP);
   fillMap(m_monoLocPLess1,m_locMonoPLess1,m_degreeP-1);
   m_numP      = numMonomials(m_degreeP);
   m_numPLess1 = numMonomials(m_degreeP-1);
   m_matrix = NULL;
   m_constraintMat = NULL;
   m_constraintVec = NULL;
}
//constructor for solving a LS problem
template<int dim> LSProblem<dim>::LSProblem(const int              & a_order,
                                            const int              & a_degreeP,
                                            const bool             & a_useConstraints,
                                            const IndexTM<Real,dim>& a_normal)
  :m_order(a_order),
   m_degreeP(a_degreeP),
   m_useConstraints(a_useConstraints),
   m_normal(a_normal)
{
  fillMap(m_monoLocP,m_locMonoP,m_degreeP);
  fillMap(m_monoLocPLess1,m_locMonoPLess1,m_degreeP-1);
  m_numP      = numMonomials(m_degreeP);
  m_numPLess1 = numMonomials(m_degreeP-1);
  setMatrix();

  m_unknowns.resize(m_numP + m_numPLess1);
  m_rhs.resize(m_numP*dim);
  m_constraintMat = NULL;
  m_constraintVec = NULL;

}

//
template<int dim> void LSProblem<dim>::invertNormalEq(const Vector<Real>& a_rhs,Vector<Real>& a_residual)
{
  m_rhs = a_rhs;

  if (!m_useConstraints)
  {
    LSquares whichMethod;
    whichMethod.LeastSquares(m_matrix,m_unknowns,m_rhs);
  }
  else
  {
    ConstrainedLS whichMethod;
    whichMethod.LeastSquares(m_matrix,m_unknowns,m_rhs,m_constraintMat,m_constraintVec,m_nInequality);
  }

  a_residual.resize(3);
  a_residual[0] = 0.0;
  a_residual[1] = 0.0;
  a_residual[2] = 0.0;
  Real maxRi = 0.0;
  for (int i = 0 ; i < m_numP*dim ; i++)
    {
      Real AtimeX = 0.0;
      for (int j = 0 ; j < m_numP + m_numPLess1 ; j++)
        {
          AtimeX += m_matrix[i][j] * m_unknowns[j];
        }
      Real ri = AtimeX - m_rhs[i];
      if (Abs(ri) > maxRi)
        {
          a_residual[0] = Abs(ri);
          maxRi = Abs(ri);
        }
      a_residual[1] += Abs(ri);
      a_residual[2] += ri * ri;
    }
  a_residual[2] = sqrt(a_residual[2]);
  //pout() << "Residual : " << a_residual[1] << ":" << a_residual[2] << endl;
}


template<int dim> int LSProblem<dim>::factorial(const int &a_n,const int &a_m)
{
  int retval = 1;
  if (a_n < 0 || a_m < 0)
    {
      MayDay::Abort("Attempting n! for n < 0");
    }
  for (int i = a_m+1; i <= a_n; ++i)
    {
      retval *= i;
    }
  return retval;
}

template<int dim> int LSProblem<dim>::numMonomials(const int& a_monoDegree)
{
  int retval = LARGEINTVAL;
  if (a_monoDegree == -1)
    {
      retval =  0;
    }
  else
    {
      int bigger;
      int smaller;
      if (dim- 1 > a_monoDegree)
        {
          bigger = dim-1;
          smaller = a_monoDegree;
    }
      else
        {
          smaller = dim-1;
          bigger = a_monoDegree;
        }
      int numerator = factorial(dim - 1 + a_monoDegree,bigger);
      int denominator = factorial(smaller);
      // pout()<<"numerator = "<<numerator<<endl;
      //pout()<<"denominator = "<<denominator<<endl;
      retval =  numerator/denominator;
    }
  return retval;
}


//uses dimension & degree create matrix for overdetermined system
template<int dim> void LSProblem<dim>::setMatrix()
{
  //solving m_matrix[x] = b

  //initializes m_matrix to zeros


  int numRows = dim*m_numP;
  int numCols = m_numP + m_numPLess1;
  allocArray(numRows,numCols,m_matrix);

  //iterate through the list of mono of Degree P
  //pout()<<"Number of mono of DegreeP = "<<m_monoLocP.size()<<endl;
  for(typename PthMomentLoc::const_iterator it = m_monoLocP.begin();
      it != m_monoLocP.end();++it)
    {
      for (int idir = 0; idir< dim; ++idir)
        {
          //this entry corresponds to the integral over the boundary
          int row = dim*(it->second) + idir;
          int pcol = it->second;
          m_matrix[row][pcol] = -m_normal[idir];

          //differentiate the mono and
          IvDim Dmono;
          IvDim mono = it->first;
          int coeff = LARGEINTVAL;
          //diff(it->first) = coeff*Dmono
          differentiate(coeff,Dmono,idir,mono);
          int pLess1Col= LARGEINTVAL;
          if (coeff == LARGEINTVAL)
            {
              MayDay::Abort("problem wth differentiate");
            }
          if (coeff > 0)
            {
              //find which mono this is in the list
              if (m_monoLocPLess1.find(Dmono) != m_monoLocPLess1.end())
                {
                  pLess1Col = m_monoLocPLess1[Dmono] + m_numP;
                  m_matrix[row][pLess1Col] = coeff;
                }
              else
                {
                  MayDay::Abort("can't find derived mono");
                }
            }
        }
    }
}





//differentiate a_mono w.r.t. x_idir. Answer = a_coeff*a_Dmono
template<int dim> void LSProblem<dim>::differentiate(int         & a_coeff,
                                                     IvDim       & a_Dmono,
                                                     int         & a_idir,
                                                     const IvDim & a_mono)
{
  if(a_mono[a_idir] > 0)
    {
      a_coeff          = a_mono[a_idir];
      a_Dmono          = a_mono;
      a_Dmono[a_idir] -= 1;
    }
  else if (a_mono[a_idir] == 0)
    {
      a_coeff = 0;
      for(int idir = 0; idir < dim; ++idir)
        {
          a_Dmono[idir] = LARGEINTVAL;
        }
    }
  else
    {
      MayDay::Abort("Monomial has negative power");
    }
}
template<int dim>void LSProblem<dim>::setRhs(const Vector<Real>& a_rhs)
{
}

template<int dim> void LSProblem<dim>::fillMap(PthMomentLoc& a_monoLoc,
                                               LocPthMoment& a_locMono,
                                               const int&     a_degree)
{
  if (a_degree >= 0)
    {
      //indexes into the map
      int index = 0;
      IndexTM<int,dim>monomial;
      for (int idir = 0; idir < dim; ++ idir)
        {
          monomial[idir] = 0;
        }
      while(true)
        {
          for (int j = 1;j < dim-1; ++j)
            {
              int t = a_degree;
              for(int k = j+1; k < dim; ++k)
                {
                  t -= monomial[k];
                }

              if (monomial[j] >  t)
                {
                  monomial[j+1] += 1;
                  monomial[j] = 0;
                }
              else
                {
                  break;
                }
            }
          if ( monomial[dim -1] > a_degree)
            {
              break;
            }
          monomial[0] = a_degree;
          for(int j=1; j < dim; ++ j)
            {
              monomial[0] -= monomial[j];
            }

          //pout()<<"monomial = "<<monomial<<endl;
          a_monoLoc[monomial] = index;
          a_locMono[index] = monomial;

          index +=1;
          monomial[1] += 1;
        }
    }
}
template<int dim> void LSProblem<dim>::allocArray(const int& rows,
                                                  const int& cols,
                                                  Real**& A)
{
  A = new Real* [rows];

  for (int i = 0; i < rows;i++)
    {
      A[i] = new Real [cols];
      Real* scanA = A[i];

      for (int j = 0; j < cols; j++)
        {
          *(scanA++) = 0.0;
        }
    }
}

template<int dim> void LSProblem<dim>::freeArray(const int& rows,
                         const int& cols,
                         Real**& A)
{
  for (int i = 0; i < rows; i++)
    {
      delete[] A[i];
    }

  delete[] A;
}
template<int dim > void LSProblem<dim>::print(ostream& out) const
{

  out <<"Dim = "<< dim<<", degree = "<< m_degreeP<<'\n';
  out <<"m_monoLocP has " << m_monoLocP.size()<<" elements"<< '\n';

   for(typename PthMomentLoc::const_iterator it = m_monoLocP.begin();
       it != m_monoLocP.end();++it)
  {
    pout()<<"Monomial =  "<<it->first<<", Loc = "<< it->second << '\n';
  }

    out <<"Dim = "<< dim<<'\n';
  out <<"m_locMonoP has " << m_locMonoP.size()<<" elements"<< '\n';

   for(typename LocPthMoment::const_iterator it = m_locMonoP.begin();
       it != m_locMonoP.end();++it)
  {
    pout()<<"Loc =  "<<it->first<<", Monomial = "<< it->second << '\n';
  }
   //one degree lower
   out <<"m_locMonoPLess1 has " << m_locMonoPLess1.size()<<" elements"<< '\n';
   for(typename PthMomentLoc::const_iterator it = m_monoLocPLess1.begin();
       it != m_monoLocPLess1.end();++it)
     {
    pout()<<"Monomial =  "<<it->first<<", Loc = "<< it->second << '\n';
  }

  out <<"m_locMonoPLess1 has " << m_locMonoPLess1.size()<<" elements"<< '\n';

   for(typename LocPthMoment::const_iterator it = m_locMonoPLess1.begin();
       it != m_locMonoPLess1.end();++it)
  {
    pout()<<"Loc =  "<<it->first<<", Monomial = "<< it->second << '\n';
  }

   pout()<<"Matrix and rhs for least squares problem of dim = "<<dim<<endl;
   outputMatrix();
   outputRhs();
   outputUnknowns();
}

template<int dim> ostream& operator<<(ostream& a_out, LSProblem<dim>& a_lsProblem)
{
  a_lsProblem.print(a_out);
  return a_out;
}
template<int dim> void LSProblem<dim>::outputMatrix()const
{
  int rows = m_locMonoP.size()*dim;
  int cols = m_locMonoPLess1.size() + m_locMonoP.size();
  pout()<<"numRows = "<< rows<<", numCols = "<< cols<<endl;
  // pout() << "outputting " << name << endl;
  for (int i = 0; i < rows; i++)
    {
      for (int j = 0; j < cols; j++)
       {
         pout()<< m_matrix[i][j] << " ";
       }
      pout()<<endl;
    }
}
template<int dim> void LSProblem<dim>::outputRhs()const
{
  pout()<<"Outputting Rhs"<<endl;
  for (int i = 0; i < m_rhs.size(); i++)
    {
      pout()<<"m_rhs["<<i<<"] = "<<m_rhs[i]<<endl;
    }
}

template<int dim> void LSProblem<dim>::outputUnknowns()const
{
  pout()<<"Outputting Rhs"<<endl;
  for (int i = 0; i < m_unknowns.size(); i++)
    {
      pout()<<"m_unknowns["<<i<<"] = "<<m_unknowns[i]<<endl;
    }
}
#include "NamespaceFooter.H"

#endif

---