13 #include "parstream.H"
14 #include "ParmParse.H"
15 #include "LayoutIterator.H"
18 #include "LoadBalance.H"
19 #include "BoxIterator.H"
20 #include "computeSum.H"
24 #include "AMRLevelPluto.H"
26 #include "NamespaceHeader.H"
29 AMRLevelPluto::AMRLevelPluto()
33 pout() <<
"AMRLevelPluto default constructor" << endl;
36 m_patchPlutoFactory = NULL;
38 m_paramsDefined =
false;
42 AMRLevelPluto::~AMRLevelPluto()
46 pout() <<
"AMRLevelPluto destructor" << endl;
50 if (m_patchPluto != NULL)
55 if (m_patchPlutoFactory != NULL)
57 delete m_patchPlutoFactory;
60 m_paramsDefined =
false;
63 void AMRLevelPluto::defineParams(
const Real& a_cfl,
64 const Real& a_domainLength,
65 const int& a_verbosity,
66 const Real& a_refineThresh,
67 const int& a_tagBufferSize,
68 const Real& a_initialDtMultiplier,
69 const PatchPluto*
const a_patchPluto)
75 m_domainLength = a_domainLength;
78 verbosity(a_verbosity);
81 m_refineThresh = a_refineThresh;
84 m_tagBufferSize = a_tagBufferSize;
87 initialDtMultiplier(a_initialDtMultiplier);
92 if (m_patchPlutoFactory != NULL)
94 delete m_patchPlutoFactory;
97 m_patchPlutoFactory = a_patchPluto->new_patchPluto();
100 if (m_patchPlutoFactory != NULL)
105 m_patchPluto = a_patchPluto->new_patchPluto();
107 m_paramsDefined =
true;
111 void AMRLevelPluto::define(AMRLevel* a_coarserLevelPtr,
112 const Box& a_problemDomain,
116 ProblemDomain physdomain(a_problemDomain);
118 MayDay::Error(
"AMRLevelPluto::define -\n\tShould never be called with a Box for a problem domain");
122 void AMRLevelPluto::define(AMRLevel* a_coarserLevelPtr,
123 const ProblemDomain& a_problemDomain,
127 if (s_verbosity >= 3)
129 pout() <<
"AMRLevelPluto::define " << a_level << endl;
133 AMRLevel::define(a_coarserLevelPtr,
139 if (a_coarserLevelPtr != NULL)
141 AMRLevelPluto* amrGodPtr =
dynamic_cast<AMRLevelPluto*
>(a_coarserLevelPtr);
143 if (amrGodPtr != NULL)
145 m_cfl = amrGodPtr->m_cfl;
146 m_domainLength = amrGodPtr->m_domainLength;
147 m_refineThresh = amrGodPtr->m_refineThresh;
148 m_tagBufferSize = amrGodPtr->m_tagBufferSize;
152 MayDay::Error(
"AMRLevelPluto::define: a_coarserLevelPtr is not castable to AMRLevelPluto*");
157 m_dx = m_domainLength / (a_problemDomain.domainBox().bigEnd(0)-a_problemDomain.domainBox().smallEnd(0)+1.);
163 CH_assert(m_patchPluto != NULL);
164 CH_assert(isDefined());
165 m_patchPluto->define(m_problem_domain,m_dx,m_level,m_numGhost);
168 m_numStates = m_patchPluto->numConserved();
169 m_ConsStateNames = m_patchPluto->ConsStateNames();
170 m_PrimStateNames = m_patchPluto->PrimStateNames();
174 Real AMRLevelPluto::advance()
176 CH_assert(allDefined());
178 if (s_verbosity >= 3)
180 pout() <<
"AMRLevelPluto::advance level " << m_level <<
" to time " << m_time + m_dt << endl;
185 DataIterator dit = m_UNew.dataIterator();
186 for(;dit.ok(); ++dit){
187 m_UOld[dit()].copy(m_UNew[dit()]);
196 LevelFluxRegister dummyFR;
199 const LevelData<FArrayBox> dummyData;
202 LevelFluxRegister* coarserFR = &dummyFR;
203 LevelFluxRegister* finerFR = &dummyFR;
205 const LevelData<FArrayBox>* coarserDataOld = &dummyData;
206 const LevelData<FArrayBox>* coarserDataNew = &dummyData;
208 Real tCoarserOld = 0.0;
209 Real tCoarserNew = 0.0;
214 AMRLevelPluto* coarserPtr = getCoarserLevel();
218 coarserFR = &coarserPtr->m_fluxRegister;
220 coarserDataOld = &coarserPtr->m_UOld;
221 coarserDataNew = &coarserPtr->m_UNew;
223 tCoarserNew = coarserPtr->m_time;
224 tCoarserOld = tCoarserNew - coarserPtr->m_dt;
232 finerFR = &m_fluxRegister;
236 LevelData<FArrayBox> flux[SpaceDim];
241 if ( m_time < tCoarserOld &&
242 fabs(m_time - tCoarserOld) < 1.e-12*tCoarserOld){
243 m_time = tCoarserOld;
249 newDt = m_levelPluto.step(m_UNew,
262 #if (TIME_STEPPING == RK2)
266 DtCool = m_levelPluto.step(m_UNew,
279 newDt = Min(newDt,DtCool);
284 Real returnDt = m_cfl * newDt;
291 Real AMRLevelPluto::getDlMin()
294 Real dlMin = m_levelPluto.getDlMin();
301 void AMRLevelPluto::postTimeStep()
303 CH_assert(allDefined());
306 static Real orig_integral = 0.0;
307 static Real last_integral = 0.0;
308 static bool first =
true;
310 if (s_verbosity >= 3)
312 pout() <<
"AMRLevelPluto::postTimeStep " << m_level << endl;
318 Real scale = -1.0/m_dx;
319 m_fluxRegister.reflux(m_UNew,scale);
322 AMRLevelPluto* amrGodFinerPtr = getFinerLevel();
324 amrGodFinerPtr->m_coarseAverage.averageToCoarse(m_UNew,
325 amrGodFinerPtr->m_UNew);
328 if (s_verbosity >= 2 && m_level == 0)
332 pout() <<
"AMRLevelPluto::postTimeStep:" << endl;
333 pout() <<
" Sums:" << endl;
334 for (
int comp = 0; comp < m_numStates; comp++)
336 Interval curComp(comp,comp);
337 Real integral = computeSum(m_UNew,NULL,nRefFine,m_dx,curComp);
339 pout() <<
" " << setw(23)
341 << setiosflags(ios::showpoint)
342 << setiosflags(ios::scientific)
344 <<
" --- " << m_ConsStateNames[comp];
346 if (comp == 0 && !first) {
347 pout() <<
" (" << setw(23)
349 << setiosflags(ios::showpoint)
350 << setiosflags(ios::scientific)
351 << (integral-last_integral)/last_integral
354 << setiosflags(ios::showpoint)
355 << setiosflags(ios::scientific)
356 << (integral-orig_integral)/orig_integral
366 orig_integral = integral;
370 last_integral = integral;
375 if (s_verbosity >= 3)
377 pout() <<
"AMRLevelPluto::postTimeStep " << m_level <<
" finished" << endl;
382 void AMRLevelPluto::tagCells(IntVectSet& a_tags)
384 CH_assert(allDefined());
386 if (s_verbosity >= 3)
388 pout() <<
"AMRLevelPluto::tagCells " << m_level << endl;
392 const DisjointBoxLayout& levelDomain = m_UNew.disjointBoxLayout();
393 IntVectSet localTags;
397 const AMRLevelPluto* amrGodCoarserPtr = getCoarserLevel();
399 PiecewiseLinearFillPluto pwl;
401 pwl.define(levelDomain,
402 amrGodCoarserPtr->m_UNew.disjointBoxLayout(),
404 amrGodCoarserPtr->m_problem_domain,
405 amrGodCoarserPtr->m_ref_ratio,
409 pwl.fillInterp(m_UNew,
410 amrGodCoarserPtr->m_UNew,
411 amrGodCoarserPtr->m_UNew,
417 m_UNew.exchange(Interval(0,m_numStates-1));
419 #if GEOMETRY != CARTESIAN
420 const LevelData<FArrayBox>& dV = m_levelPluto.getdV();
424 DataIterator dit = levelDomain.dataIterator();
425 for (dit.begin(); dit.ok(); ++dit){
426 const Box& b = levelDomain[dit()];
427 FArrayBox gradFab(b,1);
428 FArrayBox& UFab = m_UNew[dit()];
430 #if GEOMETRY != CARTESIAN
431 const FArrayBox& curdV = dV[dit()];
433 const FArrayBox curdV;
436 m_patchPluto->computeRefGradient(gradFab, UFab, curdV, b);
440 for (bit.begin(); bit.ok(); ++bit){
441 const IntVect& iv = bit();
442 if (gradFab(iv) >= m_refineThresh) {
448 localTags.grow(m_tagBufferSize);
452 Box localTagsBox = localTags.minBox();
453 localTagsBox &= m_problem_domain;
454 localTags &= localTagsBox;
460 void AMRLevelPluto::tagCellsInit(IntVectSet& a_tags)
462 CH_assert(allDefined());
464 if (s_verbosity >= 3)
466 pout() <<
"AMRLevelPolytropicGas::tagCellsInit " << m_level << endl;
473 void AMRLevelPluto::regrid(
const Vector<Box>& a_newGrids)
475 CH_assert(allDefined());
477 if (s_verbosity >= 3)
479 pout() <<
"AMRLevelPluto::regrid " << m_level << endl;
483 m_level_grids = a_newGrids;
484 m_grids = loadBalance(a_newGrids);
486 if (s_verbosity >= 4)
490 const DisjointBoxLayout& constGrids = m_grids;
492 pout() <<
"new grids: " << endl;
494 for (LayoutIterator lit = constGrids.layoutIterator(); lit.ok(); ++lit)
496 pout() << constGrids[lit()] << endl;
501 DataIterator dit = m_UNew.dataIterator();
502 for(;dit.ok(); ++dit){
503 m_UOld[dit()].copy(m_UNew[dit()]);
507 IntVect ivGhost = m_numGhost*IntVect::Unit;
508 m_UNew.define(m_grids,m_numStates,ivGhost);
517 AMRLevelPluto* amrGodCoarserPtr = getCoarserLevel();
518 m_fineInterp.interpToFine(m_UNew,amrGodCoarserPtr->m_UNew);
522 m_UOld.copyTo(m_UOld.interval(),
526 m_UOld.define(m_grids,m_numStates,ivGhost);
530 void AMRLevelPluto::initialGrid(
const Vector<Box>& a_newGrids)
532 CH_assert(allDefined());
534 if (s_verbosity >= 3)
536 pout() <<
"AMRLevelPluto::initialGrid " << m_level << endl;
540 m_level_grids = a_newGrids;
541 m_grids = loadBalance(a_newGrids);
543 if (s_verbosity >= 4)
547 const DisjointBoxLayout& constGrids = m_grids;
549 pout() <<
"new grids: " << endl;
550 for (LayoutIterator lit = constGrids.layoutIterator(); lit.ok(); ++lit)
552 pout() << constGrids[lit()] << endl;
557 IntVect ivGhost = m_numGhost*IntVect::Unit;
558 m_UNew.define(m_grids,m_numStates,ivGhost);
559 m_UOld.define(m_grids,m_numStates,ivGhost);
566 void AMRLevelPluto::initialData()
568 if (s_verbosity >= 3)
570 pout() <<
"AMRLevelPluto::initialData " << m_level << endl;
573 m_patchPluto->initialize(m_UNew);
578 void AMRLevelPluto::postInitialize()
580 CH_assert(allDefined());
582 if (s_verbosity >= 3)
584 pout() <<
"AMRLevelPluto::postInitialize " << m_level << endl;
590 AMRLevelPluto* amrGodFinerPtr = getFinerLevel();
592 amrGodFinerPtr->m_coarseAverage.averageToCoarse(m_UNew,
593 amrGodFinerPtr->m_UNew);
601 void AMRLevelPluto::writeCheckpointHeader(HDF5Handle& a_handle)
const
603 CH_assert(allDefined());
605 if (s_verbosity >= 3)
607 pout() <<
"AMRLevelPluto::writeCheckpointHeader" << endl;
611 HDF5HeaderData header;
612 header.m_int[
"num_components"] = m_numStates;
616 for (
int comp = 0; comp < m_numStates; ++comp)
618 sprintf(compStr,
"component_%d",comp);
619 header.m_string[compStr] = m_ConsStateNames[comp];
623 header.writeToFile(a_handle);
625 if (s_verbosity >= 3)
627 pout() << header << endl;
632 void AMRLevelPluto::writeCheckpointLevel(HDF5Handle& a_handle)
const
634 CH_assert(allDefined());
636 if (s_verbosity >= 3)
638 pout() <<
"AMRLevelPluto::writeCheckpointLevel" << endl;
643 sprintf(levelStr,
"%d",m_level);
644 const std::string label = std::string(
"level_") + levelStr;
646 a_handle.setGroup(label);
649 HDF5HeaderData header;
651 header.m_int [
"ref_ratio"] = m_ref_ratio;
652 header.m_int [
"tag_buffer_size"] = m_tagBufferSize;
653 header.m_real[
"dx"] = m_dx;
654 header.m_int [
"geometry"] =
GEOMETRY;
655 header.m_int [
"logr"] = CHOMBO_LOGR;
657 header.m_real[
"g_x2stretch"] = g_x2stretch;
660 header.m_real[
"g_x3stretch"] = g_x3stretch;
669 header.m_real[
"dt"] = m_dt;
671 header.m_real[
"ch"] = glm_ch_max;
673 header.m_real[
"time"] = m_time;
674 header.m_box [
"prob_domain"] = m_problem_domain.domainBox();
677 D_TERM(
if (m_problem_domain.isPeriodic(0))
679 header.m_int [
"is_periodic_0"] = 1;
683 header.m_int [
"is_periodic_0"] = 0;
686 if (m_problem_domain.isPeriodic(1))
688 header.m_int [
"is_periodic_1"] = 1;
692 header.m_int [
"is_periodic_1"] = 0;
695 if (m_problem_domain.isPeriodic(2))
697 header.m_int [
"is_periodic_2"] = 1;
701 header.m_int [
"is_periodic_2"] = 0;
705 header.writeToFile(a_handle);
707 if (s_verbosity >= 3)
709 pout() << header << endl;
713 #if GEOMETRY != CARTESIAN
714 LevelData<FArrayBox> tmp_U;
715 IntVect ivGhost = m_numGhost*IntVect::Unit;
716 tmp_U.define(m_grids,m_numStates,ivGhost);
718 const LevelData<FArrayBox>& dV = m_levelPluto.getdV();
720 DataIterator dit = m_UNew.dataIterator();
721 for(;dit.ok(); ++dit){
722 tmp_U[dit()].copy(m_UNew[dit()]);
723 FArrayBox& curU = tmp_U[dit()];
724 const FArrayBox& curdV = dV[dit()];
725 for (
int ivar = 0; ivar < curU.nComp(); ivar++) curU.divide(curdV,0,ivar);
726 #if CHOMBO_CONS_AM == YES
727 #if ROTATING_FRAME == YES
728 Box curBox = curU.box();
729 for(BoxIterator bit(curBox); bit.ok(); ++bit) {
730 const IntVect& iv = bit();
731 curU(iv,
iMPHI) /= curdV(iv,1);
735 curU.divide(curdV,1,
iMPHI);
739 write(a_handle,tmp_U.boxLayout());
740 write(a_handle,tmp_U,
"data");
742 if (g_stretch_fact != 1.) {
743 LevelData<FArrayBox> tmp_U;
744 IntVect ivGhost = m_numGhost*IntVect::Unit;
745 tmp_U.define(m_grids,m_numStates,ivGhost);
747 DataIterator dit = m_UNew.dataIterator();
748 for(;dit.ok(); ++dit){
749 tmp_U[dit()].copy(m_UNew[dit()]);
750 tmp_U[dit()] /= g_stretch_fact;
752 write(a_handle,tmp_U.boxLayout());
753 write(a_handle,tmp_U,
"data");
755 write(a_handle,m_UNew.boxLayout());
756 write(a_handle,m_UNew,
"data");
763 void AMRLevelPluto::readCheckpointHeader(HDF5Handle& a_handle)
765 if (s_verbosity >= 3)
767 pout() <<
"AMRLevelPluto::readCheckpointHeader" << endl;
771 HDF5HeaderData header;
772 header.readFromFile(a_handle);
774 if (s_verbosity >= 3)
776 pout() <<
"hdf5 header data:" << endl;
777 pout() << header << endl;
781 if (header.m_int.find(
"num_components") == header.m_int.end())
783 MayDay::Error(
"AMRLevelPluto::readCheckpointHeader: checkpoint file does not have num_components");
786 int numStates = header.m_int[
"num_components"];
787 if (numStates != m_numStates)
789 MayDay::Error(
"AMRLevelPluto::readCheckpointHeader: num_components in checkpoint file does not match solver");
793 std::string stateName;
795 for (
int comp = 0; comp < m_numStates; ++comp)
797 sprintf(compStr,
"component_%d",comp);
798 if (header.m_string.find(compStr) == header.m_string.end())
800 MayDay::Error(
"AMRLevelPluto::readCheckpointHeader: checkpoint file does not have enough component names");
803 stateName = header.m_string[compStr];
804 if (stateName != m_ConsStateNames[comp])
806 MayDay::Error(
"AMRLevelPluto::readCheckpointHeader: state_name in checkpoint does not match solver");
812 void AMRLevelPluto::readCheckpointLevel(HDF5Handle& a_handle)
814 if (s_verbosity >= 3)
816 pout() <<
"AMRLevelPluto::readCheckpointLevel" << endl;
821 sprintf(levelStr,
"%d",m_level);
822 const std::string label = std::string(
"level_") + levelStr;
825 a_handle.setGroup(label);
827 HDF5HeaderData header;
828 header.readFromFile(a_handle);
830 if (s_verbosity >= 3)
832 pout() <<
"hdf5 header data:" << endl;
833 pout() << header << endl;
837 if (header.m_int.find(
"ref_ratio") == header.m_int.end())
839 MayDay::Error(
"AMRLevelPluto::readCheckpointLevel: file does not contain ref_ratio");
841 m_ref_ratio = header.m_int[
"ref_ratio"];
843 if (s_verbosity >= 2)
845 pout() <<
"read ref_ratio = " << m_ref_ratio << endl;
849 if (header.m_int.find(
"tag_buffer_size") == header.m_int.end())
851 MayDay::Error(
"AMRLevelPluto::readCheckpointLevel: file does not contain tag_buffer_size");
853 m_tagBufferSize= header.m_int[
"tag_buffer_size"];
855 if (s_verbosity >= 2)
857 pout() <<
"read tag_buffer_size = " << m_tagBufferSize << endl;
861 if (header.m_real.find(
"dx") == header.m_real.end())
863 MayDay::Error(
"AMRLevelPluto::readCheckpointLevel: file does not contain dx");
865 m_dx = header.m_real[
"dx"];
867 if (s_verbosity >= 2)
869 pout() <<
"read dx = " << m_dx << endl;
873 if (header.m_real.find(
"dt") == header.m_real.end())
875 MayDay::Error(
"AMRLevelPluto::readCheckpointLevel: file does not contain dt");
877 m_dt = header.m_real[
"dt"];
879 if (s_verbosity >= 2)
881 pout() <<
"read dt = " << m_dt << endl;
886 if (header.m_real.find(
"g_x2stretch") == header.m_real.end())
888 MayDay::Error(
"AMRLevelPluto::readCheckpointLevel: file does not contain g_x2stretch");
890 g_x2stretch = header.m_real[
"g_x2stretch"];
892 if (s_verbosity >= 2)
894 pout() <<
"read g_x2stretch = " << g_x2stretch << endl;
898 if (header.m_real.find(
"g_x3stretch") == header.m_real.end())
900 MayDay::Error(
"AMRLevelPluto::readCheckpointLevel: file does not contain g_x3stretch");
902 g_x3stretch = header.m_real[
"g_x3stretch"];
904 if (s_verbosity >= 2)
906 pout() <<
"read g_x3stretch = " << g_x3stretch << endl;
910 #if GEOMETRY == CARTESIAN
911 g_stretch_fact = g_x2stretch*g_x3stretch;
916 if (header.m_real.find(
"ch") == header.m_real.end())
918 MayDay::Error(
"AMRLevelPluto::readCheckpointLevel: file does not contain glm_ch");
920 glm_ch_max = header.m_real[
"ch"];
922 if (s_verbosity >= 2)
924 pout() <<
"read ch max = " << glm_ch_max << endl;
929 if (header.m_box.find(
"prob_domain") == header.m_box.end())
931 MayDay::Error(
"AMRLevelPluto::readCheckpointLevel: file does not contain prob_domain");
934 Box domainBox = header.m_box[
"prob_domain"];
938 bool isPeriodic[SpaceDim];
939 D_TERM(
if (!(header.m_int.find(
"is_periodic_0") == header.m_int.end()))
940 isPeriodic[0] = (header.m_int["is_periodic_0"] == 1);
942 isPeriodic[0] = false; ,
944 if (!(header.m_int.
find("is_periodic_1") == header.m_int.end()))
945 isPeriodic[1] = (header.m_int["is_periodic_1"] == 1);
947 isPeriodic[1] = false; ,
949 if (!(header.m_int.
find("is_periodic_2") == header.m_int.end()))
950 isPeriodic[2] = (header.m_int["is_periodic_2"] == 1);
952 isPeriodic[2] = false;);
954 m_problem_domain = ProblemDomain(domainBox,isPeriodic);
958 const
int gridStatus =
read(a_handle,grids);
962 MayDay::Error(
"AMRLevelPluto::readCheckpointLevel: file does not contain a Vector<Box>");
966 m_grids = loadBalance(grids);
970 const DisjointBoxLayout& constGrids = m_grids;
972 LayoutIterator lit = constGrids.layoutIterator();
973 for (lit.begin(); lit.ok(); ++lit)
975 const Box& b = constGrids[lit()];
976 m_level_grids.push_back(b);
979 if (s_verbosity >= 4)
981 pout() <<
"read level domain: " << endl;
982 LayoutIterator lit = m_grids.layoutIterator();
983 for (lit.begin(); lit.ok(); ++lit)
985 const Box& b = m_grids[lit()];
986 pout() << lit().intCode() <<
": " << b << endl;
992 m_UNew.define(m_grids,m_numStates);
993 const int dataStatus = read<FArrayBox>(a_handle,
1000 MayDay::Error(
"AMRLevelPluto::readCheckpointLevel: file does not contain state data");
1003 m_UOld.define(m_grids,m_numStates);
1008 #if GEOMETRY != CARTESIAN
1009 const LevelData<FArrayBox>& dV = m_levelPluto.getdV();
1010 DataIterator dit = m_UNew.dataIterator();
1011 for(;dit.ok(); ++dit){
1012 FArrayBox& curU = m_UNew[dit()];
1013 const FArrayBox& curdV = dV[dit()];
1014 #if CHOMBO_CONS_AM == YES
1015 #if ROTATING_FRAME == YES
1016 Box curBox = curU.box();
1017 for(BoxIterator bit(curBox); bit.ok(); ++bit) {
1018 const IntVect& iv = bit();
1020 curU(iv,
iMPHI) *= curdV(iv,1);
1023 curU.mult(curdV,1,
iMPHI);
1026 for (
int ivar = 0; ivar < curU.nComp(); ivar++) curU.mult(curdV,0,ivar);
1029 if (g_stretch_fact != 1.) {
1030 DataIterator dit = m_UNew.dataIterator();
1031 for(;dit.ok(); ++dit) {
1032 FArrayBox& curU = m_UNew[dit()];
1033 curU *= g_stretch_fact;
1040 void AMRLevelPluto::writePlotHeader(HDF5Handle& a_handle)
const
1042 CH_assert(allDefined());
1044 if (s_verbosity >= 3)
1046 pout() <<
"AMRLevelPluto::writePlotHeader" << endl;
1050 HDF5HeaderData header;
1051 header.m_int[
"num_components"] = m_numStates;
1055 for (
int comp = 0; comp < m_numStates; ++comp)
1057 sprintf(compStr,
"component_%d",comp);
1058 header.m_string[compStr] = m_PrimStateNames[comp];
1062 header.writeToFile(a_handle);
1063 a_handle.setGroup(
"/Expressions");
1064 HDF5HeaderData expressions;
1065 DefineExpressions(expressions);
1066 expressions.writeToFile(a_handle);
1068 if (s_verbosity >= 3)
1070 pout() << header << endl;
1075 void AMRLevelPluto::writePlotLevel(HDF5Handle& a_handle)
const
1077 CH_assert(allDefined());
1081 sprintf(levelStr,
"%d",m_level);
1082 const std::string label = std::string(
"level_") + levelStr;
1084 a_handle.setGroup(label);
1087 HDF5HeaderData header;
1089 header.m_int [
"ref_ratio"] = m_ref_ratio;
1090 header.m_real[
"dx"] = m_dx;
1091 header.m_int [
"geometry"] =
GEOMETRY;
1092 header.m_int [
"logr"] = CHOMBO_LOGR;
1094 header.m_real[
"g_x2stretch"] = g_x2stretch;
1097 header.m_real[
"g_x3stretch"] = g_x3stretch;
1106 header.m_real[
"dt"] = m_dt;
1107 header.m_real[
"time"] = m_time;
1108 header.m_box [
"prob_domain"] = m_problem_domain.domainBox();
1111 header.writeToFile(a_handle);
1113 if (s_verbosity >= 3)
1115 pout() << header << endl;
1119 LevelData<FArrayBox> tmp_U;
1120 IntVect ivGhost = m_numGhost*IntVect::Unit;
1121 tmp_U.define(m_grids,m_numStates,ivGhost);
1123 #if GEOMETRY != CARTESIAN
1124 const LevelData<FArrayBox>& dV = m_levelPluto.getdV();
1127 DataIterator dit = m_UNew.dataIterator();
1128 for(;dit.ok(); ++dit){
1129 tmp_U[dit()].copy(m_UNew[dit()]);
1130 FArrayBox& curU = tmp_U[dit()];
1131 #if GEOMETRY != CARTESIAN
1132 const FArrayBox& curdV = dV[dit()];
1133 for (
int ivar = 0; ivar < curU.nComp(); ivar++) curU.divide(curdV,0,ivar);
1134 #if CHOMBO_CONS_AM == YES
1135 #if ROTATING_FRAME == YES
1136 Box curBox = curU.box();
1137 for(BoxIterator bit(curBox); bit.ok(); ++bit) {
1138 const IntVect& iv = bit();
1139 curU(iv,
iMPHI) /= curdV(iv,1);
1143 curU.divide(curdV,1,
iMPHI);
1147 if (g_stretch_fact != 1.) curU /= g_stretch_fact;
1150 m_patchPluto->convertFArrayBox(curU);
1152 write(a_handle,tmp_U.boxLayout());
1153 write(a_handle,tmp_U,
"data");
1157 void AMRLevelPluto::DefineExpressions(HDF5HeaderData& a_expressions)
const
1164 #if GEOMETRY == CARTESIAN
1165 sprintf (str_expr,
"coords(Mesh)[0]+nodal_constant(Mesh,%12.6e)",
g_domBeg[
IDIR]);
1166 a_expressions.m_string[
"scalar X"] = str_expr;
1168 sprintf (str_expr,
"nodal_constant(Mesh,%12.6e)*coords(Mesh)[1]+nodal_constant(Mesh,%12.6e)",g_x2stretch,
g_domBeg[
JDIR]);
1169 a_expressions.m_string[
"scalar Y"] = str_expr;
1170 a_expressions.m_string[
"vector Displacement"] =
"{X,Y}-coords(Mesh)";
1173 sprintf (str_expr,
"nodal_constant(Mesh,%12.6e)*coords(Mesh)[1]+nodal_constant(Mesh,%12.6e)",g_x2stretch,
g_domBeg[
JDIR]);
1174 a_expressions.m_string[
"scalar Y"] = str_expr;
1175 sprintf (str_expr,
"nodal_constant(Mesh,%12.6e)*coords(Mesh)[2]+nodal_constant(Mesh,%12.6e)",g_x3stretch,
g_domBeg[
KDIR]);
1176 a_expressions.m_string[
"scalar Z"] = str_expr;
1177 a_expressions.m_string[
"vector Displacement"] =
"{X,Y,Z}-coords(Mesh)";
1181 #if GEOMETRY == CYLINDRICAL
1182 sprintf (str_expr,
"coords(Mesh)[0]+nodal_constant(Mesh,%12.6e)",
g_domBeg[
IDIR]);
1183 a_expressions.m_string[
"scalar R"] = str_expr;
1185 sprintf (str_expr,
"nodal_constant(Mesh,%12.6e)*coords(Mesh)[1]+nodal_constant(Mesh,%12.6e)",g_x2stretch,
g_domBeg[
JDIR]);
1186 a_expressions.m_string[
"scalar Z"] = str_expr;
1187 a_expressions.m_string[
"vector Displacement"] =
"{R,Z}-coords(Mesh)";
1191 #if GEOMETRY == SPHERICAL
1192 #if CHOMBO_LOGR == NO
1193 sprintf (str_expr,
"coords(Mesh)[0]+nodal_constant(Mesh,%12.6e)",
g_domBeg[
IDIR]);
1194 a_expressions.m_string[
"scalar R"] = str_expr;
1196 sprintf (str_expr,
"exp(coords(Mesh)[0])*nodal_constant(Mesh,%12.6e)",
g_domBeg[
IDIR]);
1197 a_expressions.m_string[
"scalar R"] = str_expr;
1200 sprintf (str_expr,
"nodal_constant(Mesh,%12.6e)*coords(Mesh)[1]+nodal_constant(Mesh,%12.6e)",g_x2stretch,
g_domBeg[
JDIR]);
1201 a_expressions.m_string[
"scalar Theta"] = str_expr;
1202 a_expressions.m_string[
"scalar X"] =
"R*sin(Theta)";
1203 a_expressions.m_string[
"scalar Z"] =
"R*cos(Theta)";
1204 a_expressions.m_string[
"vector Displacement"] =
"{X,Z}-coords(Mesh)";
1207 sprintf (str_expr,
"nodal_constant(Mesh,%12.6e)*coords(Mesh)[1]+nodal_constant(Mesh,%12.6e)",g_x2stretch,
g_domBeg[
JDIR]);
1208 a_expressions.m_string[
"scalar Theta"] = str_expr;
1209 sprintf (str_expr,
"nodal_constant(Mesh,%12.6e)*coords(Mesh)[2]+nodal_constant(Mesh,%12.6e)",g_x3stretch,
g_domBeg[
KDIR]);
1210 a_expressions.m_string[
"scalar Phi"] = str_expr;
1211 a_expressions.m_string[
"scalar X"] =
"R*sin(Theta)*cos(Phi)";
1212 a_expressions.m_string[
"scalar Y"] =
"R*sin(Theta)*sin(Phi)";
1213 a_expressions.m_string[
"scalar Z"] =
"R*cos(Theta)";
1214 a_expressions.m_string[
"vector Displacement"] =
"{X,Y,Z}-coords(Mesh)";
1218 #if GEOMETRY == POLAR
1219 #if CHOMBO_LOGR == NO
1220 sprintf (str_expr,
"coords(Mesh)[0]+nodal_constant(Mesh,%12.6e)",
g_domBeg[
IDIR]);
1221 a_expressions.m_string[
"scalar R"] = str_expr;
1223 sprintf (str_expr,
"exp(coords(Mesh)[0])*nodal_constant(Mesh,%12.6e)",
g_domBeg[
IDIR]);
1224 a_expressions.m_string[
"scalar R"] = str_expr;
1227 sprintf (str_expr,
"nodal_constant(Mesh,%12.6e)*coords(Mesh)[1]+nodal_constant(Mesh,%12.6e)",g_x2stretch,
g_domBeg[
JDIR]);
1228 a_expressions.m_string[
"scalar Phi"] = str_expr;
1229 a_expressions.m_string[
"scalar X"] =
"R*cos(Phi)";
1230 a_expressions.m_string[
"scalar Y"] =
"R*sin(Phi)";
1231 a_expressions.m_string[
"vector Displacement"] =
"{X,Y}-coords(Mesh)";
1234 sprintf (str_expr,
"nodal_constant(Mesh,%12.6e)*coords(Mesh)[1]+nodal_constant(Mesh,%12.6e)",g_x2stretch,
g_domBeg[
JDIR]);
1235 a_expressions.m_string[
"scalar Phi"] = str_expr;
1236 a_expressions.m_string[
"scalar X"] =
"R*cos(Phi)";
1237 a_expressions.m_string[
"scalar Y"] =
"R*sin(Phi)";
1238 sprintf (str_expr,
"nodal_constant(Mesh,%12.6e)*coords(Mesh)[2]+nodal_constant(Mesh,%12.6e)",g_x3stretch,
g_domBeg[
KDIR]);
1239 a_expressions.m_string[
"scalar Z"] = str_expr;
1240 a_expressions.m_string[
"vector Displacement"] =
"{X,Y,Z}-coords(Mesh)";
1249 Real AMRLevelPluto::computeDt()
1251 CH_assert(allDefined());
1253 if (s_verbosity >= 3)
1255 pout() <<
"AMRLevelPluto::computeDt " << m_level << endl;
1265 Real AMRLevelPluto::computeInitialDt()
1267 CH_assert(allDefined());
1269 if (s_verbosity >= 3)
1271 pout() <<
"AMRLevelPluto::computeInitialDt " << m_level << endl;
1274 Real newDT = m_initial_dt_multiplier*m_dx/m_domainLength;
1279 const LevelData<FArrayBox>& AMRLevelPluto::getStateNew()
const
1281 CH_assert(allDefined());
1286 const LevelData<FArrayBox>& AMRLevelPluto::getStateOld()
const
1288 CH_assert(allDefined());
1293 bool AMRLevelPluto::allDefined()
const
1295 return isDefined() &&
1300 DisjointBoxLayout AMRLevelPluto::loadBalance(
const Vector<Box>& a_grids)
1302 CH_assert(allDefined());
1305 Vector<int> procMap;
1311 Vector<long long> computeLoads(a_grids.size());
1313 for (
int i = 0;
i < a_grids.size(); ++
i) {
1315 computeLoads[
i] = a_grids[
i].numPts();
1318 LoadBalance(procMap, computeLoads, a_grids, numProc());
1323 if (s_verbosity >= 4)
1325 pout() <<
"AMRLevelPluto::loadBalance (level "<<m_level<<
"): procesor map: " << endl;
1327 for (
int igrid = 0; igrid < a_grids.size(); ++igrid)
1329 pout() << igrid <<
"(load = " << computeLoads[igrid] <<
1330 "): " << procMap[igrid] <<
" " << endl;
1335 DisjointBoxLayout dbl(a_grids,procMap,m_problem_domain);
1342 void AMRLevelPluto::levelSetup()
1344 CH_assert(allDefined());
1346 if (s_verbosity >= 3)
1348 pout() <<
"AMRLevelPluto::levelSetup " << m_level << endl;
1351 AMRLevelPluto* amrGodCoarserPtr = getCoarserLevel();
1352 AMRLevelPluto* amrGodFinerPtr = getFinerLevel();
1354 m_hasCoarser = (amrGodCoarserPtr != NULL);
1355 m_hasFiner = (amrGodFinerPtr != NULL);
1358 #ifdef SKIP_SPLIT_CELLS
1360 m_split_tags.define(m_grids,1);
1361 for (DataIterator dit = m_grids.dataIterator(); dit.ok(); ++dit) {
1362 m_split_tags[dit()].setVal(1.);
1368 int nRefCrse = m_coarser_level_ptr->refRatio();
1370 m_coarseAverage.define(m_grids,
1374 m_fineInterp.define(m_grids,
1380 const DisjointBoxLayout& coarserLevelDomain = amrGodCoarserPtr->m_grids;
1382 #ifdef SKIP_SPLIT_CELLS
1384 amrGodCoarserPtr->mark_split(m_grids);
1388 m_levelPluto.define(m_grids,
1394 m_patchPlutoFactory,
1405 amrGodCoarserPtr->m_fluxRegister.define(m_grids,
1406 amrGodCoarserPtr->m_grids,
1408 amrGodCoarserPtr->m_ref_ratio,
1410 amrGodCoarserPtr->m_fluxRegister.setToZero();
1414 m_levelPluto.define(m_grids,
1415 DisjointBoxLayout(),
1420 m_patchPlutoFactory,
1427 AMRLevelPluto* AMRLevelPluto::getCoarserLevel()
const
1429 CH_assert(allDefined());
1431 AMRLevelPluto* amrGodCoarserPtr = NULL;
1433 if (m_coarser_level_ptr != NULL)
1435 amrGodCoarserPtr =
dynamic_cast<AMRLevelPluto*
>(m_coarser_level_ptr);
1437 if (amrGodCoarserPtr == NULL)
1439 MayDay::Error(
"AMRLevelPluto::getCoarserLevel: dynamic cast failed");
1443 return amrGodCoarserPtr;
1447 AMRLevelPluto* AMRLevelPluto::getFinerLevel()
const
1449 CH_assert(allDefined());
1451 AMRLevelPluto* amrGodFinerPtr = NULL;
1453 if (m_finer_level_ptr != NULL)
1455 amrGodFinerPtr =
dynamic_cast<AMRLevelPluto*
>(m_finer_level_ptr);
1457 if (amrGodFinerPtr == NULL)
1459 MayDay::Error(
"AMRLevelPluto::getFinerLevel: dynamic cast failed");
1463 return amrGodFinerPtr;
1467 void AMRLevelPluto::mark_split(
const DisjointBoxLayout& finerLevelDomain)
1469 CH_assert(allDefined());
1471 DisjointBoxLayout finerDomain;
1473 coarsen(finerDomain,finerLevelDomain,m_ref_ratio);
1475 DataIterator dit = m_grids.dataIterator();
1476 LayoutIterator lit = finerDomain.layoutIterator();
1477 for (dit.begin(); dit.ok(); ++dit)
1479 FArrayBox& splitU = m_split_tags[dit()];
1481 Box splitBox = splitU.box();
1483 for (lit.begin(); lit.ok(); ++lit)
1485 Box finBox = finerDomain.get(lit());
1486 const Box intBox = splitBox & finBox;
1488 BoxIterator bit(intBox);
1489 for (bit.begin(); bit.ok(); ++bit)
1491 const IntVect& iv = bit();
1499 #include "NamespaceFooter.H"
int g_intStage
Gives the current integration stage of the time stepping method (predictor = 0, 1st corrector = 1...
def find(fname, str, action, want)
double g_domBeg[3]
Lower limits of the computational domain.