multiadaptivetimeslab.cpp

利用C

  Alloc::realloc(&es, size_e.size, newsize);
  Alloc::realloc(&ee, size_e.size, newsize);
  Alloc::realloc(&ed, size_e.size, newsize);

  size_e.size = newsize;
}
//-----------------------------------------------------------------------------
void MultiAdaptiveTimeSlab::alloc_j(uint newsize)
{
  size_j.next = 0;

  if ( newsize <= size_j.size ) return;

  //message("Reallocating: nj = %d", newsize);

  Alloc::realloc(&jx, size_j.size, newsize);

  size_j.size = newsize;
}
//-----------------------------------------------------------------------------
void MultiAdaptiveTimeSlab::alloc_d(uint newsize)
{
  size_d.next = 0;

  if ( newsize <= size_d.size ) return;

  //message("Reallocating: nd = %d", newsize);

  Alloc::realloc(&de, size_d.size, newsize);

  size_d.size = newsize;
}
//-----------------------------------------------------------------------------
real MultiAdaptiveTimeSlab::computeEndTime(real a, real b, uint offset, uint& end)
{
  // Update partitition 
  real K = std::min(adaptivity.kmax(), b - a);
  K = partition.update(offset, end, adaptivity, K);

  //partition.debug(offset, end);
  
  // Modify time step if we're close to the end time
  if ( K < adaptivity.threshold() * (b - a) )
    b = a + K;

  // Save minimum time step
  kmin = std::min(kmin, b - a);
  
  return b;
}
//-----------------------------------------------------------------------------
real MultiAdaptiveTimeSlab::computeDataSize(real a, real b, uint offset)
{
  // Recursively compute data sizes using the same algorithm as
  // for the recursive creation of the time slab

  // Compute end time of this sub slab
  uint end = 0;
  b = computeEndTime(a, b, offset, end);

  // Use u to keep track of the latest time value for each component here
  for (uint n = offset; n < end; n++)
    u[partition.index(n)] = b;

  // Add contribution from this sub slab
  ns += 1;
  ne += end - offset;
  nj += method->nsize()*(end - offset);
  for (uint n = offset; n < end; n++)
    nd += countDependencies(partition.index(n));

  // Add contribution from all sub slabs
  real t = a;
  while ( t < b && end < partition.size() )
    t = computeDataSize(t, b, end);

  return b;
}
//-----------------------------------------------------------------------------
dolfin::uint MultiAdaptiveTimeSlab::countDependencies(uint i0)
{
  // Count the number of dependencies to components with smaller time steps
  // for the given component. This version is used before any elements are
  // created when we recursively compute the data size. We then use the
  // array u to store the latest time value for each component.

  uint n = 0;

  // Get list of dependencies for current component index
  const Array<uint>& deps = ode.dependencies[i0];
  
  // Iterate over dependencies
  for (uint pos = 0; pos < deps.size(); pos++)
  {
    // Get index of other component
    const uint i1 = deps[pos];
    
    // Use u to keep track of the latest time value for each component here
    if ( u[i0] > (u[i1] + DOLFIN_EPS) )
      n += method->nsize();
  }
  
  return n;
}
//-----------------------------------------------------------------------------
dolfin::uint MultiAdaptiveTimeSlab::countDependencies(uint i0, real b0)
{
  // Count the number of dependencies to components with smaller time steps
  // for the given component. This version is used at the time of creation
  // of elements and we may then get the time values of already created
  // elements.

  uint n = 0;

  // Get list of dependencies for current component index
  const Array<uint>& deps = ode.dependencies[i0];
  
  // Iterate over dependencies
  for (uint pos = 0; pos < deps.size(); pos++)
  {
    // Get index of other component
    const uint i1 = deps[pos];
    
    // Get last element for component
    const int e1 = elast[i1];
    
    // If we have not yet created the element, then it has not reached b0
    if ( e1 == -1 )
    {
      n += method->nsize();
      continue;
    }
    
    // Need to check end time value of element
    const uint s1 = es[e1];
    const real b1 = sb[s1];
    
    // Check if the component has reached b0
    if ( b1 < (b0 - DOLFIN_EPS) )
    {
      n += method->nsize();
    }
  }

  return n;
}
//-----------------------------------------------------------------------------
bool MultiAdaptiveTimeSlab::within(real t, real a, real b) const
{
  // Check if time is within the given interval, choosing the left interval
  // if we are close to the edge

  return (a + DOLFIN_EPS) < t && t <= (b + DOLFIN_EPS);
}
//-----------------------------------------------------------------------------
bool MultiAdaptiveTimeSlab::within(real a0, real b0, real a1, real b1) const
{
  // Check if [a0, b0] is contained in [a1, b1]

  return a1 <= (a0 + DOLFIN_EPS) && (b0 - DOLFIN_EPS) <= b1;
}
//-----------------------------------------------------------------------------
dolfin::uint MultiAdaptiveTimeSlab::coverSlab(int subslab, uint e0)
{
  // Start at e0 and step until we reach a new sub slab
  uint e = e0;
  for (; e < ne; e++)
  {
    // Check if we have reached the next sub slab
    if ( static_cast<int>(es[e]) != subslab )
      break;
    
    // Update elast 
    elast[ei[e]] = e;
  }

  // Return e1
  return e;
}
//-----------------------------------------------------------------------------
dolfin::uint MultiAdaptiveTimeSlab::coverNext(int subslab, uint element)
{
  // Check if we are still on the same sub slab
  if ( subslab == static_cast<int>(es[element]) )
    return subslab;

  // Get next sub slab
  subslab = es[element];
  
  // Update elast for all elements in the sub slab
  for (uint e = element; e < ne; e++)
  {
    // Check if we have reached the next sub slab
    if ( static_cast<int>(es[e]) != subslab )
      break;

    // Update elast 
    elast[ei[e]] = e;
  }

  return subslab;
}
//-----------------------------------------------------------------------------
void MultiAdaptiveTimeSlab::coverTime(real t)
{
  // Check if t is covered for all components
  bool ok = true;
  for (uint i = 0; i < N; i++)
  {
    // Get last covered element for the component
    const int e = elast[i];

    // Check if we need to start from the beginning
    if ( e == -1 )
    {
      emax = 0;
      ok = false;
      break;
    }

    // Get element data
    const uint s = es[e];
    const real a = sa[s];
    const real b = sb[s];

    // Check if we need to start from the beginning
    if ( t < (a + DOLFIN_EPS) )
    {
      emax = 0;
      ok = false;
      break;
    }

    // Check if we need to search forward, starting at e = emax
    if ( t > (b + DOLFIN_EPS) )
    {
      ok = false;
      break;
    }
  }

  // If ok is true, then a + DOLFIN_EPS <= t <= b + DOLFIN_EPS for all components
  if ( ok )
    return;

  // Reset sampling if necessary
  if ( emax >= ne )
    emax = 0;
  else
  {
    const uint s = es[emax];
    const real a = sa[s];
    
    if ( t < (a + DOLFIN_EPS) )
      emax = 0;
  }

  // Iterate forward until t is covered for all components
  for (uint e = emax; e < ne; e++)
  {
    // Get element data
    const uint s = es[e];
    const uint i = ei[e];
    const real a = sa[s];

    // Check if we have stepped far enough
    if ( t < (a + DOLFIN_EPS) && _a < (a - DOLFIN_EPS) )
      break;

    // Cover element
    elast[i] = e;
    emax = e;
  }
}
//-----------------------------------------------------------------------------
void MultiAdaptiveTimeSlab::cGfeval(real* f, uint s0, uint e0, uint i0, 
				    real a0, real b0, real k0)
{
  const uint& nn = method->nsize();
  const uint last = nn - 1;

  // Get list of dependencies for given component index
  const Array<uint>& deps = ode.dependencies[i0];

  // First evaluate at left end-point
  if ( a0 < (_a + DOLFIN_EPS) )
  {
    // Use previously computed value
    f[0] = f0[i0];
  }
  else
  {
    // Iterate over dependencies
    for (uint pos = 0; pos < deps.size(); pos++)
    {
      // Get other element
      const uint i1 = deps[pos];
      const int e1 = elast[i1];
      
      // Special case, component has no latest element
      if ( e1 == -1 )
      {
	u[i1] = u0[i1];
	continue;
      }
      
      // Three cases: k1 = k0, k1 < k0, k1 > k0
      const uint s1 = es[e1];
      if ( s1 == s0 )
      {
	// k1 = k0 (same sub slab)
	const int ep = ee[e1];
	const uint jp = ep * nn;
	u[i1] = ( ep != -1 ? jx[jp + last] : u0[i1] );
      }
      else
      {
	const real b1 = sb[s1];
	if ( b1 < (a0 + DOLFIN_EPS) )
	{
	  // k1 < k0 (smaller time step)
	  u[i1] = jx[e1 * nn + last];
	}
	else
	{
	  // k1 > k0 (larger time step)
	  const real a1 = sa[s1];
	  const real k1 = b1 - a1;
	  const real tau = (a0 - a1) / k1;
	  const int ep = ee[e1];
	  const uint jp = ep * nn;
	  const uint j1 = e1 * nn;
	  const real x0 = ( ep != -1 ? jx[jp + last] : u0[i1] );
	  u[i1] = method->ueval(x0, jx + j1, tau);
	}
      }
    }
    
    // Evaluate right-hand side
    f[0] = ode.f(u, a0, i0);
  }

  // Get first dependency to components with smaller time steps for element
  uint d = ed[e0];

  // Compute number of such dependencies for each nodal point
  const uint end = ( e0 < (ne - 1) ? ed[e0 + 1] : nd );
  const uint ndep = (end - d) / nn;
  dolfin_assert(ndep * nn == (end - d));

  // Evaluate the right-hand side at all quadrature points but the first
  for (uint m = 1; m < method->qsize(); m++)
  {
    // Compute quadrature point
    const real t = a0 + k0*method->qpoint(m);

    // Update values for components with larger or equal time steps
    for (uint pos = 0; pos < deps.size(); pos++)
    {
      // Get other element
      const uint i1 = deps[pos];
      const int e1 = elast[i1];
      
      // Special case, component has no latest element
      if ( e1 == -1 )
	continue;
      
      // Use fast evaluation for elements in the same sub slab
      const uint s1 = es[e1];
      const uint j1 = e1 * nn;
      if ( s0 == s1 )
      {
	u[i1] = jx[j1 + m - 1];
	continue;
      }

      // Skip components with smaller time steps
      const real b1 = sb[s1];
      if ( b1 < (a0 + DOLFIN_EPS) )
       	continue;
      
      // Interpolate value from larger element
      const real a1 = sa[s1];
      const real k1 = b1 - a1;
      const real tau = (t - a1) / k1;
      const int ep = ee[e1];
      const uint jp = ep * nn;
      const real x0 = ( ep != -1 ? jx[jp + last] : u0[i1] );
      u[i1] = method->ueval(x0, jx + j1, tau);
    }

    // Update values for components with smaller time steps
    for (uint dep = 0; dep < ndep; dep++)
    {
      // Get element
      const int e1 = de[d++];
      dolfin_assert(e1 != -1);

      // Get initial value for element
      const int ep = ee[e1];
      const uint i1 = ei[e1];
      const uint jp = ep * nn;
      const real x0 = ( ep != -1 ? jx[jp + last] : u0[i1] );
      
      // Interpolate value from smaller element
      const uint s1 = es[e1];
      const real a1 = sa[s1];
      const real b1 = sb[s1];
      const real k1 = b1 - a1;
      const real tau = (t - a1) / k1;
      const uint j1 = e1 * nn;
      u[i1] = method->ueval(x0, jx + j1, tau);
    }
    
    // Evaluate right-hand side
    f[m] = ode.f(u, t, i0);
  }
}
//-----------------------------------------------------------------------------
void MultiAdaptiveTimeSlab::dGfeval(real* f, uint s0, uint e0, uint i0, 
				  real a0, real b0, real k0)
{
  const uint& nn = method->nsize();

  // Get list of dependencies for given component index
  const Array<uint>& deps = ode.dependencies[i0];

  // Get first dependency to components with smaller time steps for element
  uint d = ed[e0];

  // Compute number of such dependencies for each nodal point
  const uint end = ( e0 < (ne - 1) ? ed[e0 + 1] : nd );
  const uint ndep = (end - d) / nn;
  dolfin_assert(ndep * nn == (end - d));

  // Evaluate the right-hand side at all quadrature points
  for (uint m = 0; m < method->qsize(); m++)
  {
    // Compute quadrature point
    const real t = a0 + k0*method->qpoint(m);

    // Update values for components with larger or equal time steps
    for (uint pos = 0; pos < deps.size(); pos++)
    {
      // Get other element
      const uint i1 = deps[pos];
      const int e1 = elast[i1];
      
      // Special case, component has no latest element
      if ( e1 == -1 )
	continue;
      
      // Use fast evaluation for elements in the same sub slab
      const uint s1 = es[e1];
      const uint j1 = e1 * nn;
      if ( s0 == s1 )
      {
	u[i1] = jx[j1 + m];
	continue;
      }

      // Skip components with smaller time steps
      const real b1 = sb[s1];
      if ( b1 < (a0 + DOLFIN_EPS) )
       	continue;
      
      // Interpolate value from larger element
      const real a1 = sa[s1];
      const real k1 = b1 - a1;
      const real tau = (t - a1) / k1;
      u[i1] = method->ueval(0.0, jx + j1, tau);
    }

    // Update values for components with smaller time steps
    for (uint dep = 0; dep < ndep; dep++)
    {
      // Get element
      const int e1 = de[d++];
      dolfin_assert(e1 != -1);

      // Interpolate value from smaller element
      const uint i1 = ei[e1];
      const uint s1 = es[e1];
      const real a1 = sa[s1];
      const real b1 = sb[s1];
      const real k1 = b1 - a1;
      const real tau = (t - a1) / k1;
      const uint j1 = e1 * nn;
      u[i1] = method->ueval(0.0, jx + j1, tau);
    }
    
    // Evaluate right-hand side
    f[m] = ode.f(u, t, i0);
  }
}
//-----------------------------------------------------------------------------
TimeSlabSolver* MultiAdaptiveTimeSlab::chooseSolver()
{
  bool implicit = ode.get("ODE implicit");
  std::string solver = ode.get("ODE nonlinear solver");

  if ( implicit )
    error("Multi-adaptive solver cannot solver implicit ODEs. Use cG(q) or dG(q) instead.");

  if ( solver == "fixed-point" )
  {
    message("Using multi-adaptive fixed-point solver.");
    return new MultiAdaptiveFixedPointSolver(*this);
  }
  else if ( solver == "newton" )
  {
    message("Using multi-adaptive Newton solver.");
    return new MultiAdaptiveNewtonSolver(*this);
  }
  else if ( solver == "default" )
  {
    message("Using multi-adaptive fixed-point solver (default for mc/dG(q)).");
    return new MultiAdaptiveFixedPointSolver(*this);
  }
  else
  {
    error("Uknown solver type: %s.", solver.c_str());
  }

  return 0;
}
//------------------------------------------------------------------------