hext_py.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2010 Werner Scholz
    Copyright (C) 2009 Amit Itagi

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    magpar is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

static char const Id[] = "$Id: hext_py.c 2983 2010-02-22 22:56:23Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "field.h"
#include "util/util.h"

static int       doinit=1;
static int       fieldon=0;
static Vec       VHthis;
static PetscReal hextamp=0.0;

int Hext_py_Init(GridData *gdata)
{
  MagparFunctionLogBegin;

  PetscTruth flg;
  char str[256];
  ierr = PetscOptionsGetString(PETSC_NULL,"-hext_py",str,256,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    fieldon=0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -hext_py not found, using default value: %i (field off)\n",
      fieldon
    );
  }
  else {
    fieldon=1;
  }

  if (!fieldon) {
    PetscPrintf(PETSC_COMM_WORLD,"Hext_py off\n");
    MagparFunctionLogReturn(0);
  }

#ifdef PYTHON
  PetscPrintf(PETSC_COMM_WORLD,"Hext_py on, using Python script: %s\n",str);

  ierr = VecDuplicate(gdata->M,&VHthis);CHKERRQ(ierr);

  /* Python stuff */

  PyObject *pFunc;
  pFunc=PyObject_GetAttrString(gdata->pModule,str);  /* new reference */
  if(!pFunc) {
    SETERRQ1(PETSC_ERR_ARG_CORRUPT,"Cannot find Python function %s!\n",str);
  }
  if (!PyCallable_Check(pFunc)) {
    SETERRQ1(PETSC_ERR_ARG_CORRUPT,"Python object %s is not callable (not a function)!\n",str);
  }

  /* Arguments for the module */
  PyObject *pTuple;
  pTuple = PyTuple_New(ND); /* new reference */

  PyObject *pVal=PETSC_NULL;

  /* make the vector of the external field accessible */
  PetscReal *ta_Hext;
  ierr = VecGetArray(VHthis,&ta_Hext);CHKERRQ(ierr);

  /* loop over all (local) nodes */
  for (int j=0; j<gdata->ln_vert; j++) {
    ierr = ProgressBar(j,gdata->ln_vert,10);
    PetscReal x,y,z;
    x=gdata->vertxyz[ND*gdata->vertl2g[j]+0];
    y=gdata->vertxyz[ND*gdata->vertl2g[j]+1];
    z=gdata->vertxyz[ND*gdata->vertl2g[j]+2];

    PyTuple_SetItem(pTuple,0,PyFloat_FromDouble(x)); /* The tuple steals the reference */
    PyTuple_SetItem(pTuple,1,PyFloat_FromDouble(y));
    PyTuple_SetItem(pTuple,2,PyFloat_FromDouble(z));

    /* Call the Python function */
    pVal = PyObject_CallObject(pFunc,pTuple);  /* new reference */

    /* set external field */
    for (int i=0; i<ND; i++) {
      ta_Hext[ND*j+i]=PyFloat_AsDouble(PyTuple_GetItem(pVal,i));
    }
  }
  ierr = ProgressBar(1,1,10);
  ierr = VecRestoreArray(VHthis,&ta_Hext);CHKERRQ(ierr);

  /* decrement references */
  Py_DECREF(pFunc);
  Py_DECREF(pTuple);
  Py_DECREF(pVal);

  int k;
  ierr = VecMax(VHthis,&k,&hextamp);CHKERRQ(ierr);
#else
  SETERRQ(PETSC_ERR_ARG_CORRUPT,"magpar not linked with Python interpreter. Cannot execute Python scripts!\n");
#endif

  MagparFunctionLogReturn(0);
}


int hext_f(GridData *gdata, PetscReal time_ns)
{
#ifdef PYTHON
  MagparFunctionInfoBegin;

  char str[256]="hext_f";
  PyObject *pFunc;
  pFunc=PyObject_GetAttrString(gdata->pModule,str);  /* new reference */
  if(!pFunc) {
    SETERRQ1(PETSC_ERR_ARG_CORRUPT,"Cannot find Python function %s!\n",str);
  }
  if (!PyCallable_Check(pFunc)) {
    SETERRQ1(PETSC_ERR_ARG_CORRUPT,"Python object %s is not callable (not a function)!\n",str);
  }

  /* Arguments for the module */
  PyObject *pTuple;
  pTuple = PyTuple_New(ND); /* new reference */

  PyObject *pVal=PETSC_NULL;

  PyTuple_SetItem(pTuple,0,PyFloat_FromDouble(time_ns)); /* The tuple steals the reference */

  /* Call the Python function */
  pVal = PyObject_CallObject(pFunc,pTuple);  /* new reference */

  PetscReal f;
  f=PyFloat_AsDouble(PyTuple_GetItem(pVal,0));

  /* decrement references */
  Py_DECREF(pFunc);
  Py_DECREF(pTuple);
  Py_DECREF(pVal);

  MagparFunctionInfoReturn(f);
#else
  SETERRQ(PETSC_ERR_ARG_CORRUPT,"magpar not linked with Python interpreter. Cannot execute Python scripts!\n");
#endif
}

int Hext_py(GridData *gdata,Vec VHextsum,PetscReal *hext)
{
  MagparFunctionInfoBegin;

  /* to be implemented by the user */

  if (doinit) {
    ierr = Hext_py_Init(gdata);CHKERRQ(ierr);
    doinit=0;
  }
  if (!fieldon) {
    MagparFunctionInfoReturn(0);
  }

  ierr = VecAXPY(VHextsum,1.0,VHthis);CHKERRQ(ierr);

  /* get time in nanoseconds */
  PetscReal nstime;
  nstime=gdata->time*gdata->tscale*1e9;

  /* get scaling factor */
  PetscReal hext2;
/* DEBUG */
#if 0
  hext2=hext_f(gdata,nstime);
#else
  hext2=0;
#endif

  ierr = VecAXPY(VHextsum,hext2,VHthis);CHKERRQ(ierr);

  *hext=hextamp;

  MagparFunctionProfReturn(0);
}