## Sampling based on constant weights

New ideas for code development

### Sampling based on constant weights

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