diff options
author | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2013-03-28 01:47:14 +0000 |
---|---|---|
committer | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2013-03-28 01:47:14 +0000 |
commit | d14686b8cb54a26334a7a66a05f5d7910203e966 (patch) | |
tree | 0b2bfb0696b64d52b1b4fd1c9c435d6d402ce57a /src/GRHydro_Con2PrimM_ptee.c | |
parent | 7d408b4489def3a4024b9e0adbba9d9832876950 (diff) |
Finish the implementation of con2prim using the entropy equation.
* Tested successfully pointwisely. It still needs to be tested on evolution.
From: Bruno Coutinho Mundim <bcmsma@astro.rit.edu>
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@500 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
Diffstat (limited to 'src/GRHydro_Con2PrimM_ptee.c')
-rw-r--r-- | src/GRHydro_Con2PrimM_ptee.c | 592 |
1 files changed, 592 insertions, 0 deletions
diff --git a/src/GRHydro_Con2PrimM_ptee.c b/src/GRHydro_Con2PrimM_ptee.c new file mode 100644 index 0000000..de58d4c --- /dev/null +++ b/src/GRHydro_Con2PrimM_ptee.c @@ -0,0 +1,592 @@ +/*********************************************************************************** + Copyright 2006 Scott C. Noble, Charles F. Gammie, Jonathan C. McKinney, + and Luca Del Zanna. + + PVS_GRMHD + + This file was derived from PVS_GRMHD. The authors of PVS_GRMHD include + Scott C. Noble, Charles F. Gammie, Jonathan C. McKinney, and Luca Del Zanna. + PVS_GRMHD is available under the GPL from: + http://rainman.astro.uiuc.edu/codelib/ + + You are morally obligated to cite the following paper in his/her + scientific literature that results from use of this file: + + [1] Noble, S. C., Gammie, C. F., McKinney, J. C., \& Del Zanna, L. \ 2006, + Astrophysical Journal, 641, 626. + + PVS_GRMHD 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. + + PVS_GRMHD 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 PVS_GRMHD; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + + If the user has any questions, please direct them to Scott C. Noble at + scn@astro.rit.edu . + +***********************************************************************************/ + + + +#include <stdio.h> +#include <stdlib.h> +#include <stdarg.h> +#include <string.h> +#include <math.h> +#include <float.h> +#include <complex.h> + +#include "cctk.h" +#include "cctk_Parameters.h" + +/* Set this to be 1 if you want debug output */ +#define DEBUG_CON2PRIMM (0) + + +/* Adiabatic index used for the state equation */ + +#define MAX_NEWT_ITER (30) /* Max. # of Newton-Raphson iterations for find_root_2D(); */ +#define NEWT_TOL (1.0e-10) /* Min. of tolerance allowed for Newton-Raphson iterations */ +#define MIN_NEWT_TOL (1.0e-10) /* Max. of tolerance allowed for Newton-Raphson iterations */ +#define EXTRA_NEWT_ITER (2) + +#define NEWT_TOL2 (1.0e-15) /* TOL of new 1D^*_{v^2} gnr2 method */ +#define MIN_NEWT_TOL2 (1.0e-10) /* TOL of new 1D^*_{v^2} gnr2 method */ + +#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed + to always be smaller than this. This + is used to detect solver failures */ + +#define FAIL_VAL (1.e30) /* Generic value to which we set variables when a problem arises */ + +/************************************************** + The following functions assume a Gamma-law EOS: +***************************************************/ + +/* Local Globals */ +struct LocGlob { + CCTK_REAL Bsq, QdotBsq, Qtsq, Qdotn, D, half_Bsq, Sc, g_o_gm1, + W_for_gnr2, rho_for_gnr2, W_for_gnr2_old, rho_for_gnr2_old, drho_dW ; +} ; + +// Declarations: + + +void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( + CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *prec, + CCTK_REAL *gamma_eos, + CCTK_REAL *dens_in, + CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, + CCTK_REAL *tau_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, + CCTK_REAL *entropycons_in, + CCTK_REAL *y_e_in, CCTK_REAL* temp_in, + CCTK_REAL *rho, + CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, + CCTK_REAL *epsilon, CCTK_REAL *pressure, + CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, + CCTK_REAL *bsq, + CCTK_REAL *w_lorentz, + CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, + CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, + CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, + CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, + CCTK_REAL *det, + CCTK_INT *epsnegative, + CCTK_REAL *retval) ; + +static CCTK_INT general_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, + struct LocGlob *lgp, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][1], CCTK_REAL *, + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ); + +static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, + CCTK_REAL gammaeos, struct LocGlob *lgp); + +/**********************************************************************/ +/********************************************************************************** + + Con2PrimM_pt(): + ----------------------------- + + -- Attempts an inversion from GRMHD conserved variables to primitive variables assuming a guess. + + -- Uses the 2D method of Noble et al. (2006): + -- Solves for two independent variables (W,v^2) via a 2D + Newton-Raphson method + -- Can be used (in principle) with a general equation of state. + + -- Minimizes two residual functions using a homemade Newton-Raphson routine. + -- It is homemade so that it can catch exceptions and handle them correctly, plus it is + optimized for this particular problem. + + -- Note that the notation used herein is that of Noble et al. (2006) except for the argument + list. + + +INPUT: (using GRHydro variable defintions) + + s[x,y,z] = scons[0,1,2] = \alpha \sqrt(\gamma) T^0_i + dens, tau = as defined in GRHydro and are assumed to be densitized (i.e. with sqrt(\gamma)) + dens = D = \sqrt(\gamma) W \rho + tau = \alpha^2 \sqrt(\gamma) T^{00} - D + g[x,y,z][x,y,x] = spatial metric corresponding to \gamma + u[x,y,z][x,y,z] = inverse of the spatial metric, g[x,y,z][x,y,x] + det = sqrt(\gamma) + B[x,y,z] = Bvec[0,1,2] + bsq = b^\mu b_\mu + + epsnegative = (integer) + = 0 if rho and epsilon are positive + != 0 otherwise + + + -- (users should set B[x,y,z] = 0 for hydrodynamic runs) + + +OUTPUT: (using GRHydro variable defintions) + rho, eps = as defined in GRHydro, primitive variables + vel[x,y,z] = as defined in GRHydro, primitive variables + + +RETURN VALUE: of retval = (i*100 + j) where + i = 0 -> Newton-Raphson solver either was not called (yet or not used) + or returned successfully; + 1 -> Newton-Raphson solver did not converge to a solution with the + given tolerances; + 2 -> Newton-Raphson procedure encountered a numerical divergence + (occurrence of "nan" or "+/-inf" ; + + j = 0 -> success + 1 -> failure: some sort of failure in Newton-Raphson; + 2 -> failure: unphysical vsq = v^2 value at initial guess; + 3 -> failure: W<0 or W>W_TOO_BIG + 4 -> failure: v^2 > 1 + ( used to be 5 -> failure: rho,uu <= 0 but now sets epsnegative to non-zero ) + +**********************************************************************************/ + +void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( + CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *prec, + CCTK_REAL *gamma_eos, + CCTK_REAL *dens_in, + CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, + CCTK_REAL *tau_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, + CCTK_REAL *entropycons_in, + CCTK_REAL *y_e_in, CCTK_REAL* temp_in, + CCTK_REAL *rho, + CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, + CCTK_REAL *epsilon, CCTK_REAL *pressure, + CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, + CCTK_REAL *bsq, + CCTK_REAL *w_lorentz, + CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, + CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, + CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, + CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, + CCTK_REAL *det, + CCTK_INT *epsnegative, + CCTK_REAL *retval) + +{ + CCTK_REAL x_1d[1]; + CCTK_REAL sx, sy, sz; + CCTK_REAL usx, usy, usz; + CCTK_REAL tau, dens, gammaeos; + CCTK_REAL QdotB; + CCTK_REAL rho0,u,p,w,gammasq,gamma,W,vsq; + CCTK_REAL g_o_WBsq, QdB_o_W; + CCTK_REAL detg = (*det); + CCTK_REAL sqrt_detg = sqrt(detg); + CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; + CCTK_REAL rho_gm1; + CCTK_REAL utsq; + + DECLARE_CCTK_PARAMETERS; + + struct LocGlob lg; + + gammaeos = *gamma_eos; + + /* Assume ok initially: */ + // BCM: But let the driver function take care of its initialization + //*retval = 0.; + *epsnegative = 0; + +#if(DEBUG_CON2PRIMM) + fprintf(stdout," *dens = %26.16e \n", *dens_in ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *tau = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); + fprintf(stdout," *entropycons = %26.16e \n", *entropycons_in ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + fprintf(stdout," *vely = %26.16e \n", *vely ); + fprintf(stdout," *velz = %26.16e \n", *velz ); + fprintf(stdout," *epsilon = %26.16e \n", *epsilon ); + fprintf(stdout," *pressure = %26.16e \n", *pressure ); + fprintf(stdout," *Bx = %26.16e \n", *Bx ); + fprintf(stdout," *By = %26.16e \n", *By ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); + fprintf(stdout," *retval = %26.16e \n", *retval ); + fflush(stdout); +#endif + + /* First undensitize all conserved variables : */ + sx = ( *sx_in) * inv_sqrt_detg; + sy = ( *sy_in) * inv_sqrt_detg; + sz = ( *sz_in) * inv_sqrt_detg; + tau = ( *tau_in) * inv_sqrt_detg; + dens = (*dens_in) * inv_sqrt_detg; + + usx = (*uxx)*sx + (*uxy)*sy + (*uxz)*sz; + usy = (*uxy)*sx + (*uyy)*sy + (*uyz)*sz; + usz = (*uxz)*sx + (*uyz)*sy + (*uzz)*sz; + + *Bx = (*Bconsx_in) * inv_sqrt_detg; + *By = (*Bconsy_in) * inv_sqrt_detg; + *Bz = (*Bconsz_in) * inv_sqrt_detg; + + lg.Sc = (*entropycons_in) * inv_sqrt_detg; + // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables: + + lg.g_o_gm1 = gammaeos/(gammaeos-1.0); + + lg.Bsq = + (*gxx) * (*Bx) * (*Bx) + + (*gyy) * (*By) * (*By) + + (*gzz) * (*Bz) * (*Bz) + + 2*( + (*gxy) * (*Bx) * (*By) + + (*gxz) * (*Bx) * (*Bz) + + (*gyz) * (*By) * (*Bz) ); + + QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ; + lg.QdotBsq = QdotB*QdotB ; + + lg.Qdotn = -(tau + dens) ; + + lg.Qtsq = (usx * sx + usy * sy + usz * sz) ; + + lg.D = dens; + + lg.half_Bsq = 0.5*lg.Bsq; + + /* calculate W from last timestep and use for guess */ + vsq = + (*gxx) * (*velx) * (*velx) + + (*gyy) * (*vely) * (*vely) + + (*gzz) * (*velz) * (*velz) + + 2*( + (*gxy) * (*velx) * (*vely) + + (*gxz) * (*velx) * (*velz) + + (*gyz) * (*vely) * (*velz) ); + + if( (vsq < 0.) && (fabs(vsq) < 1.0e-13) ) { + vsq = fabs(vsq); + } + if(vsq < 0. || vsq > 1. ) { + *retval = 2.; + return; + } + + gammasq = 1. / (1. - vsq); + gamma = sqrt(gammasq); + + // Always calculate rho from D and gamma so that using D in EOS remains consistent + // i.e. you don't get positive values for dP/d(vsq) . + rho0 = lg.D / gamma ; + rho_gm1 = pow(rho0,(gammaeos-1.)); + p = lg.Sc * rho_gm1 / gamma; + u = p / (gammaeos-1.); + w = rho0 + u + p ; + +// W_last = w*gammasq ; + + // Calculate W and vsq: + x_1d[0] = rho0; +// *retval = 1.0*twod_newton_raphson( x_2d, gammaeos, &lg, func_vsq ) ; + *retval = general_newton_raphson( x_1d, gammaeos, &lg, func_rho ) ; + rho0 = x_1d[0]; + + /* Problem with solver, so return denoting error before doing anything further */ + if( ((*retval) != 0.) || (rho0 == FAIL_VAL) ) { + *retval = *retval*100.+1.; + return; + } + else{ + if( rho0 > W_TOO_BIG) { + *retval = 3.; + return; + } + } + + // Calculate v^2: + utsq = (lg.D-rho0)*(lg.D+rho0)/(rho0*rho0); + gammasq = 1.+utsq; + gamma = sqrt(gammasq); + + if( utsq < 0. ) { + *retval = 4.; + return; + } + + + // Recover the primitive variables from the scalars and conserved variables: + rho0 = lg.D / gamma; + rho_gm1 = pow(rho0,(gammaeos-1.)); + p = lg.Sc * rho_gm1 / gamma; + u = p / (gammaeos-1.); + w = rho0 + u + p ; + W = w * gammasq; + + // User may want to handle this case differently, e.g. do NOT return upon + // a negative rho/u, calculate v^i so that rho/u can be floored by other routine: + if( (rho0 <= 0.) || (u <= 0.) ) { + *epsnegative = 1; + return; + } + + *rho = rho0; + *epsilon = u / rho0; + *w_lorentz = gamma; + *pressure = p ; + + g_o_WBsq = 1./(W+lg.Bsq); + QdB_o_W = QdotB / W; + *bsq = lg.Bsq / gammasq + QdB_o_W*QdB_o_W; + + *velx = g_o_WBsq * ( usx + QdB_o_W*(*Bx) ) ; + *vely = g_o_WBsq * ( usy + QdB_o_W*(*By) ) ; + *velz = g_o_WBsq * ( usz + QdB_o_W*(*Bz) ) ; + + if (*rho <= rho_abs_min*(1.0+GRHydro_atmo_tolerance) ) { + *rho = rho_abs_min; + *velx = 0.0; + *vely = 0.0; + *velz = 0.0; + *w_lorentz = 1.0; + } + + +#if(DEBUG_CON2PRIMM) + fprintf(stdout,"rho = %26.16e \n",*rho ); + fprintf(stdout,"epsilon = %26.16e \n",*epsilon ); + fprintf(stdout,"pressure = %26.16e \n",*pressure ); + fprintf(stdout,"w_lorentz = %26.16e \n",*w_lorentz); + fprintf(stdout,"bsq = %26.16e \n",*bsq ); + fprintf(stdout,"velx = %26.16e \n",*velx ); + fprintf(stdout,"vely = %26.16e \n",*vely ); + fprintf(stdout,"velz = %26.16e \n",*velz ); + fprintf(stdout,"gammaeos = %26.16e \n",gammaeos ); + fflush(stdout); +#endif + + /* done! */ + return; + +} + +/**********************************************************************/ +/************************************************************ + + general_newton_raphson(): + + -- performs Newton-Rapshon method on an arbitrary system. + + -- inspired in part by Num. Rec.'s routine newt(); + + Arguements: + + -- x[] = set of independent variables to solve for; + -- n = number of independent variables and residuals; + -- funcd = name of function that calculates residuals, etc.; + +*****************************************************************/ +static CCTK_INT general_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, + struct LocGlob *lgp, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][1], CCTK_REAL *, + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ) +{ + CCTK_REAL f, df, dx[1], x_old[1], resid[1], + jac[1][1]; + CCTK_REAL errx, x_orig[1]; + CCTK_INT n_iter, i_extra, doing_extra; + CCTK_REAL W,W_old; + + int keep_iterating; + + + // Initialize various parameters and variables: + errx = 1. ; + df = f = 1.; + i_extra = doing_extra = 0; + x_old[0] = x_orig[0] = x[0] ; + + W = W_old = 0.; + + n_iter = 0; + + + /* Start the Newton-Raphson iterations : */ + keep_iterating = 1; + while( keep_iterating ) { + + (*funcd) (x, dx, resid, jac, &f, &df, gammaeos, lgp); /* returns with new dx, f, df */ + + /* Save old values before calculating the new: */ + errx = 0.; + x_old[0] = x[0] ; + /* don't use line search : */ + x[0] += dx[0] ; + + /****************************************/ + /* Calculate the convergence criterion */ + /****************************************/ + + /* For the new criterion, always look at error in "W" : */ + // METHOD specific: + errx = (x[0]==0.) ? fabs(dx[0]) : fabs(dx[0]/x[0]); + + + /****************************************/ + /* Make sure that the new x[] is physical : */ + /****************************************/ + x[0] = fabs(x[0]); + + + /*****************************************************************************/ + /* If we've reached the tolerance level, then just do a few extra iterations */ + /* before stopping */ + /*****************************************************************************/ + + if( (fabs(errx) <= NEWT_TOL) && (doing_extra == 0) && (EXTRA_NEWT_ITER > 0) ) { + doing_extra = 1; + } + + if( doing_extra == 1 ) i_extra++ ; + + if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) || + (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + keep_iterating = 0; + } + + n_iter++; + + } // END of while(keep_iterating) + + + /* Check for bad untrapped divergences : */ + if( (!finite(f)) || (!finite(df)) || (!finite(x[0])) ) { + return(2); + } + + + if( fabs(errx) > MIN_NEWT_TOL){ + return(1); + } + if( (fabs(errx) <= MIN_NEWT_TOL) && (fabs(errx) > NEWT_TOL) ){ + return(0); + } + if( fabs(errx) <= NEWT_TOL ){ + return(0); + } + + return(0); + +} + +/**********************************************************************/ +/********************************************************************************* + func_rho(): + + -- residual/jacobian routine to calculate rho Qtsq equation with + the definition of W + W = ( 1 + GAMMA * K_atm * rho^(GAMMA-1)/(GAMMA-1) ) D^2 / rho + substituted in. + + Arguments: + x = current value of independent var's (on input & output); + dx = Newton-Raphson step (on output); + resid = residuals based on x (on output); + jac = Jacobian matrix based on x (on output); + f = resid.resid/2 (on output) + df = -2*f; (on output) + n = dimension of x[]; + *********************************************************************************/ +static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, + CCTK_REAL gammaeos, struct LocGlob *lgp) +{ + + CCTK_REAL t1,t2; + CCTK_REAL t12; + CCTK_REAL t17; + CCTK_REAL t3, t100,rhosq,t200,rho,W,dWdrho,dvsqdrho,vsq; + + rho = x[0]; + rhosq = rho*rho; + t200 = 1./(lgp->D*lgp->D); + t100 = lgp->g_o_gm1*lgp->Sc*pow(rho,(gammaeos-1.)); + W = lgp->D*( lgp->D + t100 ) / rho ; + dWdrho = lgp->D * ( -lgp->D + t100*(gammaeos-2.) ) / rhosq; + t1 = W*W; + t2 = lgp->Bsq+W; + // t3 = pow(Bsq+W,2.0); + t3 = t2*t2; + vsq = (lgp->D-rho)*(lgp->D+rho)*t200; + dvsqdrho = -2*rho*t200; + resid[0] = t1*(lgp->Qtsq-vsq*t3)+lgp->QdotBsq*(t2+W); + t12 = lgp->Bsq*lgp->Bsq; + t17 = dWdrho*vsq; + jac[0][0] = 2*lgp->QdotBsq*dWdrho + +((lgp->Qtsq-vsq*t12)*2*dWdrho+(-6*t17*lgp->Bsq-dvsqdrho*t12 + +(-2*dvsqdrho*lgp->Bsq-4*t17-dvsqdrho*W)*W)*W)*W; + + dx[0] = -resid[0]/jac[0][0]; + *f = 0.5*resid[0]*resid[0]; + *df = -2. * (*f); + + return; + +} + + +/****************************************************************************** + END + ******************************************************************************/ + + +#undef DEBUG_CON2PRIMM |