REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_init_bcs.cpp
Go to the documentation of this file.
1#include <AMReX_Vector.H>
2#include <AMReX_BC_TYPES.H>
3#include <AMReX_ParmParse.H>
4
5#include <REMORA.H>
6
7using namespace amrex;
8
10{
11 auto f_set_var_bc = [this] (ParmParse& pp, int bcvar_type, Orientation ori, std::string bc_type_string) {
12 if (bc_type_string == "symmetry")
13 {
14 phys_bc_type[bcvar_type][ori] = REMORA_BC::symmetry;
15 phys_bc_need_data[bdy_index[bcvar_type]][ori] = false;
16 domain_bc_type[ori] = "Symmetry";
17 }
18 else if (bc_type_string == "outflow")
19 {
20 phys_bc_type[bcvar_type][ori] = REMORA_BC::outflow;
21 phys_bc_need_data[bdy_index[bcvar_type]][ori] = false;
22 domain_bc_type[ori] = "Outflow";
23 }
24 else if (bc_type_string == "inflow")
25 {
26 phys_bc_type[bcvar_type][ori] = REMORA_BC::inflow;
27 phys_bc_need_data[bdy_index[bcvar_type]][ori] = false;
28 domain_bc_type[ori] = "Inflow";
29
30 if (bcvar_type == BCVars::xvel_bc || bcvar_type == BCVars::yvel_bc ||
31 bcvar_type == BCVars::zvel_bc) {
32 std::vector<Real> v;
33 pp.getarr("velocity", v, 0, AMREX_SPACEDIM);
34 m_bc_extdir_vals[bcvar_type][ori] = v[bcvar_type - BCVars::xvel_bc];
35 } else if (bcvar_type == BCVars::Scalar_bc_comp) {
36 Real scalar_in = 0.;
37 if (pp.query("scalar", scalar_in))
39 }
40 }
41 else if (bc_type_string == "noslipwall")
42 {
43 phys_bc_type[bcvar_type][ori] = REMORA_BC::no_slip_wall;
44 phys_bc_need_data[bdy_index[bcvar_type]][ori] = false;
45 domain_bc_type[ori] = "NoSlipWall";
46
47 if (bcvar_type == BCVars::xvel_bc || bcvar_type == BCVars::yvel_bc ||
48 bcvar_type == BCVars::zvel_bc) {
49 std::vector<Real> v;
50
51 // The values of m_bc_extdir_vals default to 0.
52 // But if we find "velocity" in the inputs file, use those values instead.
53 if (pp.queryarr("velocity", v, 0, AMREX_SPACEDIM))
54 {
55 v[ori.coordDir()] = 0.0_rt;
56 m_bc_extdir_vals[bcvar_type][ori] = v[bcvar_type - BCVars::xvel_bc];
57 }
58 }
59 }
60 else if (bc_type_string == "slipwall")
61 {
62 phys_bc_type[bcvar_type][ori] = REMORA_BC::slip_wall;
63 phys_bc_need_data[bdy_index[bcvar_type]][ori] = false;
64 domain_bc_type[ori] = "SlipWall";
65 }
66 else if (bc_type_string == "clamped")
67 {
68 phys_bc_type[bcvar_type][ori] = REMORA_BC::clamped;
69 phys_bc_need_data[bdy_index[bcvar_type]][ori] = true;
70 domain_bc_type[ori] = "Clamped";
72 }
73 else if (bc_type_string == "chapman")
74 {
75 phys_bc_type[bcvar_type][ori] = REMORA_BC::chapman;
76 phys_bc_need_data[bdy_index[bcvar_type]][ori] = true;
77 domain_bc_type[ori] = "Chapman";
79
80 if (bcvar_type != BCVars::zeta_bc) {
81 amrex::Abort("Chapman BC can only be applied to zeta");
82 }
83 }
84 else if (bc_type_string == "flather")
85 {
86 phys_bc_type[bcvar_type][ori] = REMORA_BC::flather;
87 phys_bc_need_data[bdy_index[bcvar_type]][ori] = true;
88 domain_bc_type[ori] = "Flather";
90
91 if (!(bcvar_type == BCVars::ubar_bc || bcvar_type == BCVars::vbar_bc)) {
92 amrex::Abort("Flather BC can only be applied to ubar or vbar");
93 }
94 }
95 else if (bc_type_string == "orlanski_rad")
96 {
97 phys_bc_type[bcvar_type][ori] = REMORA_BC::orlanski_rad;
98 phys_bc_need_data[bdy_index[bcvar_type]][ori] = false;
99 domain_bc_type[ori] = "Orlanski Radiation";
100 }
101 else if (bc_type_string == "orlanski_rad_nudg")
102 {
104 phys_bc_need_data[bdy_index[bcvar_type]][ori] = true;
105 domain_bc_type[ori] = "Orlanski Radiation with nudging";
107 }
108 else if (bc_type_string == "periodic")
109 {
110 if (!geom[0].isPeriodic(ori.coordDir())) {
111 amrex::Abort("Periodic boundary specified in a non-periodic direction");
112 }
113 phys_bc_type[bcvar_type][ori] = REMORA_BC::periodic;
114 phys_bc_need_data[bdy_index[bcvar_type]][ori] = false;
115 domain_bc_type[ori] = "Periodic";
116 }
117 else
118 {
119 phys_bc_type[bcvar_type][ori] = REMORA_BC::undefined;
120 phys_bc_need_data[bdy_index[bcvar_type]][ori] = false;
121 }
122
123 if (geom[0].isPeriodic(ori.coordDir()))
124 {
125 domain_bc_type[ori] = "Periodic";
126 if (phys_bc_type[bcvar_type][ori] == REMORA_BC::undefined)
127 {
128 phys_bc_type[bcvar_type][ori] = REMORA_BC::periodic;
129 phys_bc_need_data[bdy_index[bcvar_type]][ori] = false;
130 } else if (phys_bc_type[bcvar_type][ori] != REMORA_BC::periodic) {
131 amrex::Abort("Wrong BC type for periodic boundary");
132 }
133 }
134
135 if (phys_bc_type[bcvar_type][ori] == REMORA_BC::undefined)
136 {
137 amrex::Print() << "BC Type specified for fac is " << bc_type_string << std::endl;
138 amrex::Abort("This BC type is unknown");
139 }
140 };
141
142 auto f_by_side = [this, &f_set_var_bc] (std::string const& bcid, Orientation ori)
143 {
144 ParmParse pp("remora.bc."+bcid);
145 std::string bc_type_in = "null";
146 // Default z directions to slipwall
147 if (bcid=="zlo" or bcid=="zhi") bc_type_in = "slipwall";
148 pp.query("type", bc_type_in);
149 std::string bc_type = amrex::toLower(bc_type_in);
150
151 for (int icomp=0; icomp<BCVars::NumTypes; icomp++) {
152 f_set_var_bc(pp, icomp, ori, bc_type);
153 }
154 };
155
156 auto f_by_var = [this, &f_set_var_bc] (std::string const& varname, int bcvar_type)
157 {
158 amrex::Vector<Orientation> orientations = {Orientation(Direction::x,Orientation::low), Orientation(Direction::y,Orientation::high),Orientation(Direction::x,Orientation::high),Orientation(Direction::y,Orientation::low)}; // west, south, east, north [matches ROMS]
159 std::vector<std::string> bc_types = {"null","null","null","null"};
160 ParmParse pp("remora.bc."+varname);
161 std::string bc_type_in = "null";
162 // default zvel to outflow
163 if (bcvar_type == BCVars::zvel_bc) {
164 bc_types = {"outflow","outflow","outflow","outflow"};
165 for (int i=0; i<4; i++) {
166 auto ori = orientations[i];
167 if (geom[0].isPeriodic(ori.coordDir())) {
168 bc_types[i] = "periodic";
169 }
170 }
171 }
172 pp.queryarr("type", bc_types);
173 AMREX_ASSERT(bc_types.size() == 4);
174 for (int i=0; i<4; i++) {
175 std::string bc_type = amrex::toLower(bc_types[i]);
176 auto ori = orientations[i];
177 f_set_var_bc(pp, bcvar_type, ori, bc_type);
178 }
179 };
180
189
190 for (OrientationIter oit; oit; ++oit) {
191 Orientation ori = oit();
192 // These are simply defaults for Dirichlet faces -- they should be over-written below if needed
193 m_bc_extdir_vals[BCVars::Temp_bc_comp ][ori] = 1.e19_rt;
194 m_bc_extdir_vals[BCVars::Salt_bc_comp ][ori] = 1.e20_rt;
196
197 m_bc_extdir_vals[BCVars::xvel_bc][ori] = 0.0_rt; // default
198 m_bc_extdir_vals[BCVars::yvel_bc][ori] = 0.0_rt;
199 m_bc_extdir_vals[BCVars::zvel_bc][ori] = 0.0_rt;
200
201 m_bc_extdir_vals[BCVars::ubar_bc][ori] = 0.0_rt; // default
202 m_bc_extdir_vals[BCVars::vbar_bc][ori] = 0.0_rt;
205 }
206
207 // Whether to specify boundary conditions by variable (then side).
208 // Alternative is to do it by side by indicating keywords that indicate multiple variables
209 set_bcs_by_var = false;
210
211 ParmParse pp("remora");
212 pp.query("boundary_per_variable", set_bcs_by_var);
213 if (!set_bcs_by_var) {
214 f_by_side("xlo", Orientation(Direction::x,Orientation::low));
215 f_by_side("xhi", Orientation(Direction::x,Orientation::high));
216 f_by_side("ylo", Orientation(Direction::y,Orientation::low));
217 f_by_side("yhi", Orientation(Direction::y,Orientation::high));
218 } else {
219 f_by_var("temp", BCVars::Temp_bc_comp);
220 f_by_var("salt", BCVars::Salt_bc_comp);
221 f_by_var("scalar", BCVars::Scalar_bc_comp);
222 f_by_var("u", BCVars::xvel_bc);
223 f_by_var("v", BCVars::yvel_bc);
224 f_by_var("w", BCVars::zvel_bc);
225 f_by_var("ubar", BCVars::ubar_bc);
226 f_by_var("vbar", BCVars::vbar_bc);
227 f_by_var("zeta", BCVars::zeta_bc);
228 f_by_var("tke", BCVars::tke_bc);
229 }
230
231 // Always specify z direction by side keyword
232 f_by_side("zlo", Orientation(Direction::z,Orientation::low));
233 f_by_side("zhi", Orientation(Direction::z,Orientation::high));
234
235 // *****************************************************************************
236 //
237 // Here we translate the physical boundary conditions -- one type per face --
238 // into logical boundary conditions for each velocity component
239 //
240 // *****************************************************************************
241 {
242 domain_bcs_type.resize(AMREX_SPACEDIM+NCONS+8);
243 domain_bcs_type_d.resize(AMREX_SPACEDIM+NCONS+8);
244
245 for (OrientationIter oit; oit; ++oit) {
246 Orientation ori = oit();
247 int dir = ori.coordDir();
248 Orientation::Side side = ori.faceDir();
249 // only do this for xvel and yvel
250 for (int i = 0; i < 2; i++) {
251 auto const bct = phys_bc_type[BCVars::xvel_bc+i][ori];
252 if ( bct == REMORA_BC::symmetry )
253 {
254 if (side == Orientation::low) {
256 if (i==1)
258 } else {
260 if (i==1)
262 }
263 }
264 else if (bct == REMORA_BC::outflow)
265 {
266 if (side == Orientation::low) {
268 } else {
270 }
271 }
272 else if (bct == REMORA_BC::inflow)
273 {
274 if (side == Orientation::low) {
276 } else {
278 }
279 }
280 else if (bct == REMORA_BC::no_slip_wall)
281 {
282 if (side == Orientation::low) {
284 } else {
286 }
287 }
288 else if (bct == REMORA_BC::slip_wall)
289 {
290 if (side == Orientation::low) {
292 if (i==1) {
293 // Only normal direction has ext_dir
295 }
296
297 } else {
299 if (i==1) {
300 // Only normal direction has ext_dir
302 }
303 }
304 }
305 else if (bct == REMORA_BC::periodic)
306 {
307 if (side == Orientation::low) {
309 } else {
311 }
312 }
313 else if (bct == REMORA_BC::clamped)
314 {
315 if (side == Orientation::low) {
317 } else {
319 }
320 }
321 else if (bct == REMORA_BC::orlanski_rad)
322 {
323 if (side == Orientation::low) {
325 } else {
327 }
328 }
329 else if (bct == REMORA_BC::orlanski_rad_nudge)
330 {
331 if (side == Orientation::low) {
333 } else {
335 }
336 }
337 else
338 {
339 amrex::Abort("Velocity boundary condition not validly specified");
340 }
341 }
342
343 // Always set zvel_bc to foextrap
344 if (side == Orientation::low) {
346 } else {
348 }
349 }
350 }
351
352 // *****************************************************************************
353 //
354 // Here we translate the physical boundary conditions -- one type per face --
355 // into logical boundary conditions for each cell-centered variable
356 //
357 // *****************************************************************************
358 {
359 for (OrientationIter oit; oit; ++oit) {
360 Orientation ori = oit();
361 int dir = ori.coordDir();
362 Orientation::Side side = ori.faceDir();
363 for (int i = 0; i < NCONS; i++) {
364 auto const bct = phys_bc_type[BCVars::cons_bc+i][ori];
365 if ( bct == REMORA_BC::symmetry )
366 {
367 if (side == Orientation::low) {
369 } else {
371 }
372 }
373 else if ( bct == REMORA_BC::outflow)
374 {
375 if (side == Orientation::low) {
377 } else {
379 }
380 }
381 else if ( bct == REMORA_BC::no_slip_wall)
382 {
383 if (side == Orientation::low) {
385 } else {
387 }
388 }
389 else if (bct == REMORA_BC::slip_wall)
390 {
391 if (side == Orientation::low) {
393 } else {
395 }
396 }
397 else if (bct == REMORA_BC::inflow)
398 {
399 if (side == Orientation::low) {
401 } else {
403 }
404 }
405 else if (bct == REMORA_BC::periodic)
406 {
407 if (side == Orientation::low) {
409 } else {
411 }
412 }
413 else if ( bct == REMORA_BC::clamped)
414 {
415 if (side == Orientation::low) {
417 } else {
419 }
420 }
421 else if ( bct == REMORA_BC::orlanski_rad)
422 {
423 if (side == Orientation::low) {
425 } else {
427 }
428 }
429 else if ( bct == REMORA_BC::orlanski_rad_nudge)
430 {
431 if (side == Orientation::low) {
433 } else {
435 }
436 }
437 else
438 {
439 amrex::Abort("Scalar/tracer boundary condition not validly specified");
440 }
441 }
442 }
443 }
444
445 // *****************************************************************************
446 //
447 // Here we translate the physical boundary conditions -- one type per face --
448 // into logical boundary conditions for ubar and vbar. Also add simplified
449 // 2d boundary condition (corresponding to BCs in bc_2d.F
450 //
451 // *****************************************************************************
452 {
453 for (OrientationIter oit; oit; ++oit) {
454 Orientation ori = oit();
455 int dir = ori.coordDir();
456 Orientation::Side side = ori.faceDir();
457 for (int i = 0; i < 2; i++) {
458 auto const bct = phys_bc_type[BCVars::ubar_bc+i][ori];
459 if ( bct == REMORA_BC::symmetry )
460 {
461 if (side == Orientation::low) {
464 if (i==1 and dir!=2) {
467 }
468 } else {
471 if (i==1 and dir!=2) {
474 }
475 }
476 }
477 else if (bct == REMORA_BC::outflow)
478 {
479 if (side == Orientation::low) {
482 } else {
485 }
486 }
487 else if (bct == REMORA_BC::inflow)
488 {
489 if (side == Orientation::low) {
492 } else {
495 }
496 }
497 else if (bct == REMORA_BC::no_slip_wall)
498 {
499 if (side == Orientation::low) {
502 } else {
505 }
506 }
507 else if (bct == REMORA_BC::slip_wall)
508 {
509 if (side == Orientation::low) {
512 if (i==1 and dir!=2) {
513 // Only normal direction has ext_dir
516 }
517
518 } else {
521 if (i==1 and dir!=2) {
522 // Only normal direction has ext_dir
525 }
526 }
527 }
528 else if (bct == REMORA_BC::periodic)
529 {
530 if (side == Orientation::low) {
533 } else {
536 }
537 }
538 else if (bct == REMORA_BC::clamped)
539 {
540 if (side == Orientation::low) {
543 } else {
546 }
547 }
548 else if (bct == REMORA_BC::flather)
549 {
550 if (side == Orientation::low) {
553 if (i==1 and dir!=2) {
554 // Only normal direction has Flather
557 }
558
559 } else {
562 if (i==1 and dir!=2) {
563 // Only normal direction has Flather
566 }
567 }
568 }
569 else
570 {
571 amrex::Abort("ubar or vbar boundary condition not validly specified");
572 }
573 }
574 }
575 }
576
577 // *****************************************************************************
578 //
579 // Here we translate the physical boundary conditions -- one type per face --
580 // into logical boundary conditions for zeta and tke
581 //
582 // *****************************************************************************
583 {
584 for (OrientationIter oit; oit; ++oit) {
585 Orientation ori = oit();
586 int dir = ori.coordDir();
587 Orientation::Side side = ori.faceDir();
588 for (int i = 0; i < 2; i++) {
589 auto const bct = phys_bc_type[BCVars::zeta_bc+i][ori];
590 if ( bct == REMORA_BC::symmetry )
591 {
592 if (side == Orientation::low) {
594 } else {
596 }
597 }
598 else if ( bct == REMORA_BC::outflow)
599 {
600 if (side == Orientation::low) {
602 } else {
604 }
605 }
606 else if ( bct == REMORA_BC::no_slip_wall)
607 {
608 if (side == Orientation::low) {
610 } else {
612 }
613 }
614 else if (bct == REMORA_BC::slip_wall)
615 {
616 if (side == Orientation::low) {
618 } else {
620 }
621 }
622 else if (bct == REMORA_BC::inflow)
623 {
624 if (side == Orientation::low) {
626 } else {
628 }
629 }
630 else if (bct == REMORA_BC::periodic)
631 {
632 if (side == Orientation::low) {
634 } else {
636 }
637 }
638 else if (bct == REMORA_BC::chapman)
639 {
640 if (side == Orientation::low) {
642 } else {
644 }
645 }
646 else if ( bct == REMORA_BC::clamped)
647 {
648 if (side == Orientation::low) {
650 } else {
652 }
653 }
654 else
655 {
656 amrex::Abort("Free surface (zeta) boundary condition not validly specified");
657 }
658 }
659 }
660 }
661
662 // *****************************************************************************
663 //
664 // Here we define a boundary condition that will foextrap while respecting periodicity
665 // This is used as a "null BC"
666 //
667 // *****************************************************************************
668 {
669 for (OrientationIter oit; oit; ++oit) {
670 Orientation ori = oit();
671 int dir = ori.coordDir();
672 Orientation::Side side = ori.faceDir();
673 if (side == Orientation::low) {
675 } else {
677 }
678 }
679 }
680
681 // *****************************************************************************
682 //
683 // Here we define a boundary condition that will unconditionally foextrap
684 //
685 // *****************************************************************************
686 {
687 for (OrientationIter oit; oit; ++oit) {
688 Orientation ori = oit();
689 int dir = ori.coordDir();
690 Orientation::Side side = ori.faceDir();
691 if (side == Orientation::low) {
693 } else {
695 }
696 }
697 }
698
699
700#ifdef AMREX_USE_GPU
701 Gpu::htod_memcpy
702 (domain_bcs_type_d.data(), domain_bcs_type.data(),
703 sizeof(amrex::BCRec)*(NCONS+AMREX_SPACEDIM+8));
704#else
705 std::memcpy
706 (domain_bcs_type_d.data(), domain_bcs_type.data(),
707 sizeof(amrex::BCRec)*(NCONS+AMREX_SPACEDIM+8));
708#endif
709}
710
@ orlanski_rad_nudge
#define NCONS
void init_bcs()
Read in boundary parameters from input file and set up data structures.
amrex::Vector< amrex::BCRec > domain_bcs_type
vector (over BCVars) of BCRecs
Definition REMORA.H:1136
bool set_bcs_by_var
whether to set boundary conditions by variable rather than just by side
Definition REMORA.H:1121
amrex::GpuArray< amrex::GpuArray< REMORA_BC, AMREX_SPACEDIM *2 >, BCVars::NumTypes > phys_bc_type
Array holding the "physical" boundary condition types (e.g. "inflow")
Definition REMORA.H:1147
amrex::Vector< int > bdy_index
Container to connect boundary data being read in boundary condition containers.
Definition REMORA.H:1156
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1237
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NCONS+8 > m_bc_extdir_vals
Array holding the Dirichlet values at walls which need them.
Definition REMORA.H:1144
amrex::Array< std::string, 2 *AMREX_SPACEDIM > domain_bc_type
Array of strings describing domain boundary conditions.
Definition REMORA.H:1141
amrex::GpuArray< amrex::GpuArray< bool, AMREX_SPACEDIM *2 >, BdyVars::NumTypes+1 > phys_bc_need_data
These are flags that indicate whether we need to read in boundary data from file.
Definition REMORA.H:1150
amrex::Gpu::DeviceVector< amrex::BCRec > domain_bcs_type_d
GPU vector (over BCVars) of BCRecs.
Definition REMORA.H:1138