REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_Tagging.cpp
Go to the documentation of this file.
1#include <REMORA.H>
2#include <REMORA_Derive.H>
3
4using namespace amrex;
5
6/**
7 * Function to tag cells for refinement -- this overrides the pure virtual function in AmrCore
8 *
9 * @param[in] levc level of refinement (0 is coarsest level)
10 * @param[out] tags array of tagged cells
11 * @param[in] time current time
12 * @param[in] ngrow number of grow cells
13*/
14void
15REMORA::ErrorEst (int levc, TagBoxArray& tags, Real time, int /*ngrow*/)
16{
17 const int clearval = TagBox::CLEAR;
18 const int tagval = TagBox::SET;
19
20 //
21 // This mf must have ghost cells because we may take differences between adjacent values
22 //
23 std::unique_ptr<MultiFab> mf = std::make_unique<MultiFab>(grids[levc], dmap[levc], 1, 1);
24
25 for (int j=0; j < ref_tags.size(); ++j)
26 {
27 if (ref_tags[j].Field() == "tracer" || ref_tags[j].Field() == "temp" ||
28 ref_tags[j].Field() == "salt") {
29 FillPatch(levc, time, *cons_new[levc], cons_new, BCVars::cons_bc, BdyVars::t,
30 0,true,false);
31 }
32 // This allows dynamic refinement based on the value of the tracer
33 if (ref_tags[j].Field() == "tracer")
34 {
35 MultiFab::Copy(*mf,*cons_new[levc],Tracer_comp,0,1,1);
36 } else if (ref_tags[j].Field() == "temp") {
37 MultiFab::Copy(*mf,*cons_new[levc],Temp_comp,0,1,1);
38 } else if (ref_tags[j].Field() == "salt") {
39 MultiFab::Copy(*mf,*cons_new[levc],Salt_comp,0,1,1);
40 } else if (ref_tags[j].Field() == "x_velocity") {
41 FillPatch(levc, time, *xvel_new[levc], xvel_new, BCVars::xvel_bc, BdyVars::u,0,true,true);
42 MultiFab::Copy(*mf,*xvel_new[levc],0,0,1,1);
43 } else if (ref_tags[j].Field() == "y_velocity") {
44 FillPatch(levc, time, *yvel_new[levc], yvel_new, BCVars::yvel_bc, BdyVars::v,0,true,true);
45 MultiFab::Copy(*mf,*yvel_new[levc],0,0,1,1);
46 } else if (ref_tags[j].Field() == "z_velocity") {
47 FillPatch(levc, time, *zvel_new[levc], zvel_new, BCVars::zvel_bc, BdyVars::null,0,true,true);
48 MultiFab::Copy(*mf,*zvel_new[levc],0,0,1,1);
49 } else if (ref_tags[j].Field() == "vorticity") {
50 MultiFab mf_cc_vel(grids[levc],dmap[levc],3,1);
51 average_face_to_cellcenter(mf_cc_vel,0,
52 Array<const MultiFab*,3>{xvel_new[levc],
53 yvel_new[levc],
54 zvel_new[levc]});
55 // Impose bc's at domain boundaries at all levels
56 FillBdyCCVels(levc, mf_cc_vel);
57
58#ifdef _OPENMP
59#pragma omp parallel if (Gpu::notInLaunchRegion())
60#endif
61 for (MFIter mfi(*mf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
62 {
63 const Box& bx = mfi.tilebox();
64 auto& dfab = (*mf)[mfi];
65 auto& sfab = mf_cc_vel[mfi];
66 auto pm = vec_pm[levc]->const_array(mfi);
67 auto pn = vec_pn[levc]->const_array(mfi);
68 auto maskr = vec_mskr[levc]->const_array(mfi);
69 derived::remora_dervort(bx, dfab, 0, 1, sfab, pm, pn, maskr, Geom(levc), time, nullptr, levc);
70 } // mfi
71
72 mf->FillBoundary(geom[levc].periodicity());
73 //
74 // TODO: we may need to fill physical boundaries here before tagging criteria are imposed
75 //
76
77 } else if (ref_tags[j].Field() == "mask") {
78 MultiFab::Copy(*mf,*vec_mskr3d[levc],0,0,1,IntVect(1,1,0));
79#ifdef REMORA_USE_PARTICLES
80 } else {
81 //
82 // This allows dynamic refinement based on the number of particles per cell
83 //
84 // Note that we must count all the particles in levels both at and above the current,
85 // since otherwise, e.g., if the particles are all at level 1, counting particles at
86 // level 0 will not trigger refinement when regridding so level 1 will disappear,
87 // then come back at the next regridding
88 //
89 const auto& particles_namelist( particleData.getNames() );
90 mf->setVal(0.0);
91 for (ParticlesNamesVector::size_type i = 0; i < particles_namelist.size(); i++)
92 {
93 std::string tmp_string(particles_namelist[i]+"_count");
94 IntVect rr = IntVect::TheUnitVector();
95 if (ref_tags[j].Field() == tmp_string) {
96 for (int lev = levc; lev <= finest_level; lev++)
97 {
98 MultiFab temp_dat(grids[lev], dmap[lev], 1, 0); temp_dat.setVal(0);
99 particleData[particles_namelist[i]]->IncrementWithTotal(temp_dat, lev);
100
101 MultiFab temp_dat_crse(grids[levc], dmap[levc], 1, 0); temp_dat_crse.setVal(0);
102
103 if (lev == levc) {
104 MultiFab::Copy(*mf, temp_dat, 0, 0, 1, 0);
105 } else {
106 for (int d = 0; d < AMREX_SPACEDIM; d++) {
107 rr[d] *= ref_ratio[levc][d];
108 }
109 average_down(temp_dat, temp_dat_crse, 0, 1, rr);
110 MultiFab::Add(*mf, temp_dat_crse, 0, 0, 1, 0);
111 }
112 }
113 }
114 }
115
116
117#endif
118 }
119
120 ref_tags[j](tags,mf.get(),clearval,tagval,time,levc,geom[levc]);
121 }
122}
123
124/**
125 * Function to define the refinement criteria based on user input
126*/
127void
128REMORA::refinement_criteria_setup ()
129{
130 if (max_level > 0)
131 {
132 ParmParse pp(pp_prefix);
133 Vector<std::string> refinement_indicators;
134 pp.queryarr("refinement_indicators",refinement_indicators,0,pp.countval("refinement_indicators"));
135
136 for (int i=0; i<refinement_indicators.size(); ++i)
137 {
138 std::string ref_prefix = pp_prefix + "." + refinement_indicators[i];
139
140 ParmParse ppr(ref_prefix);
141 RealBox realbox;
142 int lev_for_box;
143
144 int num_real_lo = ppr.countval("in_box_lo");
145 int num_indx_lo = ppr.countval("in_box_lo_indices");
146 int num_indx_lo_crse = ppr.countval("in_box_lo_indices_crse");
147
148 int num_real_hi = ppr.countval("in_box_hi");
149 int num_indx_hi = ppr.countval("in_box_hi_indices");
150 int num_indx_hi_crse = ppr.countval("in_box_hi_indices_crse");
151
152 AMREX_ALWAYS_ASSERT( (num_real_lo == num_real_hi) && (num_real_lo == 0 || num_real_lo >= 2) );
153 AMREX_ALWAYS_ASSERT( (num_indx_lo == num_indx_hi) && (num_indx_lo == 0 || num_indx_lo >= 2) );
154 AMREX_ALWAYS_ASSERT( (num_indx_lo_crse == num_indx_hi_crse) && (num_indx_lo_crse == 0 || num_indx_lo_crse >= 2) );
155
156 // Problem low and high (in real not index space) are the same at all levels
157 if ( !((num_real_lo >= AMREX_SPACEDIM-1 && num_indx_lo == 0 && num_indx_lo_crse == 0) ||
158 (num_indx_lo >= AMREX_SPACEDIM-1 && num_real_lo == 0 && num_indx_lo_crse == 0) ||
159 (num_indx_lo == 0 && num_real_lo == 0 && num_indx_lo_crse == 0) ||
160 (num_indx_lo_crse >= AMREX_SPACEDIM-1 && num_real_lo == 0 && num_indx_lo == 0)
161 ) )
162 {
163 amrex::Abort("Must only specify box for refinement using real OR index space with fine/coarse grid indices");
164 }
165
166 if (num_real_lo > 0) {
167 std::vector<Real> box_lo(3), box_hi(3);
168 ppr.get("max_level",lev_for_box);
169 if (lev_for_box > 0 && lev_for_box <= max_level)
170 {
171 ppr.getarr("in_box_lo",box_lo,0,2);
172 ppr.getarr("in_box_hi",box_hi,0,2);
173 box_lo[2] = geom[0].ProbLo(2);
174 box_hi[2] = geom[0].ProbHi(2);
175 realbox = RealBox(&(box_lo[0]),&(box_hi[0]));
176
177 amrex::Print() << "Reading " << realbox << " at level " << lev_for_box << std::endl;
178 num_boxes_at_level[lev_for_box] += 1;
179
180 const auto* dx = geom[lev_for_box].CellSize();
181 const Real* plo = geom[lev_for_box].ProbLo();
182
183 int ilo = static_cast<int>((box_lo[0] - plo[0])/dx[0]);
184 int jlo = static_cast<int>((box_lo[1] - plo[1])/dx[1]);
185 int klo = static_cast<int>((box_lo[2] - plo[2])/dx[2]);
186 int ihi = static_cast<int>((box_hi[0] - plo[0])/dx[0]-1);
187 int jhi = static_cast<int>((box_hi[1] - plo[1])/dx[1]-1);
188 int khi = static_cast<int>((box_hi[2] - plo[2])/dx[2]-1);
189
190 Box bx_old(IntVect(ilo,jlo,klo),IntVect(ihi,jhi,khi));
191
192 int mod_ilo = ilo%ref_ratio[lev_for_box-1][0];
193 int mod_jlo = jlo%ref_ratio[lev_for_box-1][1];
194
195 int mod_ihi = (ihi+1)%ref_ratio[lev_for_box-1][0];
196 int mod_jhi = (jhi+1)%ref_ratio[lev_for_box-1][1];
197
198 if (mod_ilo != 0) {
199 ilo -= mod_ilo;
200 }
201 if (mod_jlo != 0) {
202 jlo -= mod_jlo;
203 }
204 if (mod_ihi != 0) {
205 ihi += ref_ratio[lev_for_box-1][0] - mod_ihi;
206 }
207 if (mod_jhi != 0) {
208 jhi += ref_ratio[lev_for_box-1][1] - mod_jhi;
209 }
210 Box bx(IntVect(ilo,jlo,klo),IntVect(ihi,jhi,khi));
211 if (mod_ilo !=0 || mod_jlo !=0 || mod_ihi != 0 || mod_jhi != 0) {
212 amrex::Print() << "Fine box on level " << lev_for_box << " adjusted from " << bx_old << " to " << bx << " to make it valid for refinement." << std::endl;
213 }
214 boxes_at_level[lev_for_box].push_back(bx);
215 amrex::Print() << "Saving in 'boxes at level' as " << bx << std::endl;
216 } // lev
217
218 } else if (num_indx_lo > 0) {
219
220 std::vector<int> box_lo(3), box_hi(3);
221 ppr.get("max_level",lev_for_box);
222 if (lev_for_box > 0 && lev_for_box <= max_level)
223 {
224 if (n_error_buf[0] != IntVect::TheZeroVector()) {
225 amrex::Abort("Don't use n_error_buf > 0 when setting the box explicitly");
226 }
227
228 ppr.getarr("in_box_lo_indices",box_lo,0,num_indx_lo);
229 ppr.getarr("in_box_hi_indices",box_hi,0,num_indx_hi);
230
231 if (num_indx_lo < AMREX_SPACEDIM) {
232 box_lo[2] = geom[lev_for_box].Domain().smallEnd(2);
233 box_hi[2] = geom[lev_for_box].Domain().bigEnd(2);
234 }
235
236 Box bx(IntVect(box_lo[0],box_lo[1],box_lo[2]),IntVect(box_hi[0],box_hi[1],box_hi[2]));
237 const Box& domain = geom[lev_for_box].Domain();
238
239 if (!domain.contains(bx)) {
240 amrex::Print() << "\n";
241 amrex::Print() << "Box specified is " << bx << std::endl;
242 amrex::Print() << "But domain at level is " << domain << std::endl;
243 amrex::Error("Specified box doesn't fit in the domain");
244 }
245
246 const auto* dx = geom[lev_for_box].CellSize();
247 const Real* plo = geom[lev_for_box].ProbLo();
248 realbox = RealBox(plo[0]+ box_lo[0] *dx[0], plo[1]+ box_lo[1] *dx[1], plo[2]+ box_lo[2] *dx[2],
249 plo[0]+(box_hi[0]+1)*dx[0], plo[1]+(box_hi[1]+1)*dx[1], plo[2]+(box_hi[2]+1)*dx[2]);
250
251 Print() << "Reading " << bx << " at level " << lev_for_box << std::endl;
252 num_boxes_at_level[lev_for_box] += 1;
253
254 if(box_lo[0]%ref_ratio[lev_for_box-1][0] != 0){
255 amrex::Print()<< "Requested ilo in x-direction : " << box_lo[0] << std::endl;
256 amrex::Print() << "ilo = " << box_lo[0] << " is not divisible by ref_ratio in x direction = " <<
257 ref_ratio[lev_for_box-1][0] << std::endl;
258 amrex::Error("Adjust in_box_lo_indices in x-direction to be divisible by ref_ratio and try again");
259 }
260 if((box_hi[0]+1)%ref_ratio[lev_for_box-1][0] != 0){
261 amrex::Print()<< "Requested ihi in x-direction : " << box_hi[0] << std::endl;
262 amrex::Print() << "ihi+1 = " << box_hi[0]+1 << " is not divisible by ref_ratio in x direction = " <<
263 ref_ratio[lev_for_box-1][0] << std::endl;
264 amrex::Error("Adjust in_box_hi_indices in x-direction to be divisible by ref_ratio and try again");
265 }
266 if(box_lo[1]%ref_ratio[lev_for_box-1][1] != 0){
267 amrex::Print()<< "Requested jlo in y-direction : " << box_lo[1] << std::endl;
268 amrex::Print() << "jlo = " << box_lo[1] << " is not divisible by ref_ratio in y direction = " <<
269 ref_ratio[lev_for_box-1][1] << std::endl;
270 amrex::Error("Adjust in_box_lo_indices in y-direction to be divisible by ref_ratio and try again");
271 }
272 if((box_hi[1]+1)%ref_ratio[lev_for_box-1][1] != 0){
273 amrex::Print()<< "Requested jhi in y-direction : " << box_hi[1] << std::endl;
274 amrex::Print() << "jhi+1 = " << box_hi[1]+1 << " is not divisible by ref_ratio in y direction = " <<
275 ref_ratio[lev_for_box-1][1] << std::endl;
276 amrex::Error("Adjust in_box_hi_indices in y-direction to be divisible by ref_ratio and try again");
277 }
278 if(box_lo[2]%ref_ratio[lev_for_box-1][2] != 0){
279 amrex::Print()<< "Requested klo in z-direction : " << box_lo[2] << std::endl;
280 amrex::Print() << "klo = " << box_lo[2] << " is not divisible by ref_ratio in z direction = " <<
281 ref_ratio[lev_for_box-1][2] << std::endl;
282 amrex::Error("Adjust in_box_lo_indices in z-direction to be divisible by ref_ratio and try again");
283 }
284 if((box_hi[2]+1)%ref_ratio[lev_for_box-1][2] != 0){
285 amrex::Print()<< "Requested khi in z-direction : " << box_hi[2] << std::endl;
286 amrex::Print() << "khi+1 = " << box_hi[2]+1 << " is not divisible by ref_ratio in z direction = " <<
287 ref_ratio[lev_for_box-1][2] << std::endl;
288 amrex::Error("Adjust in_box_hi_indices in z-direction to be divisible by ref_ratio and try again");
289 }
290
291 boxes_at_level[lev_for_box].push_back(bx);
292 Print() << "Saving in 'boxes at level' as " << bx << std::endl;
293 } // lev
294
295 } else if (num_indx_lo_crse > 0) {
296
297 std::vector<int> box_lo(3), box_hi(3);
298 ppr.get("max_level",lev_for_box);
299 if (lev_for_box > 0 && lev_for_box <= max_level)
300 {
301 if (n_error_buf[0] != IntVect::TheZeroVector()) {
302 amrex::Abort("Don't use n_error_buf > 0 when setting the box explicitly");
303 }
304
305 ppr.getarr("in_box_lo_indices_crse",box_lo,0,num_indx_lo_crse);
306 ppr.getarr("in_box_hi_indices_crse",box_hi,0,num_indx_hi_crse);
307
308 if (num_indx_lo_crse < AMREX_SPACEDIM) {
309 box_lo[2] = geom[lev_for_box-1].Domain().smallEnd(2);
310 box_hi[2] = geom[lev_for_box-1].Domain().bigEnd(2);
311 }
312
313 Box bx(IntVect(box_lo[0],box_lo[1],box_lo[2]),IntVect(box_hi[0],box_hi[1],box_hi[2]));
314
315 if (!geom[lev_for_box-1].Domain().contains(bx)) {
316 amrex::Print() << "\n";
317 amrex::Print() << "(Coarse) Box specified is " << bx << std::endl;
318 amrex::Print() << "But (coarse) domain at level is " << geom[lev_for_box-1].Domain() << std::endl;
319 amrex::Error("Specified box doesn't fit in the domain");
320 }
321
322 bx.refine(ref_ratio[lev_for_box-1]);
323
324 const auto* dx = geom[lev_for_box-1].CellSize();
325
326 const Real* plo = geom[lev_for_box].ProbLo();
327 realbox = RealBox(plo[0]+ box_lo[0] *dx[0], plo[1]+ box_lo[1] *dx[1], plo[2]+ box_lo[2] *dx[2],
328 plo[0]+(box_hi[0]+1)*dx[0], plo[1]+(box_hi[1]+1)*dx[1], plo[2]+(box_hi[2]+1)*dx[2]);
329
330 Print() << "Reading " << bx << " at level " << lev_for_box << std::endl;
331 num_boxes_at_level[lev_for_box] += 1;
332
333 boxes_at_level[lev_for_box].push_back(bx);
334 Print() << "Saving in 'boxes at level' as " << bx << std::endl;
335 } // lev
336 }
337
338 AMRErrorTagInfo info;
339
340 if (realbox.ok()) {
341 info.SetRealBox(realbox);
342 }
343 if (ppr.countval("start_time") > 0) {
344 Real ref_min_time; ppr.get("start_time",ref_min_time);
345 info.SetMinTime(ref_min_time);
346 }
347 if (ppr.countval("end_time") > 0) {
348 Real ref_max_time; ppr.get("end_time",ref_max_time);
349 info.SetMaxTime(ref_max_time);
350 }
351 if (ppr.countval("max_level") > 0) {
352 int ref_max_level; ppr.get("max_level",ref_max_level);
353 info.SetMaxLevel(ref_max_level);
354 }
355
356 if (ppr.countval("value_greater")) {
357 int num_val = ppr.countval("value_greater");
358 Vector<Real> value(num_val);
359 ppr.getarr("value_greater",value,0,num_val);
360 std::string field; ppr.get("field_name",field);
361 ref_tags.push_back(AMRErrorTag(value,AMRErrorTag::GREATER,field,info));
362 }
363 else if (ppr.countval("value_less")) {
364 int num_val = ppr.countval("value_less");
365 Vector<Real> value(num_val);
366 ppr.getarr("value_less",value,0,num_val);
367 std::string field; ppr.get("field_name",field);
368 ref_tags.push_back(AMRErrorTag(value,AMRErrorTag::LESS,field,info));
369 }
370 else if (ppr.countval("adjacent_difference_greater")) {
371 int num_val = ppr.countval("adjacent_difference_greater");
372 Vector<Real> value(num_val);
373 ppr.getarr("adjacent_difference_greater",value,0,num_val);
374 std::string field; ppr.get("field_name",field);
375 ref_tags.push_back(AMRErrorTag(value,AMRErrorTag::GRAD,field,info));
376 }
377 else if (realbox.ok())
378 {
379 ref_tags.push_back(AMRErrorTag(info));
380 } else {
381 Abort(std::string("Unrecognized refinement indicator for " + refinement_indicators[i]).c_str());
382 }
383 } // loop over criteria
384 // Untag anywhere we have masks
385 AMRErrorTagInfo info;
386 info.SetDerefine(1);
387 Real value = 0.5_rt;
388 ref_tags.push_back(AMRErrorTag(value,AMRErrorTag::LESS,"mask",info));
389 } // if max_level > 0
390}
391
REMORA_Derive.H
Temp_comp
#define Temp_comp
Definition REMORA_IndexDefines.H:8
Tracer_comp
#define Tracer_comp
Definition REMORA_IndexDefines.H:10
Salt_comp
#define Salt_comp
Definition REMORA_IndexDefines.H:9
REMORA::vec_pm
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:400
REMORA::cons_new
amrex::Vector< amrex::MultiFab * > cons_new
multilevel data container for current step's scalar data: temperature, salinity, passive tracer
Definition REMORA.H:241
REMORA::vec_mskr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskr
land/sea mask at cell centers (2D)
Definition REMORA.H:389
REMORA::zvel_new
amrex::Vector< amrex::MultiFab * > zvel_new
multilevel data container for current step's z velocities (largely unused; W stored separately)
Definition REMORA.H:247
REMORA::boxes_at_level
amrex::Vector< amrex::Vector< amrex::Box > > boxes_at_level
the boxes specified at each level by tagging criteria
Definition REMORA.H:1204
REMORA::ref_tags
static amrex::Vector< amrex::AMRErrorTag > ref_tags
Holds info for dynamically generated tagging criteria.
Definition REMORA.H:1431
REMORA::yvel_new
amrex::Vector< amrex::MultiFab * > yvel_new
multilevel data container for current step's y velocities (v in ROMS)
Definition REMORA.H:245
REMORA::num_boxes_at_level
amrex::Vector< int > num_boxes_at_level
how many boxes specified at each level by tagging criteria
Definition REMORA.H:1200
REMORA::xvel_new
amrex::Vector< amrex::MultiFab * > xvel_new
multilevel data container for current step's x velocities (u in ROMS)
Definition REMORA.H:243
REMORA::refinement_criteria_setup
void refinement_criteria_setup()
Set refinement criteria.
Definition REMORA_Tagging.cpp:128
REMORA::ErrorEst
virtual void ErrorEst(int lev, amrex::TagBoxArray &tags, amrex::Real time, int ngrow) override
Tag cells for refinement.
Definition REMORA_Tagging.cpp:15
REMORA::FillPatch
void FillPatch(int lev, amrex::Real time, amrex::MultiFab &mf_to_be_filled, amrex::Vector< amrex::MultiFab * > const &mfs, const int bccomp, const int bdy_var_type=BdyVars::null, const int icomp=0, const bool fill_all=true, const bool fill_set=true, const int n_not_fill=0, const int icomp_calc=0, const amrex::Real dt=amrex::Real(0.0), const amrex::MultiFab &mf_calc=amrex::MultiFab())
Fill a new MultiFab by copying in phi from valid region and filling ghost cells.
Definition REMORA_FillPatch.cpp:30
REMORA::FillBdyCCVels
void FillBdyCCVels(int lev, amrex::MultiFab &mf_cc_vel)
Fill the physical boundary conditions for cell-centered velocity (diagnostic only)
Definition REMORA_FillPatch.cpp:601
REMORA::pp_prefix
std::string pp_prefix
default prefix for input file parameters
Definition REMORA.H:229
REMORA::vec_pn
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:402
REMORA::vec_mskr3d
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskr3d
land/sea mask at cell centers, copied to all z levels (3D)
Definition REMORA.H:397
derived::remora_dervort
void remora_dervort(const amrex::Box &bx, amrex::FArrayBox &derfab, int dcomp, int ncomp, const amrex::FArrayBox &datfab, const amrex::Array4< const amrex::Real > &pm, const amrex::Array4< const amrex::Real > &pn, const amrex::Array4< const amrex::Real > &, const amrex::Geometry &, amrex::Real, const int *, const int)
Definition REMORA_Derive.cpp:42