DifDop Home

DifDop


Angular Momentum Transport Module

The angular momentum transport module that incorporates in particular the effect of $ E\times B $ shear driven residual stress, is written as an example program to the pdsol library. it can be found in the examples/angular directory.

The system of equations that are being solved can be written as:

\begin{subequations} \begin{align} \frac{\partial n}{\partial t}+\frac{1}{r}\frac{\partial}{\partial r}\left(r\Gamma_{n}\right) & =S_{n}\label{eq:mod1}\\ \frac{\partial P}{\partial t}+\frac{1}{r}\frac{\partial}{\partial r}\left(rQ\right) & =H\\ \frac{\partial L_{\phi}}{\partial t}+\frac{1}{r}\frac{\partial}{\partial r}\left(r\Pi_{r,\phi}\right) & =\tau_{\phi}\;\mbox{.} \end{align} \end{subequations}

Here n is the plasma density, P is the pressure and $L_{\phi}$ is the angular momentum density. The fluxes are defined as:

\begin{subequations} \begin{align} \Gamma_{n} & =-D_{0}\frac{\partial n}{\partial r}-D_{1}\varepsilon\left(\frac{\partial n}{\partial r}+V_{rn}n\right)\\ \Pi_{r,\phi} & =-\nu_{0}\frac{\partial L_{\phi}}{\partial r}-\nu_{1}\varepsilon\left[\frac{\partial L_{\phi}}{\partial r}+V_{rL}L_{\phi}\right]+S\\ \mbox{where}\quad S & =-\varepsilon\alpha\left(r\right)n\left(r\right)\left(1-\frac{\sigma}{P_{0}}\frac{\partial P}{\partial r}\right)\frac{\partial v_{Ey}}{\partial r}\\ Q & =-\chi_{0}\frac{\partial P}{\partial r}-\chi_{1}\varepsilon\frac{\partial P}{\partial r}\;\mbox{.} \end{align} \end{subequations}

We use a reduced force balance equation, a simple formula for turbulence reduction and relate the flux surface average of the flow to angular momentum:

\begin{subequations} \begin{align} \frac{\partial v_{Ey}}{\partial r} & =\left(\frac{\rho_{*}}{n_{0}P_{0}}\right)\frac{\partial n}{\partial r}\frac{\partial P}{\partial r}+\cdots\\ \varepsilon & =\frac{\varepsilon_{0}}{\left[1+\beta\left(\partial v_{Ey}/\partial r\right)^{2}\right]}\\ v_{\phi}\left(r,t\right) & =\left[\frac{L_{\phi}\left(r,t\right)}{n\left(r,t\right)}\right]\lambda_{2}\left(r\right)\;\mbox{.} \end{align} \end{subequations}

We also impose the boundary conditions:

\begin{align*} L_{\phi}(a) & =v_{\phi}(a)=P(a)=0\\ n(a) & =n_{a}\;\mbox{.} \end{align*}

The description of the problem requires 3 evolving fields. 3 auxillary fields are used to store the fluxes. One extra output field is used to output the velocity field (instead of angular momentum). The model can be found in the directory examples/angular. In this directory three files are found:

  • The init.xml file where the problem and integration parameters are defined.
    <?xml version="1.0"?>
    <solver nfields="3" naux="3" nexout="1" dt="1e-8" dtout="0.1" tmax="1.0" name="angular">
    
      <gsolver 
        method="rkf45"
        step="adaptive"
        eps_abs="1e-6"
        eps_rel="1e-8">
      </gsolver>
    
      <field
        index="1" 
        name="P" 
        init="P_init.h5"> 
      </field>
    
      <field
        index="2" 
        name="n" 
        init="n_init.h5">
      </field>
    
      <field
        index="3" 
        name="L" 
        init="L_init.h5">
      </field>
    
      <exout
        index="1"
        name="vphi">
      </exout>
    
      <mesh 
        x0="0.0" 
        dx="0.005"
        nx="201">
      </mesh>
    
      <params 
        D0="1.0"
        D1="4.0"
        chi0="1.0"
        chi1="4.0"
        nu0="1.0"
        nu1="4.0"
        alpha="2.0"
        rhostar="0.1"
        sigma="-0.1"
        gama_a="3.3"
        q_a="3.3"
        eta="10.0"
        eps0="1.0"
        lambda="2.0"
        Vrn="0.0"
        VrL="0.0"
        R="3.0">
      </params>
    
    </solver>
    
  • The "init_fields.m" file, which is an octave file that defines the initial condition.
    dx=0.005;
    xx=0.0:dx:1.0;
    P=1.5*(tanh(100*(-xx+0.9))+1.0);
    save('-hdf5',"P_init.h5","P");
    n=0.5*(tanh(100*(-xx+0.9))+1.0)+0.05;
    save('-hdf5',"n_init.h5","n");
    L=0.0*(tanh(100*(-xx+0.9))+1.0);
    save('-hdf5',"L_init.h5","L");
    
  • The "angular.c" file where the equation itself is defined.
    #include "mesh.h"
    #include "solver.h"
    #include <stdio.h>
    #include <mxml.h>
    #include <math.h>
    #include <string.h>
    #include "pdutils.h"
    #include <gsl/gsl_errno.h>
    #include <gsl/gsl_odeiv.h>
    
    struct pars {
      double D0,D1,chi0,chi1,nu0,nu1,alpha,rhostar,sigma,gama_a,q_a,eta,eps0,lambda,Vrn,VrL,R;
    };
    
    int fcal_angular(double t, const double y[], double dydt[], void *params){
      pdsol_rsolver* sv = (pdsol_rsolver *)params;
      pdsol_mesh *rm = sv->rm;
      const double *P, *n,*L;
      double *dPdt,*dndt,*dLdt;
      double dxp,dxm,dxc;
      double *Sn,*SP,*SL;
      double Pl,nl,Ll,dPldx,dnldx,dLldx,dvEydx,r,n_a,g_a,vareps;
      double alpar,lambda2;
      struct pars *ps=(struct pars *)sv->user;
      int l;
    
      /*setting up the fields*/
      P=&y[0];
      dPdt=&dydt[0];
    
      n=&y[rm->nx];
      dndt=&dydt[rm->nx];
    
      L=&y[2*rm->nx];
      dLdt=&dydt[2*rm->nx];
    
      SP=&sv->aux[0];
      Sn=&sv->aux[rm->nx];
      SL=&sv->aux[2*rm->nx];
    
      n_a=(n[rm->nx-1]+n[rm->nx-2])/2; // note: maybe we should just use n[rm->nx-1] ?
      g_a=(n[rm->nx-2]-n[rm->nx-1])/(rm->xx[rm->nx-1]-rm->xx[rm->nx-2]); // -dn/dx
      /*computing the fluxes */
      
      for(l=0;l<rm->nx-1;l++){
    
        r=(rm->xx[l+1]+rm->xx[l])/2;
        alpar=ps->alpha*(1+r*r/ps->R/ps->R);
        lambda2=(1-r*r/ps->R/ps->R);
        /*Spatial stepsizes:*/
        dxp=rm->xx[l+1]-rm->xx[l];
    
        sv->exout[0]->dat[l]=L[l]*lambda2/n[l];
    
        Pl=(P[l+1]+P[l])/2;
        nl=(n[l+1]+n[l])/2;
        Ll=(L[l+1]+L[l])/2;
    
        dPldx=(P[l+1]-P[l])/dxp;
        dnldx=(n[l+1]-n[l])/dxp;
        dLldx=(L[l+1]-L[l])/dxp;
    
        /*ExB shear:*/
        dvEydx=ps->rhostar*dPldx*dnldx/(Pl*nl);
    
        vareps=ps->eps0/(1+ps->lambda*pow(dvEydx,2));
    
        /* Fluxes:*/
        /* Note: S[l] is the flux at l+1/2 */
        SP[l]=-ps->chi0*dPldx-ps->chi1*vareps*dPldx-2*ps->q_a*r*(1-r*r/2);
        Sn[l]=-ps->D0*dnldx-ps->D1*vareps*dnldx+ps->D1*vareps*ps->Vrn*nl
          -ps->gama_a*exp(-ps->eta*(n_a*(1.0-r)+g_a*(1.0-r)*(1.0-r)/2.0));
        SL[l]=-ps->nu0*dLldx-ps->nu1*vareps*dLldx+ps->nu1*vareps*ps->VrL*Ll
          -alpar*vareps*nl*(1-ps->sigma*dPldx/Pl)*dvEydx;
      }
      SP[rm->nx-1]=SP[rm->nx-2];
      Sn[rm->nx-1]=Sn[rm->nx-2];
      SL[rm->nx-1]=SL[rm->nx-2];
      //SP[rm->nx-1]=0;
      //Sn[rm->nx-1]=0;
    
      /* the actual time step */
    
      for(l=1;l<rm->nx-1;l++){
        dxp=rm->xx[l+1]-rm->xx[l];
        dxm=rm->xx[l]-rm->xx[l-1];
        dxc=(dxp+dxm)/2;
        
        /*Since S[l] is at l+1/2 and S[l-1] at l-1/2, we write the equations as:*/
        dPdt[l]=-(SP[l]*(rm->xx[l+1]+rm->xx[l])/2-SP[l-1]*(rm->xx[l]+rm->xx[l-1])/2)/dxc/rm->xx[l];
        dndt[l]=-(Sn[l]*(rm->xx[l+1]+rm->xx[l])/2-Sn[l-1]*(rm->xx[l]+rm->xx[l-1])/2)/dxc/rm->xx[l];
        dLdt[l]=-(SL[l]*(rm->xx[l+1]+rm->xx[l])/2-SL[l-1]*(rm->xx[l]+rm->xx[l-1])/2)/dxc/rm->xx[l];
        //    printf("dndt[%i]=%f\n",l,dndt[l]);
      }
    
      /* the boundary conditions */
    
      dPdt[0]=dPdt[1];
      dPdt[rm->nx-1]=0;
    
      dndt[0]=dndt[1];
      dndt[rm->nx-1]=0;
    
      dLdt[0]=dLdt[1];
      dLdt[rm->nx-1]=0;
    
      if(L[2]!=L[2]){
        printf("Nan in L!");
        exit(-1);
      }
      return GSL_SUCCESS;
    }
    
    struct pars *get_params(pdsol_rsolver *sv){
      struct pars *ps;
      ps=malloc(sizeof(struct pars));
      ps->D0=pdget_par(sv,"D0");
      ps->D1=pdget_par(sv,"D1");
      ps->chi0=pdget_par(sv,"chi0");
      ps->chi1=pdget_par(sv,"chi1");
      ps->nu0=pdget_par(sv,"nu0");
      ps->nu1=pdget_par(sv,"nu1");
      ps->alpha=pdget_par(sv,"alpha");
      ps->rhostar=pdget_par(sv,"rhostar");
      ps->sigma=pdget_par(sv,"sigma");
      ps->gama_a=pdget_par(sv,"gama_a");
      ps->q_a=pdget_par(sv,"q_a");
      ps->eta=pdget_par(sv,"eta");
      ps->eps0=pdget_par(sv,"eps0");
      ps->lambda=pdget_par(sv,"lambda");
      ps->Vrn=pdget_par(sv,"Vrn");
      ps->VrL=pdget_par(sv,"VrL");
      ps->R=pdget_par(sv,"R");
      return ps;
    }
    
    int main (int argc, char *argv[]){
      double t=0;
      char *confname;
      pdsol_rsolver *sv;
      struct pars *ps;
      int l;
      if (argc==1){
        confname=strdup("init.xml");
      }
      else if (argc==2){
        if (!strcmp(argv[1],"-h")){
          printf("usage: ./angular <config.xml>\n");
          exit(-1);
        }
        confname=strdup(argv[1]);
      }
      else{
        printf("usage: ./angular <config.xml>\n");
        exit(-1);
      }
      sv=pdsol_init_rsolver(confname);
      //  sv->exout[0]->dat=&sv->aux[sv->rm->nx];
      sv->exout[0]->dat=malloc(sizeof(double)*sv->rm->nx);
      memset(sv->exout[0]->dat,0,sizeof(double)*sv->rm->nx);
      sv->gsv->sys->function=fcal_angular;
      ps=get_params(sv);
      sv->user=ps;
      pdsol_gslrun(sv);
      return 0;
    }
    

Fusion Links

Here is a list of links that you might find interesting if you want to learn more about fusion plasmas


This page is currently under construction.