Boundary conditions
Mattiads1 opened this issue · 12 comments
Is it possible to use this library with mixed BC, i.e. Dirichlet on one side of my domain and Neumann on the other side?
Thanks a lot
Hello, yes, it is possible. Only certain combinations are currently supported, though. Try it in the dirichlet-solvers branch and let me know if your combination is missing.
Also is it possible into a 3D domain to use periodic BC in two directions and a mixed BC (Dirichelet and Neumann) on the other direction?
Yes, try the dirichlet-solvers branch.
Hello. I am wondering if the mixed BCs combination (Periodic for x_min x_max, Neumann for y_min y_max) for a 2D domain can be used. I tried the following problem designed by myself in the testpoisson.cc
rhs[IND(i,j)] = cos(2*pi*x/Ls[0]) * cos(5*pi*y/Ls[1])
Exact_solution[IND(i,j)] = - 1.0 / ( pow(2*pi/Ls[0], 2.0) + pow(5*pi/Ls[1], 2.0) ) * cos(2*pi*x/Ls[0]) * cos(5*pi*y/Ls[1]) *
,but I failed.
Here is my complete test:
#include <cmath>
#include <iostream>
#include "poisfft.h"
#include <stdlib.h>
const double pi = 3.14159265358979323846;
#define IND(i,j) i*ns[1]+j
void init_rhs(const int ns[2], const double ds[2], const double Ls[2], double* a){
int i,j,k;
for (i=0;i<ns[0];i++){
for (j=0;j<ns[1];j++){
double x = ds[0]*(i+0.5);
double y = ds[1]*(j+0.5);
a[IND(i,j)] = cos(2*pi*x/Ls[0]) * cos(5*pi*y/Ls[1]);
}
}
}
void check_solution(const int ns[2], const double ds[2], const double Ls[2], double* a){
int i,j,k;
double sum = 0;
for (i=0;i<ns[0];i++){
for (j=0;j<ns[1];j++){
double x = ds[0]*(i+0.5);
double y = ds[1]*(j+0.5);
sum += pow(a[IND(i,j)] -
(- 1.0 / ( pow(2*pi/Ls[0], 2.0) +
pow(5*pi/Ls[1], 2.0) ) *
cos(2*pi*x/Ls[0]) *
cos(5*pi*y/Ls[1]))
, 2.0);
}
}
if (sum < 1e-10) {
std::cout << "OK\n" << std::endl;
}
else {
std::cout << "FAIL, residuum" << sum << std::endl;
}
}
int main()
{
// domain dimensions
const double Ls[2] = {1.1, 2.0};
// gridpoint numbers
const int ns[2] = {71, 93};
// distances between gridpoints
double ds[2];
//boundary conditions
const int BCs[4] = {PoisFFT::PERIODIC, PoisFFT::PERIODIC,
PoisFFT::NEUMANN_STAG, PoisFFT::NEUMANN_STAG};
int i;
// set the grid, depends on the boundary conditions
for (i = 0; i<2; i++){
ds[i] = Ls[i] / ns[i];
}
// allocate the arrays contiguously, you can use any other class
// from which you can get a pointer to contiguous buffer
double *arr = new double[ns[0]*ns[1]];
double *RHS = new double[ns[0]*ns[1]];
// set the right-hand side
init_rhs(ns, ds, Ls, RHS);
// create solver object, 3 dimensions, double precision
PoisFFT::Solver<2, double> S(ns, Ls, BCs);
//run the solver, can be run many times for different right-hand sides
S.execute(arr, RHS);
// check corectness
check_solution(ns, ds, Ls, arr);
delete[] RHS;
delete[] arr;
return 0;
}
It should be possible. Could you show a complete test example?
…
On Wed, 29 Jul 2020 at 11:35, baliangLt @.> wrote: Hello. I am wondering if the mixed BCs combination (Periodic for x_min x_max, Neumann for y_min y_max) for a 2D domain can be used. I tried the following problem designed by myself in the testpoisson.cc rhs[IND(i,j)] = cos(2pix/Ls[0]) * cos(5piy/Ls[1]) Exact_solution[IND(i,j)] = - 1.0 / ( pow(2pi/Ls[0], 2.0) + pow(5pi/Ls[1], 2.0) ) * cos(2pix/Ls[0]) * cos(5pi*y/Ls[1]) * ,but I failed. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#9 (comment)>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAFSIENBMDAYC2QHA3CKDDTR57UPDANCNFSM4KMT3PRA .
Something is wrong, I will have a look.
Something is wrong, I will have a look.
ok, thanks~
OK, I was wrong, it is not actually supported (yet). After seeing there is interest I will try to find time to write it.
As a workaround, you could call the 3D solver with a dummy 3D dimension of length 1, it should not be that much slower. See the code below
#include <cmath>
#include <iostream>
#include "poisfft.h"
#include <stdlib.h>
const double pi = 3.14159265358979323846;
#define IND(i,j) (i)*ns[1]+j
void init_rhs(const int ns[2], const double ds[2], const double Ls[2], double* a){
int i,j,k;
for (i=0;i<ns[0];i++){
for (j=0;j<ns[1];j++){
double x = ds[0]*(i+0.5);
double y = ds[1]*(j+0.5);
a[IND(i,j)] = cos(2*pi*x/Ls[0]) * cos(5*pi*y/Ls[1]);
}
}
}
void check_solution(const int ns[2], const double ds[2], const double Ls[2], double* a){
int i,j,k;
double sum = 0;
for (i=0;i<ns[0];i++){
for (j=0;j<ns[1];j++){
double x = ds[0]*(i+0.5);
double y = ds[1]*(j+0.5);
sum += pow(a[IND(i,j)] -
(- 1.0 / ( pow(2*pi/Ls[0], 2.0) +
pow(5*pi/Ls[1], 2.0) ) *
cos(2*pi*x/Ls[0]) *
cos(5*pi*y/Ls[1]))
, 2.0);
}
}
if (sum < 1e-10) {
std::cout << "OK\n" << std::endl;
}
else {
std::cout << "FAIL, residuum" << sum << std::endl;
}
}
int main()
{
// domain dimensions
const double Ls[2] = {1.1, 2.0};
const double Ls3D[3] = {1.1, 2.0, 1.0};
// gridpoint numbers
const int ns[2] = {71, 93};
const int ns3D[3] = {71, 93,1};
// distances between gridpoints
double ds[2];
//boundary conditions
const int BCs[6] = {PoisFFT::PERIODIC, PoisFFT::PERIODIC,
PoisFFT::NEUMANN_STAG, PoisFFT::NEUMANN_STAG,
PoisFFT::NEUMANN_STAG, PoisFFT::NEUMANN_STAG};
int i;
// set the grid, depends on the boundary conditions
for (i = 0; i<2; i++){
ds[i] = Ls[i] / ns[i];
}
// allocate the arrays contiguously, you can use any other class
// from which you can get a pointer to contiguous buffer
double *arr = new double[ns[0]*ns[1]];
double *RHS = new double[ns[0]*ns[1]];
// set the right-hand side
init_rhs(ns, ds, Ls, RHS);
// create solver object, 3 dimensions, double precision
PoisFFT::Solver<3, double> S(ns3D, Ls3D, BCs);
//run the solver, can be run many times for different right-hand sides
S.execute(arr, RHS);
// check corectness
check_solution(ns, ds, Ls, arr);
delete[] RHS;
delete[] arr;
return 0;
}