Shrirang Abhyankar avatar Shrirang Abhyankar committed 8a3e2e9

First dig at solving DAEs in semi-explicit form.

Comments (0)

Files changed (12)

include/petsc-private/tsimpl.h

               TS_STEP_COMPLETE    /* step accepted and ptime, steps, etc have been advanced */
 } TSStepStatus;
 
+PETSC_INTERN PetscErrorCode TSDestroy_DAESimple(TS);
+PETSC_INTERN PetscErrorCode TSReset_DAESimple(TS);
+PETSC_INTERN PetscErrorCode TSSetFromOptions_DAESimple(TS);
+PETSC_INTERN PetscErrorCode TSSetUp_DAESimple(TS);
+PETSC_INTERN PetscErrorCode TSSolve_DAESimple(TS);
+
 #endif

include/petscts.h

 #define TSARKIMEX         "arkimex"
 #define TSROSW            "rosw"
 #define TSEIMEX           "eimex"
+#define TSDAESIMPLE       "daesimple"
+
 /*E
     TSProblemType - Determines the type of problem this TS object is to be used to solve
 
 PETSC_EXTERN PetscErrorCode TSGetTimeStepNumber(TS,PetscInt*);
 PETSC_EXTERN PetscErrorCode TSSetTimeStep(TS,PetscReal);
 
+PETSC_EXTERN PetscErrorCode TSDAESimpleSetRHSFunction(TS,Vec,PetscErrorCode (*)(PetscReal,Vec,Vec,Vec,void*),void*);
+PETSC_EXTERN PetscErrorCode TSDAESimpleSetIFunction(TS,Vec,PetscErrorCode (*)(PetscReal,Vec,Vec,Vec,void*),void*);
+
 PETSC_EXTERN_TYPEDEF typedef PetscErrorCode (*TSRHSFunction)(TS,PetscReal,Vec,Vec,void*);
 PETSC_EXTERN_TYPEDEF typedef PetscErrorCode (*TSRHSJacobian)(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);
 PETSC_EXTERN PetscErrorCode TSSetRHSFunction(TS,Vec,TSRHSFunction,void*);

src/ts/examples/tests/ex10.c

 }
 
 #undef __FUNCT__
-#define __FUNCT__ "TSDAESimpleSetRHSFunction"
-PetscErrorCode TSDAESimpleSetRHSFunction(TSDAESimple tsdae,Vec U,PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
+#define __FUNCT__ "TS_DAESimpleSetRHSFunction"
+PetscErrorCode TS_DAESimpleSetRHSFunction(TSDAESimple tsdae,Vec U,PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
 {
   PetscErrorCode ierr;
 
 }
 
 #undef __FUNCT__
-#define __FUNCT__ "TSDAESimpleSetIFunction"
-PetscErrorCode TSDAESimpleSetIFunction(TSDAESimple tsdae,Vec V,PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
+#define __FUNCT__ "TS_DAESimpleSetIFunction"
+PetscErrorCode TS_DAESimpleSetIFunction(TSDAESimple tsdae,Vec V,PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
 {
   PetscErrorCode ierr;
 
 
   ierr = VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);CHKERRQ(ierr);
   ierr = VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&V);CHKERRQ(ierr);
-  ierr = TSDAESimpleSetRHSFunction(tsdae,U,f,NULL);CHKERRQ(ierr);
-  ierr = TSDAESimpleSetIFunction(tsdae,V,F,NULL);CHKERRQ(ierr);
+  ierr = TS_DAESimpleSetRHSFunction(tsdae,U,f,NULL);CHKERRQ(ierr);
+  ierr = TS_DAESimpleSetIFunction(tsdae,V,F,NULL);CHKERRQ(ierr);
 
   ierr = VecDuplicate(U,&Usolution);CHKERRQ(ierr);
   ierr = VecSet(Usolution,1.0);CHKERRQ(ierr);

src/ts/examples/tests/ex11.c

+static char help[] = "Simple wrapper object to solve DAE of the form:\n\
+                             \\dot{U} = f(U,V)\n\
+                             F(U,V) = 0\n\n";
+
+#include <petscts.h>
+
+#undef __FUNCT__
+#define __FUNCT__ "f"
+/*
+   Simple example:   f(U,V) = U + V
+
+*/
+PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = VecWAXPY(F,1.0,U,V);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "F"
+/*
+   Simple example: F(U,V) = U - V
+
+*/
+PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = VecWAXPY(F,-1.0,V,U);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc,char **argv)
+{
+  PetscErrorCode ierr;
+  TS             ts;
+  Vec            U,V,Usolution;
+
+  PetscInitialize(&argc,&argv,(char*)0,help);
+  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
+  ierr = TSSetType(ts,TSDAESIMPLE);CHKERRQ(ierr);
+  ierr = VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);CHKERRQ(ierr);
+  ierr = VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&V);CHKERRQ(ierr);
+
+  ierr = VecDuplicate(U,&Usolution);CHKERRQ(ierr);
+  ierr = VecSet(Usolution,1.0);CHKERRQ(ierr);
+
+  ierr = TSDAESimpleSetRHSFunction(ts,Usolution,f,NULL);CHKERRQ(ierr);
+  ierr = TSDAESimpleSetIFunction(ts,V,F,NULL);CHKERRQ(ierr);
+
+  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
+  ierr = TSSolve(ts,Usolution);CHKERRQ(ierr);
+  ierr = TSDestroy(&ts);CHKERRQ(ierr);
+
+  ierr = VecDestroy(&U);CHKERRQ(ierr);
+  ierr = VecDestroy(&Usolution);CHKERRQ(ierr);
+  ierr = VecDestroy(&V);CHKERRQ(ierr);
+  PetscFinalize();
+  return 0;
+}
+
+
+

src/ts/examples/tests/makefile

 CPPFLAGS        =
 FPPFLAGS        =
 LOCDIR          = src/ts/examples/tests/
-EXAMPLESC       = ex2.c ex3.c ex4.c ex5.c ex6.c ex7.c ex8.c ex9.c ex10.c ex25.c
+EXAMPLESC       = ex2.c ex3.c ex4.c ex5.c ex6.c ex7.c ex8.c ex9.c ex10.c ex11.c ex25.c
 EXAMPLESF       =
 EXAMPLESFH      =
 MANSEC          = TS
 	-${CLINKER} -o ex10 ex10.o ${PETSC_TS_LIB}
 	${RM} ex10.o
 
+ex11: ex11.o  chkopts
+	-${CLINKER} -o ex11 ex11.o ${PETSC_TS_LIB}
+	${RM} ex11.o
+
 ex25: ex25.o  chkopts
 	-${CLINKER} -o ex25 ex25.o ${PETSC_TS_LIB}
 	${RM} ex25.o

src/ts/impls/makefile

 
 ALL: lib
 
-DIRS     = explicit implicit pseudo python arkimex rosw eimex
+DIRS     = explicit implicit pseudo python arkimex rosw eimex semiexplicit
 LOCDIR   = src/ts/impls/
 MANSEC   = TS
 

src/ts/impls/semiexplicit/makefile

+
+ALL: lib
+
+CFLAGS   =
+FFLAGS   =
+SOURCEC  = semiexplicit.c
+SOURCEF  =
+SOURCEH  =
+LIBBASE  = libpetscts
+MANSEC   = TS
+DIRS     = reduced
+LOCDIR   = src/ts/impls/semiexplicit/
+
+include ${PETSC_DIR}/conf/variables
+include ${PETSC_DIR}/conf/rules
+include ${PETSC_DIR}/conf/test
+
+

src/ts/impls/semiexplicit/reduced/makefile

+
+ALL: lib
+
+CFLAGS   =
+FFLAGS   =
+SOURCEC  = reduced.c
+SOURCEF  =
+SOURCEH  =
+LIBBASE  = libpetscts
+MANSEC   = TS
+DIRS     = 
+LOCDIR   = src/ts/impls/semiexplicit/reduced
+
+include ${PETSC_DIR}/conf/variables
+include ${PETSC_DIR}/conf/rules
+include ${PETSC_DIR}/conf/test
+
+

src/ts/impls/semiexplicit/reduced/reduced.c

+
+#include <../src/ts/impls/semiexplicit/semiexplicit.h>
+
+typedef struct {
+  PetscReal t;
+  TS        ts;
+  SNES      snes;
+  Vec       U;
+} TS_DAESimple_Reduced;
+
+#undef __FUNCT__
+#define __FUNCT__ "TSReset_DAESimple_Reduced"
+PetscErrorCode TSReset_DAESimple_Reduced(TS ts)
+{
+  TS_DAESimple *dae=(TS_DAESimple*)ts->data;
+  TS_DAESimple_Reduced *red = (TS_DAESimple_Reduced*)dae->data;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = TSDestroy(&red->ts);CHKERRQ(ierr);
+  ierr = SNESDestroy(&red->snes);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSDestroy_DAESimple_Reduced"
+PetscErrorCode TSDestroy_DAESimple_Reduced(TS ts)
+{
+  TS_DAESimple    *dae = (TS_DAESimple*)ts->data;
+  TS_DAESimple_Reduced *red = (TS_DAESimple_Reduced*)dae->data;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = TSReset_DAESimple_Reduced(ts);CHKERRQ(ierr);
+  ierr = PetscFree(red);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSSetFromOptions_DAESimple_Reduced"
+PetscErrorCode TSSetFromOptions_DAESimple_Reduced(TS ts)
+{
+
+  PetscFunctionBegin;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSSolve_DAESimple_Reduced"
+PetscErrorCode TSSolve_DAESimple_Reduced(TS ts)
+{
+  PetscErrorCode      ierr;
+  TS_DAESimple         *tsdae = (TS_DAESimple*)ts->data;
+  TS_DAESimple_Reduced *red = (TS_DAESimple_Reduced*)tsdae->data;
+
+  PetscFunctionBegin;
+  ierr = TSSolve(red->ts,tsdae->U);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSDAESimple_Reduced_TSFunction"
+/*
+   Defines the RHS function that is passed to the time-integrator.
+
+   Solves F(U,V) for V and then computes f(U,V)
+
+*/
+PetscErrorCode TSDAESimple_Reduced_TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
+{
+  TS_DAESimple         *tsdae = (TS_DAESimple*)actx;
+  TS_DAESimple_Reduced *red = (TS_DAESimple_Reduced*)tsdae->data;
+  PetscErrorCode       ierr;
+
+  PetscFunctionBeginUser;
+  red->t = t;
+  red->U = U;
+  ierr   = SNESSolve(red->snes,NULL,tsdae->V);CHKERRQ(ierr);
+  ierr   = (*tsdae->f)(t,U,tsdae->V,F,tsdae->fctx);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSDAESimple_Reduced_SNESFunction"
+/*
+   Defines the nonlinear function that is passed to the nonlinear solver
+
+*/
+PetscErrorCode TSDAESimple_Reduced_SNESFunction(SNES snes,Vec V,Vec F,void *actx)
+{
+  TS_DAESimple         *tsdae = (TS_DAESimple*)actx;
+  TS_DAESimple_Reduced *red = (TS_DAESimple_Reduced*)tsdae->data;
+  PetscErrorCode       ierr;
+
+  PetscFunctionBeginUser;
+  ierr = (*tsdae->F)(red->t,red->U,V,F,tsdae->Fctx);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSSetUp_DAESimple_Reduced"
+PetscErrorCode TSSetUp_DAESimple_Reduced(TS ts)
+{
+  PetscErrorCode       ierr;
+  TS_DAESimple         *tsdae=(TS_DAESimple*)ts->data;
+  TS_DAESimple_Reduced *red = (TS_DAESimple_Reduced*)tsdae->data;
+  Vec                  tsrhs;
+
+  PetscFunctionBegin;
+
+  ierr = TSCreate(PetscObjectComm((PetscObject)ts),&red->ts);CHKERRQ(ierr);
+  ierr = TSSetProblemType(red->ts,TS_NONLINEAR);CHKERRQ(ierr);
+  ierr = TSSetType(red->ts,TSEULER);CHKERRQ(ierr);
+  ierr = VecDuplicate(tsdae->U,&tsrhs);CHKERRQ(ierr);
+  ierr = TSSetRHSFunction(red->ts,tsrhs,TSDAESimple_Reduced_TSFunction,tsdae);CHKERRQ(ierr);
+  ierr = TSSetFromOptions(red->ts);CHKERRQ(ierr);
+  ierr = VecDestroy(&tsrhs);CHKERRQ(ierr);
+
+  ierr = SNESCreate(PetscObjectComm((PetscObject)ts),&red->snes);CHKERRQ(ierr);
+  ierr = SNESSetOptionsPrefix(red->snes,"tsdaesimple_");CHKERRQ(ierr);
+  ierr = SNESSetFunction(red->snes,NULL,TSDAESimple_Reduced_SNESFunction,tsdae);CHKERRQ(ierr);
+  ierr = SNESSetJacobian(red->snes,NULL,NULL,SNESComputeJacobianDefault,tsdae);CHKERRQ(ierr);
+  ierr = SNESSetFromOptions(red->snes);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+
+/* ------------------------------------------------------------ */
+/*MC
+      TSDAESimple - Semi-explicit DAE solver
+
+   Level: advanced
+
+.seealso:  TSCreate(), TS, TSSetType(), TSCN, TSBEULER, TSThetaSetTheta(), TSThetaSetEndpoint()
+
+M*/
+#undef __FUNCT__
+#define __FUNCT__ "TSCreate_DAESimple_Reduced"
+PETSC_EXTERN PetscErrorCode TSCreate_DAESimple_Reduced(TS ts)
+{
+  TS_DAESimple    *tsdae;
+  TS_DAESimple_Reduced *red;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBegin;
+  ts->ops->reset          = TSReset_DAESimple_Reduced;
+  ts->ops->destroy        = TSDestroy_DAESimple_Reduced;
+  ts->ops->setup          = TSSetUp_DAESimple_Reduced;
+  ts->ops->setfromoptions = TSSetFromOptions_DAESimple;
+  ts->ops->solve          = TSSolve_DAESimple;
+
+  ierr = PetscNewLog(ts,TS_DAESimple,&tsdae);CHKERRQ(ierr);
+  ts->data = (void*)tsdae;
+
+  tsdae->setfromoptions = TSSetFromOptions_DAESimple_Reduced;
+  tsdae->solve          = TSSolve_DAESimple_Reduced;
+  tsdae->destroy        = TSDestroy_DAESimple_Reduced;
+
+  ierr = PetscMalloc(sizeof(TS_DAESimple_Reduced),&red);CHKERRQ(ierr);
+  tsdae->data = red;
+
+  PetscFunctionReturn(0);
+}

src/ts/impls/semiexplicit/semiexplicit.c

+
+#include <../src/ts/impls/semiexplicit/semiexplicit.h>            /*I  "petscts.h"  */
+
+
+
+#undef __FUNCT__
+#define __FUNCT__ "TSReset_DAESimple"
+PetscErrorCode TSReset_DAESimple(TS ts)
+{
+  TS_DAESimple *dae=(TS_DAESimple*)ts->data;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = (*dae->destroy)(ts);CHKERRQ(ierr);
+  ierr = VecDestroy(&(dae)->U);CHKERRQ(ierr);
+  ierr = VecDestroy(&(dae)->V);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSDestroy_DAESimple"
+PetscErrorCode TSDestroy_DAESimple(TS ts)
+{
+
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = TSReset_DAESimple(ts);CHKERRQ(ierr);
+  ierr = PetscFree(ts->data);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSSetFromOptions_DAESimple"
+PetscErrorCode TSSetFromOptions_DAESimple(TS ts)
+{
+  TS_DAESimple   *dae=(TS_DAESimple*)ts->data;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = (*dae->setfromoptions)(ts);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSDAESimpleSetRHSFunction"
+PetscErrorCode TSDAESimpleSetRHSFunction(TS ts,Vec U,PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
+{
+  TS_DAESimple *tsdae = (TS_DAESimple*)ts->data;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  tsdae->f    = f;
+  tsdae->U    = U;
+  ierr        = PetscObjectReference((PetscObject)U);CHKERRQ(ierr);
+  tsdae->fctx = ctx;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSDAESimpleSetIFunction"
+PetscErrorCode TSDAESimpleSetIFunction(TS ts,Vec V,PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
+{
+  TS_DAESimple *tsdae = (TS_DAESimple*)ts->data;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  tsdae->F    = F;
+  tsdae->V    = V;
+  ierr        = PetscObjectReference((PetscObject)V);CHKERRQ(ierr);
+  tsdae->Fctx = ctx;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSSolve_DAESimple"
+PetscErrorCode TSSolve_DAESimple(TS ts)
+{
+  TS_DAESimple   *tsdae=(TS_DAESimple*)ts->data;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = (*tsdae->solve)(ts);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}

src/ts/impls/semiexplicit/semiexplicit.h

+/*
+  Code for solving DAEs in semi-explicit form
+*/
+#include <petsc-private/tsimpl.h>
+
+typedef struct {
+  Vec            U,V;
+  PetscErrorCode (*setfromoptions)(TS);
+  PetscErrorCode (*solve)(TS);
+  PetscErrorCode (*destroy)(TS);
+  PetscErrorCode (*view)(TS,PetscViewer);
+  PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec,void*);
+  PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec,void*);
+  void           *fctx,*Fctx;
+  void           *data;
+}TS_DAESimple;
+
+extern PetscErrorCode TSDestroy_DAESimple(TS);
+extern PetscErrorCode TSReset_DAESimple(TS);
+extern PetscErrorCode TSSetFromOptions_DAESimple(TS);
+extern PetscErrorCode TSSetUp_DAESimple(TS);
+extern PetscErrorCode TSSolve_DAESimple(TS);

src/ts/interface/tsregall.c

 PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS);
 PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS);
 PETSC_EXTERN PetscErrorCode TSCreate_EIMEX(TS);
+PETSC_EXTERN PetscErrorCode TSCreate_DAESimple_Reduced(TS);
 
 #undef __FUNCT__
 #define __FUNCT__ "TSRegisterAll"
   ierr = TSRegister(TSARKIMEX,  TSCreate_ARKIMEX);CHKERRQ(ierr);
   ierr = TSRegister(TSROSW,     TSCreate_RosW);CHKERRQ(ierr);
   ierr = TSRegister(TSEIMEX,    TSCreate_EIMEX);CHKERRQ(ierr);
+  ierr = TSRegister(TSDAESIMPLE, TSCreate_DAESimple_Reduced);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.