multiadaptivetimeslab.cpp

利用C

// Copyright (C) 2005-2008 Anders Logg.
// Licensed under the GNU LGPL Version 2.1.
//
// Modified by Garth N. Wells
//
// First added:  2005-01-27
// Last changed: 2008-06-11

#include <string>
#include <algorithm>
#include <dolfin/common/constants.h>
#include <dolfin/log/dolfin_log.h>
#include <dolfin/parameter/parameters.h>
#include "ODE.h"
#include "Dependencies.h"
#include "Method.h"
#include "MultiAdaptiveFixedPointSolver.h"
#include "MultiAdaptiveNewtonSolver.h"
#include "MultiAdaptiveTimeSlab.h"

using namespace dolfin;

//-----------------------------------------------------------------------------
MultiAdaptiveTimeSlab::MultiAdaptiveTimeSlab(ODE& ode) : 
  TimeSlab(ode),
  sa(0), sb(0), ei(0), es(0), ee(0), ed(0), jx(0), de(0),
  ns(0), ne(0), nj(0), nd(0), solver(0), adaptivity(ode, *method), partition(N),
  elast(0), f0(0), emax(0), kmin(0)
{
  // Choose solver
  solver = chooseSolver();

  // Initialize elast
  elast = new int[N];
  for (uint i = 0; i < N; i++)
    elast[i] = -1;

  // Initialize f at left end-point for cG
  if ( method->type() == Method::cG )
    f0 = new real[N];

  // Initialize vector for u
  u.init(N);
  u.zero();

  // Initialize transpose of dependencies if necessary
  message("Computing transpose (inverse) of dependency pattern.");
  if ( ode.dependencies.sparse() && !ode.transpose.sparse() )
    ode.transpose.transp(ode.dependencies);
}
//-----------------------------------------------------------------------------
MultiAdaptiveTimeSlab::~MultiAdaptiveTimeSlab()
{
  if ( sa ) delete [] sa;
  if ( sb ) delete [] sb;
  if ( ei ) delete [] ei;
  if ( es ) delete [] es;
  if ( ee ) delete [] ee;
  if ( ed ) delete [] ed;
  if ( jx ) delete [] jx;
  if ( de ) delete [] de;

  if ( solver ) delete solver;

  if ( elast ) delete [] elast;
  if ( f0 ) delete [] f0;
}
//-----------------------------------------------------------------------------
real MultiAdaptiveTimeSlab::build(real a, real b)
{
  //cout << "Multi-adaptive time slab: building between "
  //     << a << " and " << b << endl;
  
  // Allocate data
  allocData(a, b);

  // Reset elast
  for (uint i = 0; i < N; i++)
    elast[i] = -1;

  // Create time slab recursively
  kmin = ode.endtime();
  b = createTimeSlab(a, b, 0);

  //cout << "de = "; Alloc::disp(de, nd);
  
  // Save start and end time
  _a = a;
  _b = b;

  //cout << "Multi-adaptive time slab: finished building between "
  //     << a << " and " << b << ": K = " << b - a << ", nj = " << nj << endl;

  // Update at t = 0.0
  if ( a < DOLFIN_EPS )
    ode.update(u0, a, false);

  return b;
}
//-----------------------------------------------------------------------------
bool MultiAdaptiveTimeSlab::solve()
{
  //message("Solving time slab system on [%f, %f].", _a, _b);

  // Copy u0 to u. This happens automatically in feval if user has set
  // dependencies correctly, but you never know...
  for (unsigned int i = 0; i < N; i++)
    u[i] = u0[i];

  // Compute f at left end-point for cG
  if ( method->type() == Method::cG )
  {
    for (uint i = 0; i < N; i++)
      f0[i] = ode.f(u0, _a, i);
  }

  // Solve system
  return solver->solve();

  //for (uint i = 0; i < N; i++)
  // {
  //  real endval = jx[elast[i] * method->nsize() + method->nsize() - 1];
  //  message("i = %d: u = %.16e", i, endval);
  // }
}
//-----------------------------------------------------------------------------
bool MultiAdaptiveTimeSlab::check(bool first)
{
  // Compute new time steps
  adaptivity.update(*this, _b, first);

  // Check if current solution can be accepted
  return adaptivity.accept();
}
//-----------------------------------------------------------------------------
bool MultiAdaptiveTimeSlab::shift(bool end)
{
  // Cover end time
  coverTime(_b);

  // Update the solution vector at the end time for each component
  for (uint i = 0; i < N; i++)
  {
    // Get last element of component
    const int e = elast[i];
    dolfin_assert(e != -1);
    dolfin_assert(sb[es[e]] == _b);
    
    // Get end-time value of component
    const int j = e * method->nsize();
    u[i] = jx[j + method->nsize() - 1];
  }

  // Write solution at final time if we should
  if ( save_final && end )
    write(u);

  // Let user update ODE
  if ( !ode.update(u, _b, end) )
    return false;

  // Set initial value to end-time value
  for (uint i = 0; i < N; i++)
    u0[i] = u[i];

  return true;
}
//-----------------------------------------------------------------------------
void MultiAdaptiveTimeSlab::reset()
{
  // Iterate over all elements
  uint j = 0;
  for (uint e = 0; e < ne; e++)
  {
    // Get component index
    const uint i = ei[e];

    // Iterate over degrees of freedom on element
    for (uint n = 0; n < method->nsize(); n++)
      jx[j + n] = u0[i];

    // Step to next element
    j += method->nsize();
  }
}
//-----------------------------------------------------------------------------
void MultiAdaptiveTimeSlab::sample(real t)
{
  // Cover the given time
  coverTime(t);

  //cout << "t = " << t << " elast: ";
  //for (uint i = 0; i < N; i++)
  //  cout << elast[i] << " ";
  //cout << endl;
}
//-----------------------------------------------------------------------------
real MultiAdaptiveTimeSlab::usample(uint i, real t)
{
  // Get element
  const int e = elast[i];
  dolfin_assert(e != -1);

  // Get element data
  const uint s = es[e];
  const uint j = e * method->nsize();
  const real a = sa[s];
  const real b = sb[s];
  const real k = b - a;

  // Get initial value for element (only necessary for cG)
  const int ep = ee[e];
  const uint jp = ep * method->nsize();
  const real x0 = ( ep != -1 ? jx[jp + method->nsize() - 1] : u0[i] );
  
  // Evaluate solution
  const real tau = (t - a) / k;
  const real value = method->ueval(x0, jx + j, tau);

  return value;
}
//-----------------------------------------------------------------------------
real MultiAdaptiveTimeSlab::ksample(uint i, real t)
{
  // Get element
  const int e = elast[i];
  dolfin_assert(e != -1);

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

  // Compute time step
  const real k = b - a;

  return k;
}
//-----------------------------------------------------------------------------
real MultiAdaptiveTimeSlab::rsample(uint i, real t)
{
  /*
  // Note that the residual is always sampled at the end-time
  
  // Cover end time
  coverTime(_b);
  
  // Update the solution vector at the end time for each dependent component
  
  // Get list of dependencies for component
  const Array<uint>& deps = ode.dependencies[i];
  
  // Iterate over dependencies
  for (uint pos = 0; pos < deps.size(); pos++)
  {
    // Get last element of component
    const int e = elast[pos];
    dolfin_assert(e != -1);
    dolfin_assert(sb[es[e]] == _b);
    
    // Get end-time value of component
    const int j = e * method->nsize();
    u[pos] = jx[j + method->nsize() - 1];
  }
  
  // Compute residual
  
  // Get last element of component
  const int e = elast[i];
  dolfin_assert(e != -1);
  
  // Get element data
  const uint s = es[e];
  const uint j = e * method->nsize();
  const real a = sa[s];
  const real b = sb[s];
  const real k = b - a;
  
  // Get initial value for element (only necessary for cG)
  const int ep = ee[e];
  const uint jp = ep * method->nsize();
  const real x0 = ( ep != -1 ? jx[jp + method->nsize() - 1] : u0(i) );
  
  // Evaluate right-hand side at end-point (u is already updated)
  const real f = ode.f(u, b, i);
  
  // Compute residual
  const real r = method->residual(x0, jx + j, f, k);
  */

  // Just return previously computed maximum in time slab for component
  return adaptivity.residual(i);
}
//-----------------------------------------------------------------------------
void MultiAdaptiveTimeSlab::disp() const
{
  cout << "--------------------------------------------------------" << endl;

  message("s: size = %d alloc = %d", ns, size_s.size);
  message("e: size = %d alloc = %d", ne, size_e.size);
  message("j: size = %d alloc = %d", nj, size_j.size);
  message("d: size = %d alloc = %d", nd, size_d.size);

  cout << "sa = "; Alloc::disp(sa, ns);
  cout << "sb = "; Alloc::disp(sb, ns);
 
  cout << endl;

  cout << "ei = "; Alloc::disp(ei, ne);
  cout << "es = "; Alloc::disp(es, ne);  
  cout << "ee = "; Alloc::disp(ee, ne);
  cout << "ed = "; Alloc::disp(ed, ne);

  cout << endl;

  cout << "jx = "; Alloc::disp(jx, nj);

  cout << endl;

  cout << "de = "; Alloc::disp(de, nd);
}
//-----------------------------------------------------------------------------
void MultiAdaptiveTimeSlab::allocData(real a, real b)
{ 
  // Use u to keep track of the latest time value for each component here
  for (uint i = 0; i < N; i++)
    u[i] = a;
  
  // Recursively compute data size
  ns = ne = nj = nd = 0;
  computeDataSize(a, b, 0);

  // Allocate data
  alloc_s(ns);
  alloc_e(ne);
  alloc_j(nj);
  alloc_d(nd);

  // Reset mapping de
  for (uint d = 0; d < nd; d++)
    de[d] = -1;
}
//-----------------------------------------------------------------------------
real MultiAdaptiveTimeSlab::createTimeSlab(real a, real b, uint offset)
{
  // Compute end time of this sub slab
  uint end = 0;
  b = computeEndTime(a, b, offset, end);

  //cout << "Creating sub slab between a = " << a << " and b = " << b << endl;

  // Create sub slab
  create_s(a, b, offset, end);

  // Recursively create sub slabs for components with small time steps
  real t = a;
  while ( t < b && end < partition.size() )
    t = createTimeSlab(t, b, end);
  
  return b;
}
//-----------------------------------------------------------------------------
void MultiAdaptiveTimeSlab::create_s(real a0, real b0, uint offset, uint end)
{
  dolfin_assert(size_s.next < size_s.size);
  
  // Get next available position
  uint pos = size_s.next++;

  // Create new sub slab
  sa[pos] = a0;
  sb[pos] = b0;

  // Create elements for sub slab
  for (uint n = offset; n < end; n++)
    create_e(partition.index(n), pos, a0, b0);

  // Create mapping ed
  for (uint n = offset; n < end; n++)
  {
    const uint index = partition.index(n);
    const int element = elast[index];
    dolfin_assert(element != -1);

    // Count number of dependencies from element
    size_d.next += countDependencies(index, b0);

    // Update mapping ed
    if ( element == 0 )
      ed[element] = 0;
    if ( element < static_cast<int>(ne - 1) )
      ed[element + 1] = size_d.next;
  }
}
//-----------------------------------------------------------------------------
void MultiAdaptiveTimeSlab::create_e(uint index, uint subslab, real a, real b)
{
  dolfin_assert(size_e.next < size_e.size);
  
  // Get next available position
  uint pos = size_e.next++;

  //message("  Creating element e = %d for i = %d at [%f, %f]", pos, index, a, b);
  
  //if ( index == 145 )
  //  cout << "Modified: " << b - a << endl << endl;

  // Create new element
  ei[pos] = index;
  es[pos] = subslab;
  ee[pos] = elast[index];

  // Create dofs for element
  create_j(index);

  // Create dependencies to element
  create_d(index, pos, subslab, a, b);

  // Update elast for component
  elast[index] = pos;
}
//-----------------------------------------------------------------------------
void MultiAdaptiveTimeSlab::create_j(uint index)
{
  dolfin_assert((size_j.next + method->nsize() - 1) < size_j.size);

  // Get next available position
  uint pos = size_j.next;
  size_j.next += method->nsize();

  // Create dofs
  for (uint n = 0; n < method->nsize(); n++)
    jx[pos + n] = u0[index];
}
//-----------------------------------------------------------------------------
void MultiAdaptiveTimeSlab::create_d(uint i0, uint e0, uint s0, real a0, real b0)
{
  // Add dependencies to elements that depend on the given element if the
  // depending elements use larger time steps
  
  //message("Checking dependencies to element %d (component %d)", element, index);

  // Get list of components depending on current component
  const Array<uint>& deps = ode.transpose[i0];
  
  // Iterate over dependencies
  for (uint pos = 0; pos < deps.size(); pos++)
  {
    // Get component index of other component
    const uint i1 = deps[pos];
        
    // Get other element
    const int e1 = elast[i1];
    
    //cout << "  Other element: " << e1 << endl;
    
    // Skip elements which have not been created (use smaller time steps)
    if ( e1 == -1 )
      continue;
    
    // Get data of other element
    const uint s1 = es[e1];
    const real a1 = sa[s1];
    const real b1 = sb[s1];
    const real k1 = b1 - a1;
    
    // Only add dependencies from components with larger time steps
    if ( !within(a0, b0, a1, b1) || s0 == s1 )
      continue;
    
    //message("  Checking element %d (component %d)", e1, i1);
    
    // Iterate over dofs for element
    for (uint n = 0; n < method->nsize(); n++)
    {
      //const uint j = j1 + n;
      const real t = a1 + k1*method->npoint(n);
      
      //message("    Checking dof at t = %f", t);
      
      // Check if dof is contained in the current element
      if ( within(t, a0, b0) )
      {
	// Search for an empty position
	bool found = false;
	
	//cout << "    --- Creating dependency to element = " << element << endl;
	//cout << "    --- Starting at ed = " << ed[e1] << endl;
	//cout << "    de = "; Alloc::disp(ed, ne);
	//cout << "    nd = "; Alloc::disp(de, nd);	

	for (uint d = ed[e1]; d < ed[e1 + 1]; d++)
	{
	  if ( de[d] == -1 )
	  {
	    de[d] = e0;
	    found = true;
	    break;
	  }
	}
	dolfin_assert(found);
      }
    }
  }
}
//-----------------------------------------------------------------------------
void MultiAdaptiveTimeSlab::alloc_s(uint newsize)
{
  size_s.next = 0;

  if ( newsize <= size_s.size ) return;

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

  Alloc::realloc(&sa, size_s.size, newsize);
  Alloc::realloc(&sb, size_s.size, newsize);

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

  if ( newsize <= size_e.size ) return;

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

  Alloc::realloc(&ei, size_e.size, newsize);