Sampling based on constant weights

New ideas for code development

Sampling based on constant weights

Postby jrpe » Tue Mar 13, 2018 12:48 pm

Hello All,
I have recently tried to modify the source normalization routine for critical transport in order to keep the source batch constant for all cycles by sampling random particles from the source. I found that NormalizeCritSrc() empty the bank with function FromBank(id)and loads the particles to source by AddItem(DATA_PART_PTR_SOURCE, ptr). As a quick implementation I just transfer the source back to the bank, sample duplicates to source and flush bank at the end.

To avoid weight resizing, I have replaced line 269
Code: Select all
w0 = RDB[ptr + PARTICLE_WGT]*((double)nbatch)/((double)neig)/wgtg[i];

by just:
Code: Select all
w0 = RDB[ptr + PARTICLE_WGT];


The following code have been added at line 343 (after CalculateEntropies();) to implement the sampling routine:
Code: Select all

  long PtrStore [  (long) RDB[DATA_CYCLE_BATCH_SIZE]  ];
  long IdStore [  (long) RDB[DATA_CYCLE_BATCH_SIZE]  ];
  double cWgt [ (long) RDB[DATA_CYCLE_BATCH_SIZE]  ];
  cWgt[0] = 0.0;
  long nn, cntr;
  long ngg;
  long tempptr;

   /*Step 1: get all ptr's to PtrStore while empty the source storage, and calculate cumulative weight*/
   cntr = 0;
    for (id = 0; id < (long)RDB[DATA_OMP_MAX_THREADS]; id++)
    {
      while ( ( ptr = FromSrc(id)  )  > VALID_PTR) {
         PtrStore[cntr] = ptr;
         IdStore[cntr] = id;
         if (cntr  >0)
            cWgt[cntr] = cWgt[cntr-1] + RDB[ptr + PARTICLE_WGT];

            ToBank(ptr, 0);
            /* AddItem(OMPPtr(DATA_PART_PTR_BANK, 0), ptr); */


         cntr += 1;

      }
    }

    /* step 2: duplicate particles from bank to source by sampling*/

    for ( nn = 0; nbatch+1 > nn; nn++)
    {
      double Eps = RandF(0)*wtotal;

      for (ngg = 0; RDB[DATA_CYCLE_BATCH_SIZE] > ngg; ngg++)
      {
         if (cWgt[ngg] > Eps)
            {
                ptr = -1;
                    while ((ptr = ViewBank(0, ptr)) > VALID_PTR)
                    {
                        if ( ptr == PtrStore[ngg] && ptr != tempptr )
                        {
                            ptr = DuplicateParticle(ptr, 0);
                            AddItem(DATA_PART_PTR_SOURCE, ptr);
                            break;
                        }

                    }

                break;
            }
      }
    }

    /* step 3: Flush stack/bank */
    for (ngg = 0; RDB[DATA_CYCLE_BATCH_SIZE] > ngg; ngg++)
    {
        ptr = PtrStore[ngg];
        RemoveItem(ptr);
    }

the "ViewBank(0, ptr)" function is just a copy of FromBank that only get the pointers without removing them from bank.
viewbank.c:
Code: Select all
#include "header.h"
#include "locations.h"

#define FUNCTION_NAME "ViewBank:"

/*****************************************************************************/

long ViewBank(long id, long ptr)
{

  /* Check actual thread number */

  if (OMP_THREAD_NUM != 0)
    Die(FUNCTION_NAME, "Called from an OpenMP parallel loop");

  /* Get pointer if ptr == -1*/
if (ptr == -1 )
{
  ptr = (long)RDB[OMPPtr(DATA_PART_PTR_BANK, id)];
  CheckPointer(FUNCTION_NAME, "(ptr)", DATA_ARRAY, ptr);

  /* Get pointer to first item */

  ptr = FirstItem(ptr);
  CheckPointer(FUNCTION_NAME, "(ptr)", DATA_ARRAY, ptr);
}
else if (ptr > VALID_PTR)
{
    ptr = ptr;
}
else
    Die(FUNCTION_NAME, "Error in ptr format");


if ((ptr = NextItem(ptr)) < VALID_PTR)
    {
      /* Bank is empty, return null */

      return -1;
    }

return ptr;
}


The code passes the total weight control, but fails when the TransportCycle() reaches while(FromSrc(id) > VALID_PTR) {Tracking(id);} in line 1911 of transportcycle.c and return an error requesting higher nbuf. Is there some conflicts with memory management or pointer locations that my code causes? Is there perhaps a better method than using the bank storage?
jrpe
 
Posts: 1
Joined: Wed Jan 31, 2018 1:31 pm

Return to Development

Who is online

Users browsing this forum: No registered users and 2 guests