Skip to content

Commit ff42ec6

Browse files
RevathiJambunathanRevathiJambunathan
committed
NodalPoissonSolvewithEB
1 parent 4c7a49d commit ff42ec6

File tree

7 files changed

+225
-0
lines changed

7 files changed

+225
-0
lines changed

‎ExampleCodes/EB/NodalPoisson/CMakeLists.txt‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
set(_sources main.cpp Poisson.cpp Poisson.H)
2+
set(_input_files inputs)
3+
4+
setup_tutorial(_sources _input_files)
5+
6+
unset(_sources)
7+
unset(_input_files)

‎ExampleCodes/EB/NodalPoisson/GNUmakefile‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
USE_EB = TRUE
2+
DEBUG = FALSE
3+
USE_MPI = TRUE
4+
USE_OMP = FALSE
5+
6+
USE_HYPRE = FALSE
7+
8+
COMP = gnu
9+
10+
DIM = 2
11+
12+
AMREX_HOME ?= ../../../../amrex_upstream
13+
14+
include $(AMREX_HOME)/Tools/GNUMake/Make.defs
15+
16+
include ./Make.package
17+
18+
Pdirs := Base Boundary AmrCore
19+
Pdirs += EB
20+
Pdirs += LinearSolvers/MLMG
21+
22+
Ppack += $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir)/Make.package)
23+
24+
include $(Ppack)
25+
26+
include $(AMREX_HOME)/Tools/GNUMake/Make.rules
27+

‎ExampleCodes/EB/NodalPoisson/Make.package‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
2+
CEXE_sources += main.cpp
3+
4+
CEXE_headers += Poisson.H
5+
CEXE_sources += Poisson.cpp

‎ExampleCodes/EB/NodalPoisson/Poisson.H‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
#ifndef POISSON_H_
2+
#define POISSON_H_
3+
4+
#include <AMReX_EBMultiFabUtil.H>
5+
6+
void InitData (amrex::MultiFab& State);
7+
8+
#endif

‎ExampleCodes/EB/NodalPoisson/Poisson.cpp‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
#include "Poisson.H"
2+
3+
using namespace amrex;
4+
5+
void InitData (MultiFab& State)
6+
{
7+
#ifdef AMREX_USE_OMP
8+
#pragma omp parallel
9+
#endif
10+
for (MFIter mfi(State,true); mfi.isValid(); ++mfi)
11+
{
12+
const Box& bx = mfi.growntilebox();
13+
const Array4<Real>& q = State.array(mfi);
14+
amrex::ParallelFor(bx,
15+
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
16+
{
17+
if (i==70 && j==70 && k==0) {
18+
q(i,j,k) = 1.0;
19+
} else {
20+
q(i,j,k) = 0.0;
21+
}
22+
});
23+
}
24+
}

‎ExampleCodes/EB/NodalPoisson/inputs‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
amrex.fpe_trap_invalid=1
2+
3+
verbose = 2
4+
n_cell = 128
5+
max_grid_size = 32
6+
7+
eb2.geom_type = sphere
8+
eb2.sphere_center = 0.5 0.5 0.5
9+
eb2.sphere_radius = 0.25
10+
eb2.sphere_has_fluid_inside = 1
11+

‎ExampleCodes/EB/NodalPoisson/main.cpp‎

Lines changed: 143 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,143 @@
1+
2+
#include <AMReX.H>
3+
#include <AMReX_ParmParse.H>
4+
#include <AMReX_EBMultiFabUtil.H>
5+
#include <AMReX_EB2.H>
6+
#include <AMReX_EB2_IF.H>
7+
#include <AMReX_MLEBABecLap.H>
8+
#include <AMReX_MLEBNodeFDLaplacian.H>
9+
#include <AMReX_PlotFileUtil.H>
10+
#include <AMReX_MultiFabUtil.H>
11+
#include <AMReX_MLMG.H>
12+
13+
#include "Poisson.H"
14+
15+
using namespace amrex;
16+
17+
int main (int argc, char* argv[])
18+
{
19+
amrex::Initialize(argc, argv);
20+
21+
{
22+
int verbose = 1;
23+
int n_cell = 128;
24+
int max_grid_size = 32;
25+
26+
// read parameters
27+
{
28+
ParmParse pp;
29+
pp.query("verbose", verbose);
30+
pp.query("n_cell", n_cell);
31+
pp.query("max_grid_size", max_grid_size);
32+
}
33+
34+
Geometry geom;
35+
BoxArray grids;
36+
DistributionMapping dmap;
37+
{
38+
RealBox rb({AMREX_D_DECL(0.,0.,0.)}, {AMREX_D_DECL(1.,1.,1.)});
39+
Array<int,AMREX_SPACEDIM> is_periodic{AMREX_D_DECL(1,1,1)};
40+
Box domain(IntVect{AMREX_D_DECL(0,0,0)},
41+
IntVect{AMREX_D_DECL(n_cell-1,n_cell-1,n_cell-1)});
42+
geom.define(domain, rb, CoordSys::cartesian, is_periodic);
43+
44+
grids.define(domain); // define the BoxArray to be a single grid
45+
grids.maxSize(max_grid_size); // chop domain up into boxes with length max_Grid_size
46+
47+
dmap.define(grids); // create a processor distribution mapping given the BoxARray
48+
}
49+
50+
int required_coarsening_level = 0; // typically the same as the max AMR level index
51+
int max_coarsening_level = 100; // typically a huge number so MG coarsens as much as possible
52+
// build a simple geometry using the "eb2." parameters in the inputs file
53+
EB2::Build(geom, required_coarsening_level, max_coarsening_level);
54+
55+
const EB2::IndexSpace& eb_is = EB2::IndexSpace::top();
56+
const EB2::Level& eb_level = eb_is.getLevel(geom);
57+
58+
// options are basic, volume, or full
59+
EBSupport ebs = EBSupport::full;
60+
61+
// number of ghost cells for each of the 3 EBSupport types
62+
Vector<int> ng_ebs = {2,2,2};
63+
64+
// This object provides access to the EB database in the format of basic AMReX objects
65+
// such as BaseFab, FArrayBox, FabArray, and MultiFab
66+
EBFArrayBoxFactory factory(eb_level, geom, grids, dmap, ng_ebs, ebs);
67+
68+
// charge density and electric potential are nodal
69+
const BoxArray& nba = amrex::convert(grids, IntVect::TheNodeVector());
70+
MultiFab q (nba, dmap, 1, 0, MFInfo(), factory);
71+
MultiFab phi(nba, dmap, 1, 0, MFInfo(), factory);
72+
73+
InitData(q);
74+
75+
LPInfo info;
76+
77+
// MLEBABecLap mlebabec({geom},{grids},{dmap},info,{&factory});
78+
MLEBNodeFDLaplacian linop({geom},{grids},{dmap},info,{&factory});
79+
80+
// define array of LinOpBCType for domain boundary conditions
81+
std::array<LinOpBCType,AMREX_SPACEDIM> bc_lo;
82+
std::array<LinOpBCType,AMREX_SPACEDIM> bc_hi;
83+
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
84+
bc_lo[idim] = LinOpBCType::Periodic;
85+
bc_hi[idim] = LinOpBCType::Periodic;
86+
}
87+
88+
// Boundary of the whole domain. This functions must be called,
89+
// and must be called before other bc functions.
90+
linop.setDomainBC(bc_lo,bc_hi);
91+
92+
// see AMReX_MLLinOp.H for an explanation
93+
linop.setLevelBC(0, nullptr);
94+
95+
//// operator looks like (ACoef - div BCoef grad) phi = rhs
96+
97+
//// set ACoef to zero
98+
//MultiFab acoef(grids, dmap, 1, 0, MFInfo(), factory);
99+
//acoef.setVal(0.);
100+
//mlebabec.setACoeffs(0, acoef);
101+
102+
//// set BCoef to 1.0 (and array of face-centered coefficients)
103+
//Array<MultiFab,AMREX_SPACEDIM> bcoef;
104+
//for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
105+
// bcoef[idim].define(amrex::convert(grids,IntVect::TheDimensionVector(idim)), dmap, 1, 0, MFInfo(), factory);
106+
// bcoef[idim].setVal(1.0);
107+
//}
108+
//mlebabec.setBCoeffs(0, amrex::GetArrOfConstPtrs(bcoef));
109+
110+
//// scaling factors; these multiply ACoef and BCoef
111+
//Real ascalar = 0.0;
112+
//Real bscalar = 1.0;
113+
//mlebabec.setScalars(ascalar, bscalar);
114+
115+
//// think of this beta as the "BCoef" associated with an EB face
116+
//MultiFab beta(nba, dmap, 1, 0, MFInfo(), factory);
117+
//beta.setVal(1.);
118+
119+
// set homogeneous Dirichlet BC for EB
120+
linop.setEBDirichlet(0.);
121+
122+
MLMG mlmg(linop);
123+
124+
// relative and absolute tolerances for linear solve
125+
const Real tol_rel = 1.e-10;
126+
const Real tol_abs = 0.0;
127+
128+
mlmg.setVerbose(verbose);
129+
130+
// Solve linear system
131+
phi.setVal(0.0); // initial guess for phi
132+
mlmg.solve({&phi}, {&q}, tol_rel, tol_abs);
133+
134+
//// store plotfile variables; q and phi
135+
//MultiFab plotfile_mf(grids, dmap, 2, 0, MFInfo(), factory);
136+
//MultiFab::Copy(plotfile_mf, q,0,0,1,0);
137+
//MultiFab::Copy(plotfile_mf,phi,0,1,1,0);
138+
139+
//EB_WriteSingleLevelPlotfile("plt", plotfile_mf, {"q", "phi"}, geom, 0.0, 0);
140+
}
141+
142+
amrex::Finalize();
143+
}

0 commit comments

Comments
 (0)