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
  }
  
}

Retrieved from "http://gfs.sourceforge.net/wiki/index.php/GfsPoisson"
Views
Navigation
communication