diff options
Diffstat (limited to 'src/RK87.c')
| -rw-r--r-- | src/RK87.c | 172 |
1 files changed, 72 insertions, 100 deletions
@@ -14,6 +14,7 @@ #include "cctk_Parameters.h" #include "ExternalVariables.h" +#include "Operators.h" static const char *rcsid = "$Header$"; @@ -70,17 +71,9 @@ void MoL_RK87Add(CCTK_ARGUMENTS) CCTK_INT arraydim; static CCTK_INT scratchspace_firstindex = -99; - CCTK_INT index, var, scratchstep; + CCTK_INT var, scratchstep; CCTK_INT totalsize; - CCTK_REAL * restrict UpdateVar; - CCTK_REAL const * restrict RHSVar; - CCTK_REAL * restrict ScratchVar; - CCTK_REAL * restrict ErrorVar; - CCTK_REAL const * restrict OldVar; - - CCTK_REAL beta, gamma, gamma_error; - static const CCTK_REAL beta_array[12][12] = { { 1.0/18.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, { 1.0/48.0, 1.0/16.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, @@ -131,7 +124,7 @@ void MoL_RK87Add(CCTK_ARGUMENTS) totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } if (scratchspace_firstindex == -99) @@ -145,122 +138,101 @@ void MoL_RK87Add(CCTK_ARGUMENTS) for (var = 0; var < MoLNumEvolvedVariables; var++) { - const CCTK_REAL tmp = (*Original_Delta_Time) / cctkGH->cctk_timefac; - - UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedVariableIndex[var]); - RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, - RHSVariableIndex[var]); - ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - CCTK_FirstVarIndex("MOL::SCRATCHSPACE") - + var - + MoL_Num_Evolved_Vars * - (MoL_Intermediate_Steps - - (*MoL_Intermediate_Step))); -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - ScratchVar[index] = tmp * RHSVar[index]; - } - + int const step = MoL_Intermediate_Steps - (*MoL_Intermediate_Step); + int const scratchvarindex = + scratchspace_firstindex + MoL_Num_Evolved_Vars * step + var; + + int const nsrcs = 1; + CCTK_INT const srcs[] = {RHSVariableIndex[var]}; + CCTK_INT const tls[] = {0}; + CCTK_REAL const facts[] = {(*Original_Delta_Time) / cctkGH->cctk_timefac}; + MoL_LinearCombination(cctkGH, + scratchvarindex, 1.0, + srcs, tls, facts, nsrcs); } - for (var = 0; var < MoLNumEvolvedVariables; var++) { - OldVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 1, - EvolvedVariableIndex[var]); - UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedVariableIndex[var]); - RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, - CCTK_FirstVarIndex("MOL::SCRATCHSPACE") - + var - + MoL_Num_Evolved_Vars * - (MoL_Intermediate_Steps - - (*MoL_Intermediate_Step))); - ErrorVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - CCTK_FirstVarIndex("MOL::ERRORESTIMATE") - + var); + int const error_firstindex = CCTK_FirstVarIndex("MOL::ERRORESTIMATE"); + int const step = MoL_Intermediate_Steps - (*MoL_Intermediate_Step); if (*MoL_Intermediate_Step - 1) - { - -#pragma omp parallel for - for (index = 0; index < totalsize; index++) + { + int const nsrcs = step + 2; + CCTK_INT srcs[nsrcs]; + CCTK_INT tls[nsrcs]; + CCTK_REAL facts[nsrcs]; + srcs[0] = EvolvedVariableIndex[var]; + tls[0] = 1; + facts[0] = 1.0; + + for (scratchstep = 0; scratchstep < step + 1; scratchstep++) { - UpdateVar[index] = OldVar[index]; + int const scratchvarindex = + scratchspace_firstindex + MoL_Num_Evolved_Vars * scratchstep + var; + srcs [scratchstep+1] = scratchvarindex; + tls [scratchstep+1] = 0; + facts[scratchstep+1] = beta_array[step][scratchstep]; } - for (scratchstep = 0; - scratchstep < MoL_Intermediate_Steps - (*MoL_Intermediate_Step) + 1; - scratchstep++) + MoL_LinearCombination(cctkGH, + EvolvedVariableIndex[var], 0.0, + srcs, tls, facts, nsrcs); + } + else + { { + int const nsrcs = 14; + CCTK_INT srcs[nsrcs]; + CCTK_INT tls[nsrcs]; + CCTK_REAL facts[nsrcs]; + srcs[0] = EvolvedVariableIndex[var]; + tls[0] = 1; + facts[0] = 1.0; - ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - CCTK_FirstVarIndex("MOL::SCRATCHSPACE") - + var - + MoL_Num_Evolved_Vars * scratchstep); - - beta = beta_array[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)][scratchstep]; - - if ( (beta > MoL_Tiny)||(beta < -MoL_Tiny) ) + for (scratchstep = 0; scratchstep < 13; scratchstep++) { -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateVar[index] += beta * ScratchVar[index]; - } + int const scratchvarindex = + scratchspace_firstindex + MoL_Num_Evolved_Vars * scratchstep + var; + srcs [scratchstep+1] = scratchvarindex; + tls [scratchstep+1] = 0; + facts[scratchstep+1] = gamma_array[scratchstep]; } - - } - - } - else - { -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateVar[index] = OldVar[index]; - ErrorVar[index] = 0; + MoL_LinearCombination(cctkGH, + EvolvedVariableIndex[var], 0.0, + srcs, tls, facts, nsrcs); } - - for (scratchstep = 0; scratchstep < 13; scratchstep++) + { + int const nsrcs = 13; + CCTK_INT srcs[nsrcs]; + CCTK_INT tls[nsrcs]; + CCTK_REAL facts[nsrcs]; - ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - CCTK_FirstVarIndex("MOL::SCRATCHSPACE") - + var - + MoL_Num_Evolved_Vars * scratchstep); - - gamma = gamma_array[scratchstep]; - gamma_error = gamma - gammastar_array[scratchstep]; - - if ( (gamma > MoL_Tiny)||(gamma < -MoL_Tiny) ) + for (scratchstep = 0; scratchstep < 13; scratchstep++) { -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateVar[index] += gamma * ScratchVar[index]; - ErrorVar[index] += gamma_error * ScratchVar[index]; - } + int const scratchvarindex = + scratchspace_firstindex + MoL_Num_Evolved_Vars * scratchstep + var; + srcs [scratchstep] = scratchvarindex; + tls [scratchstep] = 0; + facts[scratchstep] = + gamma_array[scratchstep] - gammastar_array[scratchstep]; } + MoL_LinearCombination(cctkGH, + error_firstindex + var, 0.0, + srcs, tls, facts, nsrcs); } - + } - + } /* Real arrays */ - for (var = 0; var < MoLNumEvolvedArrayVariables; var++) - { - - CCTK_WARN(0, "Peter has been too lazy to write the RK87 routine " - "out for array variables. Better send him an email..."); - - } + CCTK_WARN(0, "Peter has been too lazy to write the RK87 routine " + "out for array variables. Better send him an email..."); return; } |