GfsPoisson

From Gerris

Jump to: navigation, search

GfsPoisson solves the Laplace equation ∇² P = 0, by default, for a scalar variable P; however, if a quantity Div is defined, typically in a GfsInit, this becomes the right-hand side Θ of the Poisson equation ∇² P = Θ.

Boundary conditions can be set on an embedded GfsSurface with GfsSurfaceBc, or on a GfsBoundary with a GfsBc; in either case the default is homogeneous Neumann.


Examples

  • Convergence of the Poisson solver
  • 1 0 GfsPoisson GfsBox GfsGEdge {} {
      Time { iend = 1 }
      Refine LEVEL
    
      GModule hypre
      ApproxProjectionParams { tolerance = 1e-30 nitermin = CYCLE nitermax = CYCLE }
    
      Init {} {
        Div = {
          int k = 3, l = 3;
          return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
        }
      }
      OutputTime { istep = 1 } {
        awk '{if ($2 == 1) print CYCLE, $8;}' >> time
      }
      OutputProjectionStats { start = end } {
        awk '{if ($1 == "residual.infty:") print CYCLE, $3, $4;}' >> proj
      }
      OutputErrorNorm { start = end } {
        awk '{print LEVEL, $5, $7, $9}' >> error 
      } { v = P } {
        s = (sin (M_PI*3.*x)*sin (M_PI*3.*y))
        unbiased = 1
      }
      OutputSimulation { start = end } end-SOLVER.gfs
    }
    

  • Convergence with a refined circle
  • 1 0 GfsPoisson GfsBox GfsGEdge {} {
      Time { iend = 1 }
      Refine (x*x + y*y <= 0.25*0.25 ? LEVEL + 2 : LEVEL)
    
      GModule hypre
      ApproxProjectionParams { tolerance = 1e-30 nitermin = CYCLE nitermax = CYCLE }
    
      Init {} {
        Div = {
          int k = 3, l = 3;
          return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
        }
      }
      OutputTime { istep = 1 } {
        awk '{if ($2 == 1) print CYCLE, $8;}' >> time
      }
      OutputProjectionStats { start = end } {
        awk '{if ($1 == "residual.infty:") print CYCLE, $3, $4;}' >> proj
      }
      OutputErrorNorm { start = end } {
        awk '{print LEVEL, $5, $7, $9}' >> error 
      } { v = P } {
        s = (sin (M_PI*3.*x)*sin (M_PI*3.*y))
        unbiased = 1
      }
      OutputSimulation { start = end } end-SOLVER.gfs
    }
    

  • Dirichlet boundary condition
  • 1 0 GfsPoisson GfsBox GfsGEdge {} {
        Time { iend = 1 }
        Refine LEVEL
    
        Solid ({
                double theta = atan2 (y, x);
                double radius = 0.30 + 0.15*cos (6.*theta);
                return radius*radius - x*x - y*y;
        })
        SurfaceBc P Dirichlet {
    	double theta = atan2 (y, x);
    	double r2 = x*x + y*y;
    	return r2*r2*cos (3.*theta);
        }
    
        GModule hypre
        ApproxProjectionParams { 
    	tolerance = 1e-30 
    	erelax = 2 
    	nitermin = CYCLE nitermax = CYCLE 
        }
    
        Init {} {
    	Div = {
    	    double theta = atan2 (y, x);
    	    double r2 = x*x + y*y;
    	    return 7.*r2*cos (3.*theta);
    	}
        }
        OutputTime { istep = 1 } {
    	awk '{if ($2 == 1) print CYCLE, $8;}' >> time
        }
        OutputProjectionStats { start = end } {
    	awk '{if ($1 == "residual.infty:") print CYCLE, $3, $4;}' >> proj
        }
        OutputErrorNorm { start = end } {
    	awk '{print LEVEL, $5, $7, $9}' >> error 
        } { v = P } {
    	s = {
    	    double theta = atan2 (y, x);
    	    double r2 = x*x + y*y;
    	    return r2*r2*cos (3.*theta);
    	}
        }
        OutputSimulation { start = end } end-SOLVER.gfs
    }
    

  • Convergence of the Poisson solver with solid boundaries
  • 1 0 GfsPoisson GfsBox GfsGEdge {} {
      Time { iend = 1 }
      Refine LEVEL
      Solid (ellipse (0, 0, 0.25, 0.25))
    
      GModule hypre
      ApproxProjectionParams { tolerance = 1e-30 erelax = 2 nitermin = CYCLE nitermax = CYCLE }
    
      Init {} {
        Div = {
          int k = 3, l = 3;
          return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
        }
      }
      OutputTime { istep = 1 } {
        awk '{if ($2 == 1) {print "CYCLE " $8}}' >> time
      }
      OutputProjectionStats { start = end } {
        awk '{
          if ($1 == "residual.infty:") print "CYCLE "$3 " " $4;
        }' >> proj
      }
      OutputSimulation { start = end } sim-LEVEL-SOLVER { variables = P }
    }
    GfsBox {}
    

  • Star-shaped solid boundary with refinement
  • 1 0 GfsPoisson GfsBox GfsGEdge {} {
        Time { iend = 1 }
        Refine LEVEL
        RefineSolid (LEVEL + 2)
        Solid ({
                double dr = 0.1;
                double theta = atan2 (y, x);
                double radius = 0.79*(0.45 - dr + dr*cos (6.*theta));
                return x*x + y*y - radius*radius;
              })
        GModule hypre
        ApproxProjectionParams { tolerance = 1e-30 nitermin = CYCLE nitermax = CYCLE }
        Init {} {
            Div = {
                int k = 3, l = 3;
                return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
            }
        }
        OutputTime { istep = 1 } {
          awk '{if ($2 == 1) {print "CYCLE " $8}}' >> time
        }
        OutputProjectionStats { start = end } {
            awk '{
          if ($1 == "residual.infty:") print "CYCLE "$3 " " $4;
        }' >> proj
        }
        OutputSimulation { start = end } sim-LEVEL-SOLVER { variables = P }
    }
    GfsBox {}
    

  • Star-shaped solid boundary
  • 1 0 GfsPoisson GfsBox GfsGEdge {} {
        Time { iend = 1 }
        Refine LEVEL
        Solid ({
                double dr = 0.1;
                double theta = atan2 (y, x);
                double radius = 0.79*(0.45 - dr + dr*cos (6.*theta));
                return x*x + y*y - radius*radius;
              })
        GModule hypre
        ApproxProjectionParams { erelax = 2 tolerance = 1e-30 nitermin = CYCLE nitermax = CYCLE }
        Init {} {
            Div = {
                int k = 3, l = 3;
                return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
            }
        }
        OutputTime { istep = 1 } {
          awk '{if ($2 == 1) {print "CYCLE " $8}}' >> time
        }
        OutputProjectionStats { start = end } {
            awk '{
          if ($1 == "residual.infty:") print "CYCLE "$3 " " $4;
        }' >> proj
        }
        OutputSimulation { start = end } sim-LEVEL-SOLVER { variables = P }
    }
    GfsBox {}
    

  • Thin wall at box boundary
  • 4 3 GfsPoisson GfsBox GfsGEdge {} {
      Time { iend = 1 }
      Refine LEVEL
      GModule hypre
      ApproxProjectionParams { erelax = 2 tolerance = 1e-30 nitermin = CYCLE nitermax = CYCLE }
      Init {} {
        Div = {
          int k = 3, l = 3;
          x = (x - 0.5)/2.;
          y = (y + 0.5)/2.;
          return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
        }
      }
      OutputTime { istep = 1 } {
        awk '{if ($2 == 1) {print "CYCLE " $8}}' >> time
      }
      OutputProjectionStats { start = end } {
        awk '{
          if ($1 == "residual.infty:") print "CYCLE "$3 " " $4;
        }' >> proj
      }
      OutputSimulation { start = end } sim-LEVEL-SOLVER { variables = P }
    }
    GfsBox {}
    GfsBox {}
    GfsBox {}
    GfsBox {}
    1 2 right
    2 3 bottom
    3 4 left
    

  • Poisson solution in a dumbbell-shaped domain
  • 1 0 GfsPoisson GfsBox GfsGEdge {} {
      Refine 3
      ApproxProjectionParams { nitermax = 1000 minlevel = 1 tolerance = 1e-30 }
      Solid dumbell.gts
      Init {} { Div = y }
      OutputProjectionStats { istep = 1 } stdout
    }
    

  • Poisson problem with a pure spherical harmonics solution
  • 6 12 GfsPoisson GfsBox GfsGEdge {} {
      Time { iend = 1 }
      PhysicalParams { L = 2.*M_PI/4. }
      MetricCubed M LEVEL
    
      # Use alternative implementation
      # MetricCubed1 M
      # MapFunction {
      #     x = atan2 (X, Z)*180./M_PI
      #     y = asin (CLAMP(Y,-1.,1.))*180./M_PI
      # }
    
      Refine LEVEL
    
      GModule hypre
      ApproxProjectionParams { tolerance = 1e-30 nitermin = CYCLE nitermax = CYCLE }
    
      Global {
          #include <gsl/gsl_sf.h>
          @link -lgsl -lgslcblas
    
          double fact (double n) {
            if (n <= 1)
          	  return 1.;
        	else
    	  return n*fact(n - 1.);
          }
    
          double legendre (int l, int m, double x) {
    	  if (m < 0) 
    	      return pow(-1.,fabs(m))*fact(l-fabs(m))/fact(l+fabs(m))*
    	      gsl_sf_legendre_Plm (l, fabs(m), x);
    	  else
    	      return gsl_sf_legendre_Plm (l, m, x);
          }
    
          double spherical_harmonic_re (int l, int m, double X, double Y) {
    	  return sqrt((2*l+1)/(4*M_PI)*fact(l-m)/fact(l+m))*
    	  legendre (l, m, sin(Y/180*M_PI))*cos(m*X/180*M_PI);
          }
      }
    
      Init { } { 
          Div = -4*(4+1)*spherical_harmonic_re (4, 2, x, y) 
          Sol = spherical_harmonic_re (4, 2, x, y)
      }
    
      OutputTime { istep = 1 } {
        awk '{if ($2 == 1) print CYCLE, $8;}' >> time
      }
      OutputProjectionStats { start = end } {
        awk '{if ($1 == "residual.infty:") print CYCLE, $3, $4;}' >> proj
      }
      OutputErrorNorm { start = end } {
        awk '{print LEVEL, $5, $7, $9}' >> error 
      } { v = P } {
        s = Sol
        v = E
        unbiased = 1
      }
      OutputSimulation { start = end } end-SOLVER-LEVEL.gfs
    }
    

  • Spherical harmonics with longitude-latitude coordinates
  • 2 2 GfsPoisson GfsBox GfsGEdge {} {
      Time { iend = 1 }
      PhysicalParams { L = M_PI }
      MetricLonLat M 1.
      Refine LEVEL
    
      GModule hypre
      ApproxProjectionParams { tolerance = 1e-30 nitermin = CYCLE nitermax = CYCLE }
    
      Global {
          #include <gsl/gsl_sf.h>
          @link -lgsl -lgslcblas
    
          double fact (double n) {
            if (n <= 1)
          	  return 1.;
        	else
    	  return n*fact(n - 1.);
          }
    
          double legendre (int l, int m, double x) {
    	  if (m < 0) 
    	      return pow(-1.,fabs(m))*fact(l-fabs(m))/fact(l+fabs(m))*
    	      gsl_sf_legendre_Plm (l, fabs(m), x);
    	  else
    	      return gsl_sf_legendre_Plm (l, m, x);
          }
    
          double spherical_harmonic_re (int l, int m, double X, double Y) {
    	  return sqrt((2*l+1)/(4*M_PI)*fact(l-m)/fact(l+m))*
    	  legendre (l, m, sin(Y/180*M_PI))*cos(m*X/180*M_PI);
          }
      }
    
      Init { } { 
          Div = -4*(4+1)*spherical_harmonic_re (4, 2, x, y) 
          Sol = spherical_harmonic_re (4, 2, x, y)
      }
    
      OutputTime { istep = 1 } {
        awk '{if ($2 == 1) print CYCLE, $8;}' >> time
      }
      OutputProjectionStats { start = end } {
        awk '{if ($1 == "residual.infty:") print CYCLE, $3, $4;}' >> proj
      }
      OutputErrorNorm { start = end } {
        awk '{print LEVEL, $5, $7, $9}' >> error 
      } { v = P } {
        s = Sol
        v = E
        unbiased = 1
      }
      OutputSimulation { start = end } end-SOLVER-LEVEL.gfs
    }
    

  • Poisson equation on a sphere with Gaussian forcing
  • 6 12 GfsPoisson GfsBox GfsGEdge {} {
      PhysicalParams { L = 2.*M_PI/4. }
      MetricCubed M LEVEL
      Time { iend = 1 }
      Refine LEVEL
    
      ProjectionParams { tolerance = 1e-12 }
      ApproxProjectionParams { tolerance = 1e-12 }
    
      Init {} { Div = exp(-2.*5.*5.*(1.-cos((y + 90.)/180.*M_PI))) }
    
      OutputLocation { start = end } prof.dat profile
    
      OutputSimulation { start = end } end.gfs
      OutputProjectionStats { istep = 1 } stderr
    
      EventScript { start = end } {
          gnuplot <<EOF
          set term pos enhanced eps color solid 20 lw 3 
          set out 'profile.eps'
    
          set key bottom right
          set xl "Latitude"
          set yl "{/Symbol F}"
          set xr [-90:90]
    
          plot './prof.dat' u 3:5 every 5 ps 2 t "Gerris",\
               'prof.ref'u 1:2 w l t "(Boyd and Zhou, 2009)"
    EOF
      }
      
    }
    

  • Gaussian forcing using longitude-latitude coordinates
  • 2 2 GfsPoisson GfsBox GfsGEdge {} {
      PhysicalParams { L = M_PI }
      MetricLonLat M 1.
      Time { iend = 1 }
      Refine LEVEL
    
      ProjectionParams { tolerance = 1e-12 }
      ApproxProjectionParams { tolerance = 1e-12 }
    
      Init {} { Div = exp(-2.*5.*5.*(1.-cos((y + 90.)/180.*M_PI))) }
    
      OutputLocation { start = end } prof.dat profile
    
      OutputSimulation { start = end } end.gfs
      OutputProjectionStats { istep = 1 } stderr
    
      EventScript { start = end } {
          gnuplot <<EOF
          set term pos enhanced eps color solid 20 lw 3 
          set out 'profile.eps'
    
          set key bottom right
          set xl "Latitude"
          set yl "{/Symbol F}"
          set xr [-90:90]
    
          plot './prof.dat' u 3:(\$5-0.0078011) every 5 ps 2 t "Gerris",\
               'prof.ref'u 1:2 w l t "(Boyd and Zhou, 2009)"
    EOF
      }
      
    }
    

Personal tools
communication