Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

solve_pde function for PDE problems #931

Closed
yizhang-yiz opened this issue Jul 10, 2018 · 36 comments
Closed

solve_pde function for PDE problems #931

yizhang-yiz opened this issue Jul 10, 2018 · 36 comments
Labels
feature
Milestone
2.19.0

Comments

@yizhang-yiz
Copy link
Contributor

yizhang-yiz commented Jul 10, 2018 •
edited
Loading

Summary:

Use solve_pde to solve PDEs.

Description:

The design involves cmdstan, stan, and math. See also stan-dev/stan#2567. This issue is about Math repo only.

In order to solve PDE, similar to ODE solver in Stan, the PDE is provided through a functor. Unlike ODE solver, we cannot ship a PDE library with Stan Math. Alternatively, the user is asked to provide his own PDE solver. His workflow could be the following.

  • provide a functor that uses external PDE solver. The functor has the following signature
  inline std::vector<std::vector<double> >
  operator()(const std::vector<double>& theta,
             const bool need_sens,
             const std::vector<double>& x_r,
             const std::vector<int>& x_i,
             std::ostream* msgs = nullptr);

The functor returns PDE solution(quantity of interest) when need_sens=0, and returns PDE solution together with its gradient with respect to theta(parameters) when need_sens=1.

  • call solve_pde, with the above functor as the first argument. Within the solve_pde function, the above functor calculates PDE solution and sensitivity, andvar variables are constructed based on these values.
  • provide Makefile for the external PDE solver. See separate issue Making mechanism for external library/user-supplied code #934 .

The user/PDE developer does not need to know utilize stan::math::var to write the above functor. He just needs to decide the behavior of the PDE solver: how the solution and the sensitivity are calculated. solve_pde should be able to take these values and construct proper var items, possibly using precomputed_gradient.

Additional Information:

Provide any additional information here.

Current Version:

v2.17.0
@yizhang-yiz yizhang-yiz mentioned this issue Jul 10, 2018
Closed
3 tasks
@syclik
Copy link
Member

syclik commented Jul 10, 2018 via email

@yizhang-yiz
Copy link
Contributor Author

Will there be a solve_pde function in the math library?

Yes.

What does that implementation need to do?

For prim version, just use the PDE functor to solve PDE. For rev version, use the PDE functor to solve PDE and sensitivity, then construct var using precomputed_gradient.

Since it sounds like it needs to be linked against something to solve, how do we provide a solver ?

We don't.

Or is the user supposed to provide a solver?

Yes.

Is there a reason we don’t handle the two different behaviors with two different functions?

I'm not sure I understand the question. If by "two different behaviors" you mean solve /w and w/o sensitivity, then yes, we can certainly do two versions of functions, something like solve_pde and solve_pde_with_sensitivity. But that would require later on user know which one to use.

@yizhang-yiz
Copy link
Contributor Author

Within Math library the solution would be similar to integrate_ode. The one proposed in #930 takes a PDE functor argument and use it to calculate PDE solutions. The difference between the solve_pde proposed in #930 and integrate_ode is that

  • The PDE functor passed to solve_pde is actually a PDE solution functor, while the ODE functor passed to integrate_ode simply calculates the ODE's right-hand-side. A PDE solution functor operates on parameters and returns PDE solutions(and sensitivities). This way in Math library we don't have to include a PDE solver.
  • In the above sense the proposed solve_pde function serves as Math libraries' API for external PDE solver developers. Users are asked to provide a functor that utilizes certain PDE solver, as well as makefile for that solver.

@yizhang-yiz yizhang-yiz mentioned this issue Jul 11, 2018
Open
@bob-carpenter
Copy link
Contributor

bob-carpenter commented Jul 11, 2018 via email

@yizhang-yiz
Copy link
Contributor Author

yizhang-yiz commented Jul 11, 2018 via email

@syclik
Copy link
Member

syclik commented Jul 11, 2018

I'm still quite a bit confused. The user is supposed to provide a functor like this:

inline std::vector<std::vector<double>> operator()(const std::vector<double>& theta,
             const bool need_sens,
             const std::vector<double>& x_r,
             const std::vector<int>& x_i,
             std::ostream* msgs = nullptr);

And passes that into solve_pde()?

That makes sense, but then why do we need any changes to the build process? Is that necessary?

@yizhang-yiz
Copy link
Contributor Author

yizhang-yiz commented Jul 11, 2018 via email

@syclik
Copy link
Member

syclik commented Jul 11, 2018

Shouldn't the user just put in make/local:

CXX += -I <path to source>
LDFLAGS += -l<library>

?

I don't understand why we need to define new variables. We did that for MPI and GPU, but I don't want to get in the habit of making our builds more complicated without there being a very good reason.

@yizhang-yiz
Copy link
Contributor Author

In the most simple case, yes. But there are more complicated scenarios, especially for HPC users, they may want to adjust optimizations etc for tuning purpose, or may want to turn on/off certain functionality in his solver. Note that we're not asking user to replace Stan Math makefile and dive into our make process, we are asking him to provide only switches that are related to his solver. In the minimalist case, the makefile could only contain what you mentioned in the local. Of course another option is just let them use stuff local. The two approaches are the same after make's include does its job.

@syclik
Copy link
Member

syclik commented Jul 11, 2018

Can you provide a concrete example? I don't quite see what the complication is yet.

@yizhang-yiz
Copy link
Contributor Author

Here's small of part of my petsc makefile rules. In general it goes through many decisions, such as if user needs to rebuild the library, or download another library, or if user is building fortran vs build c, if he's using another compiler, etc. Of course we can just ask user to include this file in local. Either way works for me.

# 3. Check if the shared libs are out of date
chkopts: chk_upgrade
	@for LIBNAME in ${SHLIBS}; do \
	  library=${INSTALL_LIB_DIR}/$$LIBNAME.a; \
	sharedlibrary=${INSTALL_LIB_DIR}/$$LIBNAME.${SL_LINKER_SUFFIX}; \
	flag=""; \
	if [ -f $$library ]; then \
	  if [ -f $$sharedlibrary ]; then \
	    flag=`find ${INSTALL_LIB_DIR} -type f -name $$LIBNAME.a -newer $$sharedlibrary -print`; \
	  fi; \
	fi; \
	if [ "$$flag" != "" ]; then \
	  echo "Shared libs in ${INSTALL_LIB_DIR} are out of date, attempting to rebuild."; \
	  if [ -w ${INSTALL_LIB_DIR} ]; then \
	    ${OMAKE}  PETSC_ARCH=${PETSC_ARCH} shared; \
	  else \
	    printf ${PETSC_TEXT_HILIGHT}"*********************** ERROR ************************\n"; \
	    echo "Unable to rebuild shared libraries; you do not have write permission."; \
	    user=`ls -l ${INSTALL_LIB_DIR}/$$LIBNAME.${SL_LINKER_SUFFIX}  | tr -s ' ' | cut -d ' ' -f 3`; \
	    echo "Libraries were built by user $$user; please contact him/her to have them rebuilt."; \
	    printf "******************************************************"${PETSC_TEXT_NORMAL}"\n" ; \
	    false; \
	  fi; \
	fi; \
	done

#
# uses the cmake infrastructure to build/rebuild the libraries
ccmake:
	@echo "=========================================="
	+@cd ${PETSC_DIR}/${PETSC_ARCH} && MAKEFLAGS="-j$(MAKE_NP) -l$(NPMAX) $(MAKEFLAGS)" ${OMAKE} VERBOSE=${VERBOSE}
	@if [ "${BUILDSHAREDLIB}" = "yes" -a "${DSYMUTIL}" != "true" ]; then \
        echo "Running ${DSYMUTIL} on ${SHLIBS}";\
        for LIBNAME in ${SHLIBS}; do ${DSYMUTIL} ${INSTALL_LIB_DIR}/$$LIBNAME.${SL_LINKER_SUFFIX}; done; fi
	@echo "========================================="

cmake:
	+@MAKEFLAGS="-j$(MAKE_NP) -l$(NPMAX) $(MAKEFLAGS)" ${OMAKE} ccmake VERBOSE=1

gnumake:
	+@cd ${PETSC_DIR} && MAKEFLAGS="-j$(MAKE_NP) -l$(NPMAX) $(MAKEFLAGS)" ${OMAKE_PRINTDIR} -f gmakefile ${MAKE_PAR_OUT_FLG} V=${V}

# Does nothing; needed for some rules that require actions.
foo:

# Builds library
lib:  ${SOURCE}
	@${OMAKE}  PETSC_ARCH=${PETSC_ARCH} chk_petscdir
	@-if [ "${SPECIALLIB}" = "yes" ] ; then \
	   ${OMAKE}  PETSC_ARCH=${PETSC_ARCH}  speciallib;  \
	  else \
            if [ "${SOURCECU}" != "" ] ; then \
              ${OMAKE}  PETSC_ARCH=${PETSC_ARCH}  libcu; fi ; \
            if [ "${SOURCEC}" != "" ] ; then \
              ${OMAKE}  PETSC_ARCH=${PETSC_ARCH}  libc; fi ; \
            if [ "${SOURCECXX}" != "" ] ; then \
              ${OMAKE}  PETSC_ARCH=${PETSC_ARCH}  libcxx; fi ; \
            if [ "${SOURCEF}" != "" ] ; then \
		${OMAKE}  PETSC_ARCH=${PETSC_ARCH}  libf; fi ; \
            if [ "${OBJS}" != " " ] ; then \
		${RANLIB}  ${LIBNAME}; \
		${RM} ${OBJS}; \
	    fi;\
          fi
#
#  Does not work for some machines with .F fortran files.
#
# Builds library - fast versiong

libfast:  ${SOURCE}
	-@if [ "${SPECIALFASTLIB}" = "yes" ] ; then \
	    ${OMAKE}  PETSC_ARCH=${PETSC_ARCH}  specialfastlib;  \
	else \
	    ${OMAKE}  PETSC_ARCH=${PETSC_ARCH}  libfastcompile;  \
	fi

libfastcompile:
	-@if [ "${SOURCECU}" != "" ]; then \
                    ${PETSC_CUCOMPILE}; \
	            ${AR} ${FAST_AR_FLAGS} ${LIBNAME} ${OBJSCU}; \
		    ${RM} ${OBJSCU} ${OBJSCU:.o=.lo}; \
	    fi; \
          if [ "${SOURCEC}" != "" ]; then \
                    ${PETSC_COMPILE}; \
	            ${AR} ${FAST_AR_FLAGS} ${LIBNAME} ${OBJSC}; \
		    ${RM} ${OBJSC} ${OBJSC:.o=.lo}; \
	    fi; \
          if [ "${SOURCECXX}" != "" ]; then \
                    ${PETSC_CXXCOMPILE}; \
	            ${AR} ${FAST_AR_FLAGS} ${LIBNAME} ${OBJSCXX}; \
		    ${RM} ${OBJSCXX} ${OBJSCXX:.o=.lo}; \
	    fi; \
	if [ "${SOURCEF}" != "" ]; then \
                    ${PETSC_FCOMPILE}; \
	            ${AR} ${FAST_AR_FLAGS} ${LIBNAME} ${OBJSF}; \
		    ${RM} ${OBJSF}; \
	    fi

# Build f90 .mod from .F source, and add .o to the corresponding library
buildmod:
	-@${OMAKE} clean-legacy
	@${OMAKE}  libf
	@${OMAKE}  modcopy
	@${OMAKE}  clean-legacy

buildmodfast:
	-@${OMAKE} clean-legacy
	@${OMAKE}  libfastcompile
	@${OMAKE}  modcopy
	@${OMAKE}  clean-legacy

# copy modules to the include dir
modcopy:
	@${CP} -f *.mod ${PETSC_DIR}/${PETSC_ARCH}/include

@syclik
Copy link
Member

syclik commented Jul 11, 2018 via email

@yizhang-yiz
Copy link
Contributor Author

yizhang-yiz commented Jul 11, 2018 •
edited
Loading

I think the bottom line is that user are looking at different optimization options, and different architectures that he may want to play with his PDE solver. Also his source code might be C or C++ or fortran or F90, all asking for special treatment.

It does add more complexity either way, IMO. Because one can either bloat up local or include the external makefile. But I think it's conceptually easier to understand to use a designated makefile for external library.

@syclik
Copy link
Member

syclik commented Jul 11, 2018 via email

@yizhang-yiz
Copy link
Contributor Author

yizhang-yiz commented Jul 11, 2018 •
edited
Loading

Sorry for being a bad interpreter.

https://github.com/yizhang-cae/cmdstan/tree/forward_pde/examples/laplace_pde
Here's the example I mentioned in stan-dev/stan#2567.
The makefile for libmesh used in this example requires the follow minimum setting:

CXX = $(libmesh_CXX)
CC = $(libmesh_CXX)
CPPFLAGS += $(libmesh_CPPFLAGS)
CXXFLAGS += $(libmesh_CXXFLAGS)
CXXFLAGS += -isystem $(LIBMESH_DIR)/include
CXXFLAGS += $(libmesh_INCLUDE)
CXXFLAGS :=$(filter-out -std=gnu++11,$(CXXFLAGS))
CXXFLAGS :=$(filter-out -Wunused,$(CXXFLAGS))
CXXFLAGS :=$(filter-out -Wunused-parameter,$(CXXFLAGS))
CXXFLAGS :=$(filter-out -Qunused-arguments,$(CXXFLAGS))
CXXFLAGS +=-DLIBMESH_HAVE_CXX14_MAKE_UNIQUE

LDFLAGS += $(libmesh_LIBS)
LDFLAGS += $(libmesh_LDFLAGS)
LDFLAGS += $(EXTERNAL_FLAGS)

# Useful rules.
dust:
	@echo "Deleting old output and runtime files"
	@rm -f out*.m job_output.txt output.txt* *.gmv.* *.plt.* out*.xdr* out*.xda* PI* *.o *.libs

gmv:
	@$(MAKE) -C $(LIBMESH_DIR)/roy/meshplot/ meshplot-$(METHOD)
	@for file in out.mesh.*; do ${LIBMESH_RUN} $(LIBMESH_DIR)/roy/meshplot/meshplot-$(METHOD) $$file out.soln.$${file##out.mesh.} out.gmv.$${file:9:4}; done

What I meant is

  • optimization options

Questions like use -O3 vs -O2? Link withmetis(a partitioner) or scotch(another partitioner)? etc.

  • different architectures

Questions like Build with MPICH or OPENMPI? linear solver being SUPERLU or MUMPS? These are related to how a solver library is built. Some libraries (like petsc) allow the user to build different versions with different setups, aka "architectures".

  • play with his PDE solver

Large problems often requires some tuning of the aforementioned building parameters before sent out to the cluster for days-long run. What I meant by "play" is this tuning process, which involves change switches in the building process.

  • special treatment

One example is that if the solver is written in fortran, user needs to supply mpif or gfortran compiler/linker information.

@bbbales2
Copy link
Member

There's nothing here that's PDE specific, so is the goal to interface with external software in general?

It seems like there are three big components of all these pull reqs:

  1. The high level interface (that would appear in .stan files)
  2. The mechanism for working with the autodiff library without having to worry too much about the details of vars
  3. Whatever is needed to link in external libraries

3 seems like a support nightmare to me, and 1 seems easy to get around (if someone can build petsc, I'll hazard they can figure out how to get the Rstan/cmdstan external C++ stuff working, even if those are hacky solutions). Seems to me like if someone has some big fancy weird external solver they need to incorporate in Stan they could maintain it as a fork? That's what I do for my regular job -- not that my development patterns are anything to emulate...

I'm working with Bob now on 2, for other reasons, but it's the same goal: #924 . I'll be full time on this until September. We think it's a good idea too -- and someone else coming up with it independently seems like a good sign :P. Fingers crossed that should progress quickly.

Maybe this should be an item for the engineering part of the Thursday meeting?

@yizhang-yiz
Copy link
Contributor Author

As I mentioned in stan-dev/stan#2567, this is indeed true to any external library. for 3 I don't think asking PDE developer to maintain a fork is a better option(I'm doing that with Torsten, and I'm desperately looking for alternatives) then providing a good API.

@yizhang-yiz
Copy link
Contributor Author

yizhang-yiz commented Jul 11, 2018 •
edited
Loading

Add @bob-carpenter for Thursday agenda: mechanism to include external library.

@yizhang-yiz
Copy link
Contributor Author

A related issue is stan-dev/stan#2572: to have a way to support user-supplied higher-order function, like regular function through user_header.hpp. Though this is more a Stan issue than Math issue.

@syclik
Copy link
Member

syclik commented Jul 11, 2018

As I mentioned in stan-dev/stan#2567, this is indeed true to any external library. for 3 I don't think asking PDE developer to maintain a fork is a better option(I'm doing that with Torsten, and I'm desperately looking for alternatives) then providing a good API.

Re: Torsten. We should put this on another thread. There was an active decision to keep that separate due to lack of resources for design, documentation, testing, and maintenance. We can revisit that if you think that's appropriate.

@yizhang-yiz
Copy link
Contributor Author

yizhang-yiz commented Jul 11, 2018 •
edited
Loading

I didn't mean to bring Torsten into Stan or involve Torsten in the current design, what I meant is I'd prefer a way to use Stan with external libraries through API(which is the topic of this thread) to include Stan as a fork. This is to respond to @bbbales2 's comment. Sorry for the confusion.

@syclik
Copy link
Member

syclik commented Jul 11, 2018

@yizhang-cae, can we split the issue into two? I think there are two very different things going on:

  1. solve_pde(). That looks fine and easy to review and merge. There's nothing different about that.
  2. Building external tools and linking them. That's a completely separate, orthogonal discussion.

For 1, we can move forward. For 2, it would really help to have a list of things you want (a spec would be great, but even a list of wants and needs would be sufficient), then we can discuss whether that's a good thing and how we can make that happen.

@yizhang-yiz yizhang-yiz mentioned this issue Jul 11, 2018
Closed
@yizhang-yiz
Copy link
Contributor Author

I've edited the issue and added #934 . Let's focus on solve_pde in the current issue.

@yizhang-yiz
Copy link
Contributor Author

yizhang-yiz commented Jul 11, 2018 •
edited
Loading

@syclik , for the current issue, please allow me to pitch the name forward_pde again. I understand solve_pde is a more popular name, but it's also a misleading name in the context of statistical inference: solve_pde gives an impression that the return would be the solution of the PDE, but that's not the intended use case. In PDE inference the data are often not collected from PDE solution directly but some functionals based on the PDE solution, usually involving spatial and/or time integrals. This is the "forward" process from PDE parameters to the observable functionals. It involves solving the PDE, but the numerical solution of PDE is often hidden behind.

@bob-carpenter
Copy link
Contributor

bob-carpenter commented Jul 11, 2018 via email

@syclik
Copy link
Member

syclik commented Jul 11, 2018

Btw, what you've described as "optimization options" are already available. In make/local, just set:

O = 3

or

O = 2

For the other ones, I'm not sure? These are all different use cases. Knowing what you're trying to do would really help. (It's still not clear enough for me to understand exactly what you're trying to do to help design.)

@syclik
Copy link
Member

syclik commented Jul 11, 2018

And yes, forward_pde sounds much more appropriate given your last comment.

@syclik
Copy link
Member

syclik commented Jul 11, 2018

That said Matlab calls it solvepde: https://www.mathworks.com/help/pde/pde-solvers.html.

What about simulate_pde? TensorFlow's doc talks about simulating from a PDE: https://www.tensorflow.org/tutorials/non-ml/pdes

So does iPython Cookbook: https://ipython-books.github.io/124-simulating-a-partial-differential-equation-reaction-diffusion-systems-and-turing-patterns/

And just to add more confusion, FEniCS calls it solve (they do PDE only): https://fenicsproject.org/

@yizhang-yiz
Copy link
Contributor Author

yizhang-yiz commented Jul 11, 2018 •
edited
Loading

That's exactly what I'm talking about.

Matlab:

This solver returns a StationaryResults or TimeDependentResults object whose properties contain the solution and its gradient at the mesh nodes.

FEniCS

Compute solution
w = Function(W)
solve(a == L, w, [bc1, bc0])

TensorFlow

Discretized PDE update rules
U_ = U + eps * Ut
Ut_ = Ut + eps * (laplace(U) - damping * Ut)

All the above examples solve for the numerical solution. None of the above example ask for QoI(quantity of interest).

BTW, the adjoint QoI solver is treated separately in FEniCS by DOLFIN, so the example above is a misquote.

@syclik
Copy link
Member

syclik commented Jul 11, 2018

I missed the point about the QoI. In the words / notation from Matlab, FEniCS, and TensorFlow, what's different about what you're proposing the Math library handles?

@syclik
Copy link
Member

syclik commented Jul 11, 2018

btw, would this go a lot quicker if we jumped on a video call?

@yizhang-yiz
Copy link
Contributor Author

Sure. Hangout? yiz@metrumrg.com

@yizhang-yiz
Copy link
Contributor Author

In the words / notation from Matlab, FEniCS, and TensorFlow, what's different about what you're proposing the Math library handles?

Don't remember how Matlab & tensorflow treat it, but for FEniCS the adjoint sensitivity is treated in a sublibrary http://www.dolfin-adjoint.org/. What we are asking user to provide in the PDE functor is the combination of the PDE solution like quoted above and the calculation of QoI like in http://www.dolfin-adjoint.org/. Of course since user is responsible for the behavior of the functor, nothing stops him from returning a full-blown PDE solution and stuff the cache. That's why hopefully the name forward would mean something to someone who's familiar with PDE inverse problems.

@syclik
Copy link
Member

syclik commented Jul 11, 2018 via email

@yizhang-yiz
Copy link
Contributor Author

Would it be clearer if we didn’t call if “pde” but something like “pde_qoi”? Since the functor is not a pde, I don’t think we should call it a pde?

@syclik , good point. Changes made at #930 (comment)

@yizhang-yiz yizhang-yiz mentioned this issue Aug 3, 2018
Closed
3 tasks
@yizhang-yiz
Copy link
Contributor Author

@syclik I'd like to have it merged into Torsten first so I'm closing it for the moment.

@syclik syclik added this to the 2.19.0 milestone Mar 18, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
feature
Projects
None yet
Milestone
2.19.0
Development

No branches or pull requests
4 participants
@syclik @bob-carpenter @bbbales2 @yizhang-yiz