Merge branch 'FD_solver' of https://git.mpi-cbg.de/openfpm/openfpm_pdata into FD_solver

.gitignore

@@ -127,3 +127,4 @@ projectId.sh
 
 **/.gitignore
 .gitignore
+/doxygen/

example/Numerics/2D_ActiveFluid/Makefile

+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

example/Numerics/Stoke_flow/2_2D_LidDrivenCavity_PC/Makefile

+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

example/Numerics/Stoke_flow/2_2D_LidDrivenCavity_PC/mainDCPSE.cpp

+//
+// 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"
+
+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;
+        double sampling_factor = 3.1;
+        double sampling_factor2 = 1.9;
+        
+        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

example/Numerics/Stoke_flow/2_2D_LidDrivenCavity_PC/mainFD.cpp

+#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

example/Numerics/Stoke_flow/3_3D_StokesFlowBall/Makefile

+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

example/Numerics/Stoke_flow/3_3D_StokesFlowBall/main.cpp

+//
+// 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);
+        }
+        Particles.map();
+        Particles.ghost_get<0>();
+
+
+        std::unordered_map<const lm,double,key_hash,key_equal> Vr;
+        std::unordered_map<const lm,double,key_hash,key_equal> V1;
+        std::unordered_map<const lm,double,key_hash,key_equal> V2;
+        //Setting max mode l_max
+        constexpr int K = 2;
+        //Setting amplitudes to 0
+        for(int l=0;l<=K;l++){
+            for(int m=-l;m<=l;m++){
+                Vr[std::make_tuple(l,m)]=0.0;
+                V1[std::make_tuple(l,m)]=0.0;
+                V2[std::make_tuple(l,m)]=0.0;
+            }
+
+
+        }
+        //Setting some amplitude for boundary velocity
+        V1[std::make_tuple(2,0)]=1.0;
+
+        auto it2 = Particles.getDomainIterator();
+        while (it2.isNext()) {
+            auto p = it2.get();
+            Point<3, double> xp = Particles.getPos(p);
+            Point<3, double> xP = Particles.getProp<8>(p);
+            Particles.getProp<0>(p) =0;
+            if (xP[0]==1.0) {
+                Particles.getProp<0>(p) =  0;
+                std::vector<double> SVel;
+                SVel=openfpm::math::sumY<K>(xP[0],xP[1],xP[2],Vr,V1,V2);
+                double SP=openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Vr);
+                Particles.getProp<2>(p)[0] = SVel[0];
+                Particles.getProp<2>(p)[1] = SVel[1];
+                Particles.getProp<2>(p)[2] = SVel[2];
+                Particles.getProp<9>(p)[0] = SVel[0];
+                Particles.getProp<9>(p)[1] = SVel[1];
+                Particles.getProp<9>(p)[2] = SVel[2];
+                Particles.getProp<5>(p) = SP;
+                Particles.setSubset(p,1);
+
+            }
+            else {
+                Particles.setSubset(p,0);
+                Particles.getProp<0>(p) =  0;
+                Particles.getProp<1>(p)[0] =  0;
+                Particles.getProp<1>(p)[1] =  0;
+                Particles.getProp<1>(p)[2] =  0;
+            }
+            ++it2;
+        }
+
+        vector_dist_subset<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_bulk(Particles,0);
+        vector_dist_subset<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_surface(Particles,1);
+
+        auto & bulk = Particles_bulk.getIds();
+        auto & Surface = Particles_surface.getIds();
+
+        for (int j = 0; j < bulk.size(); j++) {
+            auto p = bulk.get<0>(j);
+            Point<3, double> xp = Particles.getPos(p);
+            Point<3, double> xP = Particles.getProp<8>(p);
+
+            std::unordered_map<const lm,double,key_hash,key_equal> Ur;
+            std::unordered_map<const lm,double,key_hash,key_equal> U2;
+            std::unordered_map<const lm,double,key_hash,key_equal> U1;
+            std::unordered_map<const lm,double,key_hash,key_equal> Plm;
+
+            for (int l = 0; l <= K; l++) {
+                for (int m = -l; m <= l; m++) {
+                    auto Er= Vr.find(std::make_tuple(l,m));
+                    auto E1= V1.find(std::make_tuple(l,m));
+                    auto E2= V2.find(std::make_tuple(l,m));
+                    std::vector<double> Sol=openfpm::math::sph_anasol_u(nu,l,m,Er->second,E1->second,E2->second,xP[0]);
+                    Ur[std::make_tuple(l,m)]=Sol[0];
+                    U1[std::make_tuple(l,m)]=Sol[1];
+                    U2[std::make_tuple(l,m)]=Sol[2];
+                    Plm[std::make_tuple(l,m)]=Sol[3];
+                }
+
+            }
+
+            if(fabs(xP[0])>=1e-5 && xP[1]>1e-5 && (M_PI-xP[1])>=1e-5)
+            {
+                std::vector<double> SVel = openfpm::math::sumY<K>(xP[0], xP[1], xP[2], Ur, U1, U2);
+                Particles.getProp<9>(p)[0] = SVel[0];
+                Particles.getProp<9>(p)[1] = SVel[1];
+                Particles.getProp<9>(p)[2] = SVel[2];
+                Particles.getProp<5>(p) = openfpm::math::sumY_Scalar<K>(xP[0], xP[1], xP[2], Plm);
+            }
+        }
+
+        auto P = getV<0>(Particles);
+        auto V = getV<1>(Particles);
+        auto V_B = getV<2>(Particles);
+        V.setVarId(0);
+        auto DIV = getV<3>(Particles);
+        auto V_t = getV<4>(Particles);
+        auto P_anal = getV<5>(Particles);
+        auto temp=getV<6>(Particles);
+        auto RHS = getV<7>(Particles);
+        auto P_bulk = getV<0>(Particles_bulk);
+        auto RHS_bulk = getV<7>(Particles_bulk);
+        auto V_anal = getV<9>(Particles);
+
+        V_t=V;
+        P=0;
+        P_bulk=0;
+        eq_id vx,vy,vz;
+
+        vx.setId(0);
+        vy.setId(1);
+        vz.setId(2);
+
+        double sampling=3.1;
+        double sampling2=1.9;
+        double rCut2=3.9*spacing;
+
+        Derivative_x Dx(Particles, 2, rCut,sampling, support_options::RADIUS),B_Dx(Particles_bulk, 2, rCut,sampling, support_options::RADIUS);
+        Derivative_y Dy(Particles, 2, rCut,sampling, support_options::RADIUS),B_Dy(Particles_bulk, 2, rCut,sampling, support_options::RADIUS);
+        Derivative_z Dz(Particles, 2, rCut,sampling, support_options::RADIUS),B_Dz(Particles_bulk, 2, rCut,sampling, support_options::RADIUS);
+        Derivative_xx Dxx(Particles, 2, rCut2,sampling2,support_options::RADIUS);
+        Derivative_yy Dyy(Particles, 2, rCut2,sampling2,support_options::RADIUS);
+        Derivative_zz Dzz(Particles, 2, rCut2,sampling2,support_options::RADIUS);
+
+        //std::cout << "DCPSE KERNELS DONE" << std::endl;
+        petsc_solver<double> solverPetsc;
+        solverPetsc.setPreconditioner(PCNONE);
+        timer tt;
+        double sum=0,sum1=0;
+        V_t=V;
+        //double V_err_eps = 1e-5;
+
+        double V_err = 1, V_err_old;
+        int n = 0;
+        int nmax = 30;
+        int ctr = 0, errctr, Vreset = 0;
+        V_err = 1;
+        n = 0;
+        tt.start();
+        while (V_err >= V_err_eps && n <= nmax) {
+            //Particles.write_frame("StokesSphere",n);
+            Particles.ghost_get<0>(SKIP_LABELLING);
+            RHS_bulk[0] = B_Dx(P);
+            RHS_bulk[1] = B_Dy(P);
+            RHS_bulk[2] = B_Dz(P);
+            DCPSE_scheme<equations3d3, decltype(Particles)> Solver(Particles);
+            auto Stokes1 = nu * (Dxx(V[0])+Dyy(V[0])+Dzz(V[0]));
+            auto Stokes2 = nu * (Dxx(V[1])+Dyy(V[1])+Dzz(V[1]));
+            auto Stokes3 = nu * (Dxx(V[2])+Dyy(V[2])+Dzz(V[2]));
+            Solver.impose(Stokes1, bulk, RHS[0], vx);
+            Solver.impose(Stokes2, bulk, RHS[1], vy);
+            Solver.impose(Stokes3, bulk, RHS[2], vz);
+            Solver.impose(V[0], Surface, V_B[0], vx);
+            Solver.impose(V[1], Surface, V_B[1], vy);
+            Solver.impose(V[2], Surface, V_B[2], vz);
+            Solver.solve_with_solver(solverPetsc, V[0], V[1], V[2]);
+            Particles.ghost_get<1>();
+            DIV = -(Dx(V[0])+Dy(V[1])+Dz(V[2]));
+            P_bulk = P + DIV;
+            sum = 0;
+            sum1 = 0;
+            for (int j = 0; j < bulk.size(); j++) {
+                auto p = bulk.get<0>(j);
+                sum += (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) *
+                       (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) +
+                       (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) *
+                       (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) +
+                       (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]) *
+                       (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]);
+                sum1 += Particles.getProp<1>(p)[0] * Particles.getProp<1>(p)[0] +
+                        Particles.getProp<1>(p)[1] * Particles.getProp<1>(p)[1] +
+                        Particles.getProp<1>(p)[2] * Particles.getProp<1>(p)[2];
+            }
+            sum = sqrt(sum);
+            sum1 = sqrt(sum1);
+            v_cl.sum(sum);
+            v_cl.sum(sum1);
+            v_cl.execute();
+            V_t = V;
+            Particles.ghost_get<1>(SKIP_LABELLING);
+            V_err_old = V_err;
+            V_err = sum / sum1;
+            if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-14) {
+                errctr++;
+            } else {
+                errctr = 0;
+            }
+            if (n > 3) {
+                if (errctr > 1) {
+                    std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN Divergence" << std::endl;
+                    Vreset = 1;
+                    break;
+                } else {
+                    Vreset = 0;
+                }
+            }
+            n++;
+
+        }
+        Particles.write("StokesSphere");
+    }
+    openfpm_finalize();
+
+}
\ No newline at end of file

openfpm_data @ c5d76187

-Subproject commit d78de3919144686d1e88a2c46a88ad3fa2a79043
+Subproject commit c5d7618790eb69bf9fc5ca38f8bd1c0afdcc2476

openfpm_io @ b2b9d762

-Subproject commit 1a8e1bb010f3562596263fea2c01e8b3514ef015
+Subproject commit b2b9d76267afb6548a82e4b8eea541879d2555f0

openfpm_numerics @ b0954918

-Subproject commit 525fd4682cbafca1a15276cd616de77ca777ba8d
+Subproject commit b09549183092eb01315d1e6ca9cf6213ee792b1b