Commit 4802e0d1 authored by Abhinav Singh's avatar Abhinav Singh
Browse files

Adding examples to develop

parent e4ce975f
include ../../example.mk
CC=mpic++
LDIR =
OBJ = main.o
%.o: %.cpp
$(CC) -O0 -g -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
Active2d: $(OBJ)
$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
all: Active2d
run: all
mpirun -np 4 /Active2d
.PHONY: clean all run
clean:
rm -f *.o *~ core Active2d
This diff is collapsed.
include ../../../example.mk
CC=mpic++
LDIR =
OBJ = mainDCPSE.o
OBJ2 = mainFD.o
%.o: %.cpp
$(CC) -O3 -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
DcpseLid: $(OBJ)
$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
FDLid: $(OBJ2)
$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
all: DcpseLid FDLid
run: all
mpirun -np 1 /DcpseLid
.PHONY: clean all run
clean:
rm -f *.o *~ core DcpseLid FDLid
[pack]
files = mainDCPSE.cpp mainFD.cpp Makefile
//
// Created by Abhinav Singh on 15.11.2021.
//
#include "config.h"
#include <iostream>
#include "DCPSE/DCPSE_OP/DCPSE_op.hpp"
#include "DCPSE/DCPSE_OP/DCPSE_Solver.hpp"
#include "Operators/Vector/vector_dist_operators.hpp"
#include "Vector/vector_dist_subset.hpp"
int main(int argc, char* argv[])
{
{ openfpm_init(&argc,&argv);
timer tt2;
tt2.start();
size_t gd_sz = 81;
double Re=100;
double V_err_eps = 1e-2;
double alpha=0.0125;
constexpr int x = 0;
constexpr int y = 1;
const size_t sz[2] = {gd_sz,gd_sz};
Box<2, double> box({0, 0}, {1,1});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing;
spacing = 1.0 / (sz[0] - 1);
double rCut = 3.1 * spacing;
int ord = 2;
Ghost<2, double> ghost(rCut);
auto &v_cl = create_vcluster();
typedef aggregate<double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double,double> LidCavity;
vector_dist_ws<2, double, LidCavity> Particles(0, box,bc,ghost);
Particles.setPropNames({"00-Pressure","01-Velocity","02-RHS","03-dV","04-Divergence","05-Divergence","06-H","07-dP"});
double x0, y0, x1, y1;
x0 = box.getLow(0);
y0 = box.getLow(1);
x1 = box.getHigh(0);
y1 = box.getHigh(1);
auto it = Particles.getGridIterator(sz);
while (it.isNext())
{
Particles.add();
auto key = it.get();
mem_id k0 = key.get(0);
double xp = k0 * spacing;
Particles.getLastPos()[0] = xp;
mem_id k1 = key.get(1);
double yp = k1 * spacing;
Particles.getLastPos()[1] = yp;
Particles.getLastProp<1>()[0]=0.0;
Particles.getLastProp<1>()[1]=0.0;
if (xp != x0 && yp != y0 && xp != x1 && yp != y1)
{
Particles.getLastSubset(0);
}
else if(yp==y1 && xp>x0 && xp<x1)
{
Particles.getLastSubset(1);
Particles.getLastProp<1>()[0]=1.0;
}
else if(xp==0)
{
Particles.getLastSubset(3);
}
else if(xp==x1)
{
Particles.getLastSubset(4);
}
else
{
Particles.getLastSubset(2);
}
++it;
}
Particles.map();
Particles.ghost_get<0>();
vector_dist_subset<2, double, LidCavity> Particles_bulk(Particles,0);
vector_dist_subset<2, double, LidCavity> Particles_up(Particles,1);
vector_dist_subset<2, double, LidCavity> Particles_down(Particles,2);
vector_dist_subset<2, double, LidCavity> Particles_left(Particles,3);
vector_dist_subset<2, double, LidCavity> Particles_right(Particles,4);
auto &bulk=Particles_bulk.getIds();
auto &up_p=Particles_up.getIds();
auto &dw_p=Particles_down.getIds();
auto &l_p=Particles_left.getIds();
auto &r_p=Particles_right.getIds();
auto P = getV<0>(Particles);
auto V = getV<1>(Particles);
auto RHS = getV<2>(Particles);
auto dV = getV<3>(Particles);
auto div = getV<4>(Particles);
auto V_star = getV<5>(Particles);
auto H = getV<6>(Particles);
auto dP = getV<7>(Particles);
auto P_bulk = getV<0>(Particles_bulk);
auto V_bulk = getV<1>(Particles_bulk);
auto RHS_bulk =getV<2>(Particles_bulk);
auto V_star_bulk = getV<5>(Particles_bulk);
auto dP_bulk = getV<7>(Particles_bulk);
P_bulk = 0;
Derivative_x Dx(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
Derivative_xx Dxx(Particles, 2, rCut,sampling_factor2, support_options::RADIUS);
Derivative_yy Dyy(Particles, 2, rCut,sampling_factor2, support_options::RADIUS);
Derivative_y Dy(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
Derivative_x Bulk_Dx(Particles_bulk, 2, rCut,sampling_factor, support_options::RADIUS);
Derivative_y Bulk_Dy(Particles_bulk, 2, rCut,sampling_factor, support_options::RADIUS);
int n = 0, nmax = 300, ctr = 0, errctr=1, Vreset = 0;
double V_err=1;
if (Vreset == 1) {
P_bulk = 0;
P = 0;
Vreset = 0;
}
P=0;
eq_id vx,vy;
vx.setId(0);
vy.setId(1);
/* W.setVarId(0);
Sf.setVarId(1);
*/
double sum, sum1, sum_k,V_err_old;
auto StokesX=(V[x]*Dx(V_star[x])+V[y]*Dy(V_star[x]))-(1.0/Re)*(Dxx(V_star[x])+Dyy(V_star[x]));
auto StokesY=(V[x]*Dx(V_star[y])+V[y]*Dy(V_star[y]))-(1.0/Re)*(Dxx(V_star[y])+Dyy(V_star[y]));
petsc_solver<double> solverPetsc;
solverPetsc.setSolver(KSPGMRES);
RHS[x] = 0.0;
RHS[y] = 0.0;
dV=0;
P=0;
timer tt;
while (V_err >= V_err_eps && n <= nmax) {
if (n%5==0){
Particles.ghost_get<0,1,2,3,4,5,6,7>(SKIP_LABELLING);
Particles.deleteGhost();
Particles.write_frame("LID",n,BINARY);
Particles.ghost_get<0>();
}
tt.start();
Particles.ghost_get<0>(SKIP_LABELLING);
RHS_bulk[x] = -Bulk_Dx(P);
RHS_bulk[y] = -Bulk_Dy(P);
DCPSE_scheme<equations2d2, decltype(Particles)> Solver(Particles);
Solver.impose(StokesX, bulk, RHS[x], vx);
Solver.impose(StokesY, bulk, RHS[y], vy);
Solver.impose(V_star[x], up_p, 1.0, vx);
Solver.impose(V_star[y], up_p, 0.0, vy);
Solver.impose(V_star[x], l_p, 0.0, vx);
Solver.impose(V_star[y], l_p, 0.0, vy);
Solver.impose(V_star[x], r_p, 0.0, vx);
Solver.impose(V_star[y], r_p, 0.0, vy);
Solver.impose(V_star[x], dw_p, 0.0, vx);
Solver.impose(V_star[y], dw_p, 0.0, vy);
Solver.solve_with_solver(solverPetsc,V_star[x], V_star[y]);
Particles.ghost_get<5>(SKIP_LABELLING);
div = (Dx(V_star[x]) + Dy(V_star[y]));
DCPSE_scheme<equations2d1E,decltype(Particles)> SolverH(Particles,options_solver::LAGRANGE_MULTIPLIER);
auto Helmholtz = Dxx(H)+Dyy(H);
SolverH.impose(Helmholtz,bulk,prop_id<4>());
SolverH.impose(Dy(H), up_p,0);
SolverH.impose(Dx(H), l_p,0);
SolverH.impose(Dx(H), r_p,0);
SolverH.impose(Dy(H), dw_p,0);
//SolverH.solve_with_solver(solverPetsc2,H);
SolverH.solve(H);
Particles.ghost_get<6>(SKIP_LABELLING);
Particles.ghost_get<6>(SKIP_LABELLING);
P_bulk = P - alpha*(div-0.5*(V[x]*Bulk_Dx(H)+V[y]*Bulk_Dy(H)));
V_star_bulk[0] = V_star[0] - Bulk_Dx(H);
V_star_bulk[1] = V_star[1] - Bulk_Dy(H);
sum = 0;
sum1 = 0;
for (int j = 0; j < bulk.size(); j++) {
auto p = bulk.get<0>(j);
sum += (Particles.getProp<5>(p)[0] - Particles.getProp<1>(p)[0]) *
(Particles.getProp<5>(p)[0] - Particles.getProp<1>(p)[0]) +
(Particles.getProp<5>(p)[1] - Particles.getProp<1>(p)[1]) *
(Particles.getProp<5>(p)[1] - Particles.getProp<1>(p)[1]);
sum1 += Particles.getProp<5>(p)[0] * Particles.getProp<5>(p)[0] +
Particles.getProp<5>(p)[1] * Particles.getProp<5>(p)[1];
}
V = V_star;
v_cl.sum(sum);
v_cl.sum(sum1);
v_cl.execute();
sum = sqrt(sum);
sum1 = sqrt(sum1);
V_err_old = V_err;
V_err = sum / sum1;
if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-8) {
errctr++;
//alpha_P -= 0.1;
} else {
errctr = 0;
}
if (n > 3) {
if (errctr > 5) {
Vreset = 1;
errctr=0;
//break;
} else {
Vreset = 0;
}
}
n++;
tt.stop();
if (v_cl.rank() == 0) {
std::cout << "Rel l2 cgs err in V = " << V_err << " at " << n << " and took " <<tt.getwct() <<"("<<tt.getcputime()<<") seconds(CPU)." << std::endl;
}
}
Particles.deleteGhost();
Particles.write("LID");
tt2.stop();
if (v_cl.rank() == 0) {
std::cout << "The simulation took " << tt2.getcputime() << "(CPU) ------ " << tt2.getwct()
<< "(Wall) Seconds.";
}
}
openfpm_finalize();
}
\ No newline at end of file
#include "config.h"
#include <iostream>
#include "FiniteDifference/FD_Solver.hpp"
#include "FiniteDifference/FD_op.hpp"
#include "FiniteDifference/FD_expressions.hpp"
#include "FiniteDifference/util/EqnsStructFD.hpp"
int main(int argc, char* argv[])
{
{ openfpm_init(&argc,&argv);
using namespace FD;
timer tt2;
tt2.start();
size_t gd_sz = 81;
constexpr int x = 0;
constexpr int y = 1;
const size_t szu[2] = {gd_sz,gd_sz};
int sz[2]={int(gd_sz),int(gd_sz)};
Box<2, double> box({0, 0}, {1,1});
periodicity<2> bc = {NON_PERIODIC, NON_PERIODIC};
double spacing;
spacing = 1.0 / (sz[0] - 1);
Ghost<2,long int> ghost(1);
auto &v_cl = create_vcluster();
typedef aggregate<double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double,double> LidCavity;
grid_dist_id<2, double, LidCavity> domain(szu, box,ghost,bc);
double x0, y0, x1, y1;
x0 = box.getLow(0);
y0 = box.getLow(1);
x1 = box.getHigh(0);
y1 = box.getHigh(1);
auto it = domain.getDomainIterator();
while (it.isNext())
{
auto key = it.get();
auto gkey = it.getGKey(key);
double x = gkey.get(0) * domain.spacing(0);
double y = gkey.get(1) * domain.spacing(1);
if (y==1)
{domain.get<1>(key) = 1.0;}
else
{domain.get<1>(key) = 0.0;}
++it;
}
domain.ghost_get<0>();
Derivative_x Dx;
Derivative_y Dy;
Derivative_xx Dxx;
Derivative_yy Dyy;
auto P = getV<0>(domain);
auto V = getV<1>(domain);
auto RHS = getV<2>(domain);
auto dV = getV<3>(domain);
auto div = getV<4>(domain);
auto V_star = getV<5>(domain);
auto H = getV<6>(domain);
auto dP = getV<7>(domain);
eq_id x_comp,y_comp;
x_comp.setId(0);
y_comp.setId(1);
double sum, sum1, sum_k,V_err_old;
auto StokesX=(V[x]*Dx(V_star[x])+V[y]*Dy(V_star[x]))-(1.0/100.0)*(Dxx(V_star[x])+Dyy(V_star[x]));
auto StokesY=(V[x]*Dx(V_star[y])+V[y]*Dy(V_star[y]))-(1.0/100.0)*(Dxx(V_star[y])+Dyy(V_star[y]));
RHS=0;
dV=0;
double V_err=1;
int n;
while (V_err >= 0.05 && n <= 50) {
if (n%5==0){
domain.ghost_get<0,1>(SKIP_LABELLING);
domain.write_frame("LID",n,BINARY);
domain.ghost_get<0>();
}
domain.ghost_get<0>(SKIP_LABELLING);
RHS[x] = -Dx(P);
RHS[y] = -Dy(P);
FD_scheme<equations2d2,decltype(domain)> Solver(ghost,domain);
Solver.impose(StokesX, {1,1},{sz[0]-2,sz[1]-2}, RHS[x], x_comp);
Solver.impose(StokesY, {1,1},{sz[0]-2,sz[1]-2}, RHS[y], y_comp);
Solver.impose(V_star[x], {0,sz[1]-1},{sz[0]-1,sz[1]-1}, 1.0, x_comp);
Solver.impose(V_star[y], {0,sz[1]-1},{sz[0]-1,sz[1]-1}, 0.0, y_comp);
Solver.impose(V_star[x], {0,1},{0,sz[1]-2}, 0.0, x_comp);
Solver.impose(V_star[y], {0,1},{0,sz[1]-2}, 0.0, y_comp);
Solver.impose(V_star[x], {sz[0]-1,1},{sz[0]-1,sz[1]-2}, 0.0, x_comp);
Solver.impose(V_star[y], {sz[0]-1,1},{sz[0]-1,sz[1]-2}, 0.0, y_comp);
Solver.impose(V_star[x], {0,0},{sz[0]-1,0}, 0.0, x_comp);
Solver.impose(V_star[y], {0,0},{sz[0]-1,0}, 0.0, y_comp);
Solver.solve(V_star[x], V_star[y]);
domain.ghost_get<5>(SKIP_LABELLING);
div = (Dx(V_star[x]) + Dy(V_star[y]));
FD_scheme<equations2d1E,decltype(domain)> SolverH(ghost,domain);
auto Helmholtz = Dxx(H)+Dyy(H);
SolverH.impose(Helmholtz,{1,1},{sz[0]-2,sz[1]-2},div);
SolverH.impose(Dy(H), {0,sz[0]-1},{sz[0]-1,sz[1]-1},0);
SolverH.impose(Dx(H), {sz[0]-1,1},{sz[0]-1,sz[1]-2},0);
SolverH.impose(Dx(H), {0,0},{sz[0]-1,0},0,x_comp,true);
SolverH.impose(H, {0,0},{0,0},0);
SolverH.impose(Dy(H), {0,1},{0,sz[1]-2},0);
SolverH.solve(H);
domain.ghost_get<1,4,6>(SKIP_LABELLING);
P = P - 0.01*(div-0.5*(V[x]*Dx(H)+V[y]*Dy(H)));
V_star[0] = V_star[0] - Dx(H);
V_star[1] = V_star[1] - Dy(H);
sum = 0;
sum1 = 0;
auto it2 = domain.getDomainIterator();
while (it2.isNext()) {
auto p = it2.get();
sum += (domain.getProp<5>(p)[0] - domain.getProp<1>(p)[0]) *
(domain.getProp<5>(p)[0] - domain.getProp<1>(p)[0]) +
(domain.getProp<5>(p)[1] - domain.getProp<1>(p)[1]) *
(domain.getProp<5>(p)[1] - domain.getProp<1>(p)[1]);
sum1 += domain.getProp<5>(p)[0] * domain.getProp<5>(p)[0] +
domain.getProp<5>(p)[1] * domain.getProp<5>(p)[1];
++it2;
}
V[x] = V_star[x];
V[y] = V_star[y];
v_cl.sum(sum);
v_cl.sum(sum1);
v_cl.execute();
sum = sqrt(sum);
sum1 = sqrt(sum1);
V_err_old = V_err;
V_err = sum / sum1;
n++;
if (v_cl.rank() == 0) {
std::cout << "Rel l2 cgs err in V = " << V_err << "."<< std::endl;
}
}
domain.write("LID_final");
tt2.stop();
if (v_cl.rank() == 0) {
std::cout << "The entire solver took " << tt2.getcputime() << "(CPU) ------ " << tt2.getwct()
<< "(Wall) Seconds.";
}
}
openfpm_finalize();
}
\ No newline at end of file
include ../../../example.mk
CC=mpic++
LDIR =
OBJ = main.o
%.o: %.cpp
$(CC) -O3 -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
Stokes3dBall: $(OBJ)
$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
all: Stokes3dBall
run: all
mpirun -np 4 /Stokes3dBall
.PHONY: clean all run
clean:
rm -f *.o *~ core Stokes3dBall
//
// Created by Abhinav Singh on 15.06.20.
//
#include "config.h"
#include <iostream>
#include "DCPSE/DCPSE_OP/DCPSE_op.hpp"
#include "DCPSE/DCPSE_OP/DCPSE_Solver.hpp"
#include "Operators/Vector/vector_dist_operators.hpp"
#include "Vector/vector_dist_subset.hpp"
#include "util/SphericalHarmonics.hpp"
int main(int argc, char* argv[])
{
{ openfpm_init(&argc,&argv);
timer tt2;
tt2.start();
size_t grd_sz = 18;
double V_err_eps = 5e-4;
double nu=1.0;
const size_t sz[3] = {grd_sz,grd_sz,grd_sz};
Box<3, double> box({-1.0, -1.0,-1.0}, {1.0,1.0,1.0});
size_t bc[3] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
double spacing = 2.0 / (sz[0] - 1);
double rCut = 3.9*spacing;
double R=1.0;
Ghost<3, double> ghost(rCut);
// P V v_B RHS V_t P_anal RHS2 Polar cord
vector_dist_ws<3, double, aggregate<double,VectorS<3, double>,VectorS<3, double>,double,VectorS<3, double>,double,double,VectorS<3, double>,VectorS<3, double>,VectorS<3, double>>> Particles(0, box, bc, ghost);
Particles.setPropNames({"Pressure","Velocity","Velocity_Boundary","Divergence","V_T","P_analytical","TEMP" ,"RHS","PolarCoord","V_analytical"});
auto &v_cl = create_vcluster();
auto it = Particles.getGridIterator(sz);
while (it.isNext()) {
auto key = it.get();
double x = -1.0+key.get(0) * it.getSpacing(0);
double y = -1.0+key.get(1) * it.getSpacing(1);
double z = -1.0+key.get(2) * it.getSpacing(2);
double r=sqrt(x*x+y*y+z*z);
if (r<R-spacing/2.0) {
Particles.add();
Particles.getLastPos()[0] = x;
Particles.getLastPos()[1] = y;
Particles.getLastPos()[2] = z;
Particles.getLastProp<8>()[0] = r;
if (r==0){
Particles.getLastProp<8>()[1] = 0.0;
}
else{
Particles.getLastProp<8>()[1] = std::atan2(sqrt(x*x+y*y),z);
}
Particles.getLastProp<8>()[2] = std::atan2(y,x);
}
++it;
}
int n_sp=int(grd_sz)*int(grd_sz)*3;
double Golden_angle=M_PI * (3.0 - sqrt(5.0));
for(int i=1;i<=n_sp;i++)
{
double y = 1.0 - (i /double(n_sp - 1.0)) * 2.0;
double radius = sqrt(1 - y * y);
double Golden_theta = Golden_angle * i;
double x = cos(Golden_theta) * radius;
double z = sin(Golden_theta) * radius;
if (acos(z)==0 || acos(z)==M_PI){
std::cout<<"Theta 0/Pi "<<std::endl;
continue;
}
Particles.add();
Particles.getLastPos()[0] = x;
Particles.getLastPos()[1] = y;
Particles.getLastPos()[2] = z;
Particles.getLastProp<8>()[0] = 1.0 ;
Particles.getLastProp<8>()[1] = std::atan2(sqrt(x*x+y*y),z);
Particles.getLastProp<8>()[2] = std::atan2(y,x);
}