Skip to content

Commit

Permalink
Initial sink term (ORNL-Fusion#112)
Browse files Browse the repository at this point in the history
  • Loading branch information
Varun Shah authored and Varun Shah committed Jul 6, 2023
1 parent db865ea commit 74ad1ef
Show file tree
Hide file tree
Showing 3 changed files with 22 additions and 2 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,7 @@ class CustomFitFluxHandler : public FluxHandler
xGrid = grid;

// Read the parameter file
fs::ifstream paramFile(profileFilePath);
std::ifstream paramFile(profileFilePath.c_str());

// Gets the process ID
int procId;
Expand Down
17 changes: 16 additions & 1 deletion xolotl/core/include/xolotl/core/network/impl/PSIReaction.tpp
Original file line number Diff line number Diff line change
Expand Up @@ -247,6 +247,8 @@ PSISinkReaction<TSpeciesEnum>::computeRate(IndexType gridIndex, double time)
{
auto cl = this->_clusterData->getCluster(this->_reactant);
double dc = cl.getDiffusionCoefficient(gridIndex);

double grainSize = 50000.0; // 50 um

auto rate = this->_clusterData->extraData.leftSideRates;
double s_m = 0.0;
Expand All @@ -260,8 +262,21 @@ PSISinkReaction<TSpeciesEnum>::computeRate(IndexType gridIndex, double time)
}
s_m = rate(i) / dc;
}

if (util::equal(s_m,0)){
return 0;
}

double s_sk = 0; // Define the gb sink strength
double tan = 0; // Define the tanh term
double cot = 0; // Defin the coth term
tan = tanh(sqrt(s_m)*grainSize/2);
double term1 = (sqrt(s_m)*grainSize)/2*(1/tan); // Repetitive term
// Calculate the sink term
s_sk = s_m*(term1 - 1)*(1/(1+(s_m*(grainSize*grainSize)/12) - term1));
//std::cout<<tan<<" "<<term1<<" "<<s_sk<<std::endl;

double strength = this->getSinkBias() * this->getSinkStrength() * dc;
double strength = s_sk*dc;

return strength;
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -781,6 +781,11 @@ PSIReactionGenerator<TSpeciesEnum>::addSinks(IndexType i, TTag tag) const
if (lo[Species::V] == 1)
this->addSinkReaction(tag, {i, NetworkType::invalidIndex()});
}

//He
if (clReg.isSimplex() && lo.isOnAxis(Species::He)) {
this->addSinkReaction(tag, {i, NetworkType::invalidIndex()});
}
}

template <typename TSpeciesEnum>
Expand Down

0 comments on commit 74ad1ef

Please sign in to comment.