How to choose what copula to use for a certain application? I'm using the copula package in R for modelling dependece using copulas.
1)What is the suggested course of action for choosing a copula model?
2)Should I use the function gofCopula() to assess whether the model is a good fit? My main reason for not using it is that it is very slow, and I'm not 100% sure how the hypothesis testing is formulated and therefore how to interpret the p-value (the docs don't help much). 
Right now what I'm doing is fitting the copula to the data and generate simulated observations from it, if the simulated data has proprieties which are similar to the observed data then I use the model, otherwise I don't.
This is a toy example where a normal copula looks to be appropriate:
# Generate 10000 correlated observation with cor=0.8
set.seed(5)
x <- rnorm(10000)
y <- correlatedValue(x=x, r=.8)

# Check cor(x,y) == 0.8 (approx)
cor(x,y)

# Generate Empirical copula assuming normal marginals (safe assumption here I guess)
u <- pnorm(x,0,1)
v <- pnorm(y,0,1)
plot(u,v,pch='.',col='blue',main='Empirical copula')

# Use the copula package to fit and generate the copula
library(copula)

# We use a normal copula. Copula fitting
m <- as.matrix(cbind(u,v)) # Or pobs(as.matrix(cbind(x,y)))
normal.cop <- normalCopula(dim=2)
fit.cop<- fitCopula(normal.cop,m,method="ml")

# Coefficients of the Copula
rho <- coef(fit.cop)
print(rho)

# Simulate data using the fited copula
u1 = rCopula(10000,normalCopula(coef(fit.cop),dim=2))
points(u1[,1],u1[,2],col="red",pch='.',main='Simulated copula')

# Get simulated data back into its original scale
u_sim <- u1[,1]

v_sim <- u1[,2]

# Getting back simulated x and y
x1 <- qnorm(u_sim,0,1)
y1 <- qnorm(v_sim,0,1)

# Cor should be approximately 0.8
cor(x1,y1)

# Plot original observation against simulated
plot(x,y,pch='.',col='blue')
points(x1,y1,pch='.',col='red')

# Add a legend
legend("topleft",c("Simulated","Observed"), col=c("red", "blue"), pch=20)

 A: Yes, you should use a goodness-of-fit test for model selection in conjunction with cross-validation train/test splits to avoid over-fitting and model biasing. 
In general, it is indeed true that performing these gof tests can be time-consuming. Especially if one chooses to also compare fitness of the rotated copulas, the number of candidate models can quickly go beyond 30-40 types. One suggestion I can make is computing some non-linear dependence statistics such as kendall's tau etc. which may give an indication of what kind of dependency your data is exhibiting. Using this information you can narrow down your search to only a subset of all available families.
Another way you can test is to fit all your candidate models, simulate from them and then perform a two-sample test between original and sampled data (in the oringal domain). You susequently choose the model with the best results. An overview of multivariate two-sample tests can be found here:
https://normaldeviate.wordpress.com/2012/07/14/modern-two-sample-tests/
