TimeGLMGEE.c

Go to the documentation of this file.

#include <string.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <simulation_object.h>
#include <timeintegration.h>

int TimeGLMGEE(void *ts )
{
  TimeIntegration* TS = (TimeIntegration*) ts;
  SimulationObject* sim = (SimulationObject*) TS->simulation;
  GLMGEEParameters* params = (GLMGEEParameters*) sim[0].solver.msti;
  int ns, stage, i, j, nsims = TS->nsims;
  _DECLARE_IERR_;

  int    s  = params->nstages;
  int    r  = params->r;
  double dt = TS->dt;
  double  *A =params->A, 
          *B =params->B, 
          *C=params->C, 
          *D=params->D, 
          *c=params->c,
          **U = TS->U, 
          **Udot = TS->Udot, 
          **Uaux = &TS->U[r];

  /* Calculate stage values */
  for (j=0; j<s; j++) {

    double stagetime = TS->waqt + c[j]*dt;

    for (ns = 0; ns < nsims; ns++) {
      _ArrayScaleCopy1D_( sim[ns].solver.u,
                          C[j*r+0],
                          (U[0] + TS->u_offsets[ns]),
                          (TS->u_sizes[ns]) );
    }

    for (i=1;i<r;i++) _ArrayAXPY_(Uaux[i-1], C[j*r+i]   , U[0], TS->u_size_total);
    for (i=0;i<j;i++) _ArrayAXPY_(Udot[i]  , dt*A[j*s+i], U[0], TS->u_size_total); 

    for (ns = 0; ns < nsims; ns++) {
      if (sim[ns].solver.PreStage) { 
        fprintf(stderr,"Call to solver->PreStage() commented out in TimeGLMGEE()!\n");
        return 1;
//        IERR sim[ns].solver.PreStage( stage,
//                                      (U),
//                                      &(sim[ns].solver),
//                                      &(sim[ns].mpi),
//                                      stagetime ); CHECKERR(ierr); 
      }
    }

    for (ns = 0; ns < nsims; ns++) {
      IERR TS->RHSFunction( (Udot[j] + TS->u_offsets[ns]),
                            (U[0] + TS->u_offsets[ns]),
                            &(sim[ns].solver),
                            &(sim[ns].mpi),
                            stagetime);
    }

    for (ns = 0; ns < nsims; ns++) {
      if (sim[ns].solver.PostStage) { 
        IERR sim[ns].solver.PostStage(  (U[j] + TS->u_offsets[ns]),
                                        &(sim[ns].solver),
                                        &(sim[ns].mpi),
                                        stagetime); CHECKERR(ierr); 
      }
    }

    _ArraySetValue_(TS->BoundaryFlux[j], TS->bf_size_total, 0.0);
    for (ns = 0; ns < nsims; ns++) {
      _ArrayCopy1D_(  sim[ns].solver.StageBoundaryIntegral,
                      (TS->BoundaryFlux[j] + TS->bf_offsets[ns]),
                      TS->bf_sizes[ns] );
    }
  }

  /* Step completion */
  for (j=0; j<r; j++) {

    for (ns = 0; ns < nsims; ns++) {

      _ArrayScaleCopy1D_( sim[ns].solver.u,
                          D[j*r+0],
                          (U[j]+TS->u_offsets[ns]),
                          (TS->u_sizes[ns]) );

    }

    for (i=1; i<r; i++) _ArrayAXPY_(Uaux[i-1], D[j*r+i]   , U[j], TS->u_size_total);
    for (i=0; i<s; i++) _ArrayAXPY_(Udot[i]  , dt*B[j*s+i], U[j], TS->u_size_total);

  }

  for (ns = 0; ns < nsims; ns++) {
    for (i=0; i<s; i++) {
      _ArrayAXPY_(  (TS->BoundaryFlux[i] + TS->bf_offsets[ns]),
                    dt*B[0*s+i],
                    sim[ns].solver.StepBoundaryIntegral,
                    TS->bf_sizes[ns] );
    }

  }
    
  for (ns = 0; ns < nsims; ns++) {
    _ArrayCopy1D_(  (U[0] + TS->u_offsets[ns]),
                    sim[ns].solver.u,
                    TS->u_sizes[ns] );
  }
  for (i=1; i<r; i++) {
    _ArrayCopy1D_(U[i],Uaux[i-1],TS->u_size_total);
  }

  return(0);
}