// This codes shows how to find the static solution of // Sine-Gordon equation on a pseudo 1d domain. // -\nabla^2 u+sin(u) = 0 // // Finite energy solutions satisfy asymptotic behaviour // u(x->\infty) = 2k\pi (k\in Z) // Kink solutions interpolate between two different // topological vacua (e.g 0 <- u -> 2\pi ) // // (see e.g Topological Solitons by N. Manton and P. Sutcliffe) // // Minimization is done with TAO routines compared with custom nlcg // // Usage : // mpirun --oversubscribe -np 4 FreeFem++-mpi -wg hpddm-petsc-NLCG.edp /**************************************/ /* Load PETSc macros */ /**************************************/ load "PETSc" // PETSc plugin macro partitioner()metis // End Of Macro // metis, scotch, or parmetis macro dimension( )2 // End Of Macro // 2D or 3D include "macro_ddm.idp" // Additional DDM functions include "ddm-NLCG.idp"; // NLCG DDM functions /**************************************/ /* Macros for distributed solver */ /**************************************/ macro def(i)i // End Of Macro // Definition of scalar fields macro init(i)i // End Of Macro // Initialization of scalar fields func Pk = P1; // Finite-element space macro Grad(u)[dx(u), dy(u)] // End Of Macro ;; /***************************************/ /* Options for geometry and solver */ /***************************************/ int Npts = getARGV("-npts" , 400); // Number of points on the perimeter real Lx = getARGV("-lx" , 20.0); // Dimension of the domain real Ly = getARGV("-ly" , 1.0); // Dimension of the domain // int iMAX = getARGV("-iMAX", 200); // Max iter of NLCG real tol = getARGV("-tol" , 1.0e-08); // Dimension of the domain // real mass = getARGV("-mass", 1.0); // Mass Parameter of Sine-Gordon equatiom /***************************************/ /* ############################### */ /***************************************/ meshN ThNo, Th = square(1, 1); // Local mesh real[int] D; int[int][int] restriction; /***************************************/ /* Finite Element space */ /***************************************/ fespace Vh(Th, Pk); // local finite element space /***************************************/ /* Geometry */ /***************************************/ int[int] Labels=[1,2,3,4]; // labels : bottom, right, top, left sides { // Construction of the rectangular domain int Xpts= int(0.5*Npts*Lx/(Lx+Ly)); // pts on the x-axis sides int Ypts= int(0.5*Npts*Ly/(Lx+Ly)); // pts on the y-axis sides Th = square(Xpts,Ypts,[Lx*(x-0.5),Ly*(y-0.5)],label=Labels); fespace Ph(Th, P0); Ph part; if(mpirank == 0) partitionerSeq(part[], Th, mpisize); partitionerPar(part[], Th, mpiCommWorld, size); buildWithPartitioning(Th, part[], 1, restriction, D, Pk, mpiCommWorld); part = part; ThNo = trunc(Th, abs(part - mpirank) < 1e-1); } Mat H(Vh.ndof, restriction, D); /***************************************/ /* Problem definition */ /***************************************/ // Energy of Sine-Gordon func real Energy(real[int] & XX){ Vh u; u[]=XX; // Copy the d.o.f to local d.o.f u real vEtot,vE= intN(ThNo) // Energy of the problem ( 0.5*Grad(u)'*Grad(u) + mass^2*(1.0-cos(u)) ); mpiAllReduce(vE,vEtot,mpiCommWorld,mpiSUM); return vEtot; } // Sine-Gordon equations of motion func real[int] EOM(real[int] & XX){ Vh u; u[]=XX; varf vDJ(uh,vh) =intN(Th) // Equations of motion (Euler-Lagrange) ( Grad(u)'*Grad(vh) // variation of the kinetic term +mass^2*sin(u)*vh // variation of the potential term ); real[int] AuxVec = vDJ( 0,Vh); exchange(H, AuxVec, scaled = true); return AuxVec; } /***************************************/ /* Problem resolution */ /***************************************/ Vh def(u); // Define the local solution int converged=0; u=(x>0?2*pi:0.0); ddmNlcg(real,H,mpiCommWorld,"FR",converged,EOM,u[],iMAX,-1.0e-08); real E = Energy(u[]); plotMPI(Th,u,Pk,def,real, {cmm = "Final value, E="+E,dim=3,value=1,fill=1,wait=1} ) ;; // Analytical solution u=4*atan(exp(mass*x)); E = Energy(u[]); plotMPI(Th,u,Pk,def,real, {cmm = "Analytical solution "+E,dim=3,value=1,fill=1,wait=1} ) ;;